tetmesh_faceted_surfaces.h

Go to the documentation of this file.

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented, 
//LIC// multi-physics finite-element library, available 
//LIC// at http://www.oomph-lib.org.
//LIC// 
//LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
//LIC// 
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// License as published by the Free Software Foundation; either
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC// 
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//LIC// Lesser General Public License for more details.
//LIC// 
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301  USA.
//LIC// 
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC// 
//LIC//====================================================================
 
 
//Include guard to prevent multiple inclusions of the header
#ifndef OOMPH_TETMESH_FACETED_SURFACES_HEADER
#define OOMPH_TETMESH_FACETED_SURFACES_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
  #include <oomph-lib-config.h>
#endif
 
 
namespace oomph
{
 
//================================================================
//================================================================
namespace HelperForSortingVectorOfPairsOfVertexPtAndDouble
{
 bool less_than_based_on_double(std::pair<TetMeshVertex*,double> i,
                                std::pair<TetMeshVertex*,double> j) 
 {
  return (i.second<j.second); 
 }
 
}
 
 
 
 
 
//=============================================================================
//=============================================================================
class CubicTetMeshFacetedSurface : public virtual TetMeshFacetedClosedSurface
{
 
 
public: 
 
 CubicTetMeshFacetedSurface(const double& box_half_width, 
                            const double& box_half_length,
                            const unsigned& one_based_boundary_id_offset=0)
  {
   unsigned one_based_region_id=0;
   Vector<double> offset(3,0.0);
   build_it(one_based_region_id,
            box_half_width,
            box_half_length,
            offset,
            one_based_boundary_id_offset);
  }
   
 CubicTetMeshFacetedSurface(const double& box_half_width, 
                            const double& box_half_length,
                            const Vector<double>& offset,
                            const unsigned& one_based_boundary_id_offset=0)
  {   
   unsigned one_based_region_id=0;
   build_it(one_based_region_id,
            box_half_width,
            box_half_length,
            offset,
            one_based_boundary_id_offset);
  }
 CubicTetMeshFacetedSurface(const unsigned& one_based_region_id,
                            const double& box_half_width, 
                            const double& box_half_length,
                            const unsigned& one_based_boundary_id_offset=0)
  {
   Vector<double> offset(3,0.0);
   build_it(one_based_region_id,
            box_half_width,
            box_half_length,
            offset,
            one_based_boundary_id_offset);
  }
   
 CubicTetMeshFacetedSurface(const unsigned& one_based_region_id,
                            const double& box_half_width, 
                            const double& box_half_length,
                            const Vector<double>& offset,
                            const unsigned& one_based_boundary_id_offset=0)
  {   
   build_it(one_based_region_id,
            box_half_width,
            box_half_length,
            offset,
            one_based_boundary_id_offset);
  }
 
 
private:
 
 void build_it(const unsigned& one_based_region_id,
               const double& box_half_width, 
               const double& box_half_length,
               const Vector<double>& offset,
               const unsigned& one_based_boundary_id_offset)
  {
   // Make vertices
   unsigned n_vertex=8;
   Vertex_pt.resize(n_vertex);
   Vector<double> box_point(3);
 
   box_point[0] = offset[0]-box_half_width;
   box_point[1] = offset[1]-box_half_width;
   box_point[2] = offset[2]-box_half_length;
   Vertex_pt[0]=new TetMeshVertex(box_point);
  
   box_point[0] = offset[0]-box_half_width;
   box_point[1] = offset[1]+box_half_width;
   box_point[2] = offset[2]-box_half_length;
   Vertex_pt[1]=new TetMeshVertex(box_point);
 
   box_point[0] = offset[0]-box_half_width;
   box_point[1] = offset[1]+box_half_width;
   box_point[2] = offset[2]+box_half_length;
   Vertex_pt[2]=new TetMeshVertex(box_point);
   
   box_point[0] = offset[0]-box_half_width;
   box_point[1] = offset[1]-box_half_width;
   box_point[2] = offset[2]+box_half_length;
   Vertex_pt[3]=new TetMeshVertex(box_point);
 
   box_point[0] = offset[0]+box_half_width;
   box_point[1] = offset[1]-box_half_width;
   box_point[2] = offset[2]-box_half_length;
   Vertex_pt[4]=new TetMeshVertex(box_point);
   
   box_point[0] = offset[0]+box_half_width;
   box_point[1] = offset[1]+box_half_width;
   box_point[2] = offset[2]-box_half_length;
   Vertex_pt[5]=new TetMeshVertex(box_point);
   
   box_point[0] = offset[0]+box_half_width;
   box_point[1] = offset[1]+box_half_width;
   box_point[2] = offset[2]+box_half_length;
   Vertex_pt[6]=new TetMeshVertex(box_point);
   
   box_point[0] = offset[0]+box_half_width;
   box_point[1] = offset[1]-box_half_width;
   box_point[2] = offset[2]+box_half_length;
   Vertex_pt[7]=new TetMeshVertex(box_point);
   
 
   // Make facets
   unsigned n_facet=6;
   Facet_pt.resize(n_facet);
   
   unsigned n_vertex_on_facet=4;
   Facet_pt[0]=new TetMeshFacet(n_vertex_on_facet);
   unsigned one_based_boundary_id=1+one_based_boundary_id_offset;
   Facet_pt[0]->set_one_based_boundary_id(one_based_boundary_id);   
   Facet_pt[0]->set_vertex_pt(0,Vertex_pt[0]);
   Facet_pt[0]->set_vertex_pt(1,Vertex_pt[4]);
   Facet_pt[0]->set_vertex_pt(2,Vertex_pt[7]);
   Facet_pt[0]->set_vertex_pt(3,Vertex_pt[3]);
 
   Facet_pt[1]=new TetMeshFacet(n_vertex_on_facet);
   one_based_boundary_id=2+one_based_boundary_id_offset;
   Facet_pt[1]->set_one_based_boundary_id(one_based_boundary_id); 
   Facet_pt[1]->set_vertex_pt(0,Vertex_pt[4]);
   Facet_pt[1]->set_vertex_pt(1,Vertex_pt[5]);
   Facet_pt[1]->set_vertex_pt(2,Vertex_pt[6]);
   Facet_pt[1]->set_vertex_pt(3,Vertex_pt[7]);
 
   // top
   Facet_pt[2]=new TetMeshFacet(n_vertex_on_facet);
   one_based_boundary_id=3+one_based_boundary_id_offset;
   Facet_pt[2]->set_one_based_boundary_id(one_based_boundary_id); 
   Facet_pt[2]->set_vertex_pt(0,Vertex_pt[3]);
   Facet_pt[2]->set_vertex_pt(1,Vertex_pt[7]);
   Facet_pt[2]->set_vertex_pt(2,Vertex_pt[6]);
   Facet_pt[2]->set_vertex_pt(3,Vertex_pt[2]);
 
   Facet_pt[3]=new TetMeshFacet(n_vertex_on_facet);
   one_based_boundary_id=4+one_based_boundary_id_offset;
   Facet_pt[3]->set_one_based_boundary_id(one_based_boundary_id);   
   Facet_pt[3]->set_vertex_pt(0,Vertex_pt[6]);
   Facet_pt[3]->set_vertex_pt(1,Vertex_pt[5]);
   Facet_pt[3]->set_vertex_pt(2,Vertex_pt[1]);
   Facet_pt[3]->set_vertex_pt(3,Vertex_pt[2]);
 
   Facet_pt[4]=new TetMeshFacet(n_vertex_on_facet);
   one_based_boundary_id=5+one_based_boundary_id_offset;
   Facet_pt[4]->set_one_based_boundary_id(one_based_boundary_id);  
   Facet_pt[4]->set_vertex_pt(0,Vertex_pt[3]);
   Facet_pt[4]->set_vertex_pt(1,Vertex_pt[2]);
   Facet_pt[4]->set_vertex_pt(2,Vertex_pt[1]);
   Facet_pt[4]->set_vertex_pt(3,Vertex_pt[0]);
 
   // bottom
   Facet_pt[5]=new TetMeshFacet(n_vertex_on_facet);
   one_based_boundary_id=6+one_based_boundary_id_offset;
   Facet_pt[5]->set_one_based_boundary_id(one_based_boundary_id);   
   Facet_pt[5]->set_vertex_pt(0,Vertex_pt[5]);
   Facet_pt[5]->set_vertex_pt(1,Vertex_pt[1]);
   Facet_pt[5]->set_vertex_pt(2,Vertex_pt[0]);
   Facet_pt[5]->set_vertex_pt(3,Vertex_pt[4]);
 
 
   // Specify the adjacent region that is bounded by these facets
   for (unsigned f=0;f<6;f++)
    {
     Facet_pt[f]->set_one_based_adjacent_region_id(one_based_region_id);
    }
 
   // Is it a genuine region?
   if (one_based_region_id>0)
    {
     set_region_for_tetgen(one_based_region_id-1,offset);
    }
   // Otherwise declare it to be a hole
   else
    {
     set_hole_for_tetgen(offset);
     enable_faceted_volume_represents_hole_for_gmsh();
    }
 
  }
 
};
 
 
 
 
 
 
//=============================================================================
//=============================================================================
class RectangularTetMeshFacetedSurface : public virtual TetMeshFacetedSurface
{
 
 
public: 
 
 
 RectangularTetMeshFacetedSurface(const double& half_x_width, 
                                  const double& half_y_length,
                                  const Vector<double>& offset,
                                  const unsigned& one_based_boundary_id)
  {
   build_it(half_x_width,
            half_y_length,
            offset,
            one_based_boundary_id);
  }
   
 
 
private:
 
 void build_it(const double& half_x_width, 
               const double& half_y_length,
               const Vector<double>& offset,
               const unsigned& one_based_boundary_id)
  {
   // Make vertices
   unsigned n_vertex=4;
   Vertex_pt.resize(n_vertex);
   Vector<double> box_point(3);
 
   box_point[0] = offset[0]-half_x_width;
   box_point[1] = offset[1]-half_y_length;
   box_point[2] = offset[2];
   Vertex_pt[0]=new TetMeshVertex(box_point);
  
   box_point[0] = offset[0]+half_x_width;
   box_point[1] = offset[1]-half_y_length;
   box_point[2] = offset[2];
   Vertex_pt[1]=new TetMeshVertex(box_point);
 
   box_point[0] = offset[0]+half_x_width;
   box_point[1] = offset[1]+half_y_length;
   box_point[2] = offset[2];
   Vertex_pt[2]=new TetMeshVertex(box_point);
   
   box_point[0] = offset[0]-half_x_width;
   box_point[1] = offset[1]+half_y_length;
   box_point[2] = offset[2];
   Vertex_pt[3]=new TetMeshVertex(box_point);
 
 
   // Make facets
   unsigned n_facet=1;
   Facet_pt.resize(n_facet);
   
   unsigned n_vertex_on_facet=4;
   Facet_pt[0]=new TetMeshFacet(n_vertex_on_facet);
   Facet_pt[0]->set_one_based_boundary_id(one_based_boundary_id);   
   Facet_pt[0]->set_vertex_pt(0,Vertex_pt[0]);
   Facet_pt[0]->set_vertex_pt(1,Vertex_pt[1]);
   Facet_pt[0]->set_vertex_pt(2,Vertex_pt[2]);
   Facet_pt[0]->set_vertex_pt(3,Vertex_pt[3]);
 
  }
 
};
 
 
 
 
 
//================================================================
//================================================================
class DiskTetMeshFacetedSurface : public virtual TetMeshFacetedSurface
{
 
 
public:
 
 DiskTetMeshFacetedSurface(
  DiskLikeGeomObjectWithBoundaries* 
  disk_parametrised_by_nonsingular_coordinates_pt,
  const unsigned& half_nsegment,
  const unsigned& first_one_based_boundary_id_for_disk,
  unsigned& last_one_based_boundary_id_for_disk) 
 
  {
   Geom_object_with_boundaries_pt=
    disk_parametrised_by_nonsingular_coordinates_pt;
 
   // Provide storage for pointers to the two parts of the curvilinear boundary
   Vector<TriangleMeshCurveSection*> outer_curvilinear_boundary_pt(2);
   
   // First bit
   GeomObject* outer_boundary_ellipse0_pt=
    disk_parametrised_by_nonsingular_coordinates_pt->
    boundary_parametrising_geom_object_pt(0);
   double zeta_start=disk_parametrised_by_nonsingular_coordinates_pt->
    zeta_boundary_start(0);
   double zeta_end=disk_parametrised_by_nonsingular_coordinates_pt->
    zeta_boundary_end(0);
   unsigned nsegment=half_nsegment;
   unsigned boundary_id=0;
   outer_curvilinear_boundary_pt[0]=new TriangleMeshCurviLine(
    outer_boundary_ellipse0_pt,zeta_start,zeta_end,nsegment,boundary_id);
   
   // Second bit
   GeomObject* outer_boundary_ellipse1_pt=
    disk_parametrised_by_nonsingular_coordinates_pt->
    boundary_parametrising_geom_object_pt(1);
   zeta_start=disk_parametrised_by_nonsingular_coordinates_pt->
    zeta_boundary_start(1);
   zeta_end=disk_parametrised_by_nonsingular_coordinates_pt->
    zeta_boundary_end(1);
   nsegment=half_nsegment;
   boundary_id=1;
   outer_curvilinear_boundary_pt[1]=new TriangleMeshCurviLine(
    outer_boundary_ellipse1_pt,zeta_start,zeta_end,nsegment,boundary_id);
   
   // Combine to curvilinear boundary and define the
   // outer boundary
   TriangleMeshClosedCurve* closed_curve_pt=
    new TriangleMeshClosedCurve(outer_curvilinear_boundary_pt);
     
   // Use the TriangleMeshParameters object for helping on the manage of the
   // TriangleMesh parameters
   TriangleMeshParameters triangle_mesh_parameters(closed_curve_pt);
   
   // Do we have an annnular internal boundary?
   if (disk_parametrised_by_nonsingular_coordinates_pt->nboundary()==4)
    {
     
     // Provide storage for pointers to the two parts of the curvilinear boundary
     Vector<TriangleMeshCurveSection*> inner_curvilinear_boundary_pt(2);
   
     // First bit
     GeomObject* inner_boundary_ellipse0_pt=
      disk_parametrised_by_nonsingular_coordinates_pt->
      boundary_parametrising_geom_object_pt(2);
     double zeta_start=disk_parametrised_by_nonsingular_coordinates_pt->
      zeta_boundary_start(2);
     double zeta_end=disk_parametrised_by_nonsingular_coordinates_pt->
      zeta_boundary_end(2);
     unsigned nsegment=half_nsegment;
     unsigned boundary_id=2;
     inner_curvilinear_boundary_pt[0]=new TriangleMeshCurviLine(
      inner_boundary_ellipse0_pt,zeta_start,zeta_end,nsegment,boundary_id);
     
     // Second bit
     GeomObject* inner_boundary_ellipse1_pt=
      disk_parametrised_by_nonsingular_coordinates_pt->
      boundary_parametrising_geom_object_pt(3);
     zeta_start=disk_parametrised_by_nonsingular_coordinates_pt->
      zeta_boundary_start(3);
     zeta_end=disk_parametrised_by_nonsingular_coordinates_pt->
      zeta_boundary_end(3);
     nsegment=half_nsegment;
     boundary_id=3;
     inner_curvilinear_boundary_pt[1]=new TriangleMeshCurviLine(
      inner_boundary_ellipse1_pt,zeta_start,zeta_end,nsegment,boundary_id);
 
     // Combine to curvilinear boundary and define the
     // inner boundary
     TriangleMeshClosedCurve* inner_curve_pt=
      new TriangleMeshClosedCurve(inner_curvilinear_boundary_pt);
 
 
     Vector<TriangleMeshClosedCurve *> inner_boundaries_pt(1);
     inner_boundaries_pt[0]=inner_curve_pt;
 
     // Specify the internal closed boundaries
     triangle_mesh_parameters.internal_closed_curve_pt() = inner_boundaries_pt;
 
 
 
     // Let's get a map of region coordinates (if any!)
     std::map<unsigned, Vector<double> > zeta_in_region_map=
      disk_parametrised_by_nonsingular_coordinates_pt->zeta_in_region();
     for (std::map<unsigned, Vector<double> >::iterator it=
           zeta_in_region_map.begin();it!=zeta_in_region_map.end();it++)
      {
       // Pass information about the defined regions
       unsigned r=(*it).first;       
       Vector<double> region_coords=(*it).second;
       triangle_mesh_parameters.add_region_coordinates(r,region_coords);
      }
 
    }
   
   // Specify the element area so that it matches the boundary discretisation
   double uniform_element_area=sqrt(3.0)/4.0*
    pow(2.0*MathematicalConstants::Pi/(2.0*double(half_nsegment)),2);
   triangle_mesh_parameters.element_area() = uniform_element_area;
 
   // Create the mesh 
   Tri_mesh_pt = new
    TriangleMesh<TPoissonElement<2,2> >(triangle_mesh_parameters);
 
   // Loop over all boundary elements and rotate their nodes so that
   // the first two nodes are on the boundary
   std::map<FiniteElement*,bool> is_on_boundary;
   for (unsigned b=0;b<2;b++)
    {
     unsigned nel=Tri_mesh_pt->nboundary_element(b);
     for (unsigned e=0;e<nel;e++)
      {
       FiniteElement* el_pt=Tri_mesh_pt->boundary_element_pt(b,e);
       unsigned count=0;
       for (unsigned j=0;j<3;j++)
        {
         Node* nod_pt=el_pt->node_pt(j);
         if (nod_pt->is_on_boundary()) count++;
        }
       if (count==2)
        {
         is_on_boundary[el_pt]=true;
         if ((el_pt->node_pt(0)->is_on_boundary())&&
             (el_pt->node_pt(1)->is_on_boundary()))
          {
           // fine
          }
         else
          {
           // Reorder nodes so that the first two nodes are on the boundary
           Node* nod0_pt=el_pt->node_pt(0);
           Node* nod1_pt=el_pt->node_pt(1);
           Node* nod2_pt=el_pt->node_pt(2);
           if (!el_pt->node_pt(0)->is_on_boundary())
            {
             el_pt->node_pt(0)=nod1_pt;
             el_pt->node_pt(1)=nod2_pt;
             el_pt->node_pt(2)=nod0_pt;
            }
           else if (!el_pt->node_pt(1)->is_on_boundary())
            {
             el_pt->node_pt(0)=nod2_pt;
             el_pt->node_pt(1)=nod0_pt;
             el_pt->node_pt(2)=nod1_pt;
            }
           else
            {
             std::ostringstream error_message;
             error_message << "This doesn't make sense!";
             throw OomphLibError(error_message.str(),
                                 OOMPH_CURRENT_FUNCTION,
                                 OOMPH_EXCEPTION_LOCATION);
            }
          }
        }
       else
        {
         std::ostringstream error_message;
         error_message <<
          "Boundary simplex element doesn't have two nodes on boundary!";
         throw OomphLibError(error_message.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
        }
      }
    }
 
 
   // Now loop over all nodes and turn them into vertices
   unsigned nnod=Tri_mesh_pt->nnode();
   Vertex_pt.resize(nnod);
   for (unsigned j=0;j<nnod;j++)
    {
     Node* nod_pt=Tri_mesh_pt->node_pt(j);
     Vector<double> sheet_point(3,0.0);
     Vector<double> zeta(nod_pt->position());
     Geom_object_with_boundaries_pt->position(zeta,sheet_point);
     Vertex_pt[j]=new TetMeshVertex(sheet_point);
     Vertex_pt[j]->set_zeta_in_geom_object(zeta);
     Equivalent_vertex_pt[nod_pt]=Vertex_pt[j];
    }
 
 
   // Quick count to see how many elements/facets are on the boundary
   Nelement_on_disk_boundary=0;
   unsigned nel=Tri_mesh_pt->nelement();
   for (unsigned e=0;e<nel;e++)
    {
     FiniteElement* el_pt=Tri_mesh_pt->finite_element_pt(e);
     if (is_on_boundary[el_pt])
      {
       Nelement_on_disk_boundary++;
      }
    }
 
   // Angles of right/left node on boundary (for snapping)
   Left_boundary_coordinate.resize(nel,0.0);
   Right_boundary_coordinate.resize(nel,0.0);
   Boundary_id.resize(nel,0.0);
 
   // Now loop over all elements
   Nfacet=nel;
   Facet_pt.resize(Nfacet);
   Facet_is_on_boundary.resize(Nfacet,false);
   unsigned n_vertex_on_facet=3;
   for (unsigned e=0;e<nel;e++)
    {
     FiniteElement* el_pt=Tri_mesh_pt->finite_element_pt(e);
     if (is_on_boundary[el_pt])
      {
       Node* left_node_pt=el_pt->node_pt(0);
       Node* right_node_pt=el_pt->node_pt(1);
       Vector<double> boundary_zeta(1);
       if (left_node_pt->is_on_boundary(0)&&
           right_node_pt->is_on_boundary(0))
        {
         unsigned b=0;
         left_node_pt->get_coordinates_on_boundary(b,boundary_zeta);
         Left_boundary_coordinate[e]=boundary_zeta[0];
         right_node_pt->get_coordinates_on_boundary(b,boundary_zeta);
         Right_boundary_coordinate[e]=boundary_zeta[0]; 
         Boundary_id[e]=b;
        }
       else if (left_node_pt->is_on_boundary(1)&&
                right_node_pt->is_on_boundary(1))
        {
         unsigned b=1;
         left_node_pt->get_coordinates_on_boundary(b,boundary_zeta);
         Left_boundary_coordinate[e]=boundary_zeta[0];
         right_node_pt->get_coordinates_on_boundary(b,boundary_zeta);
         Right_boundary_coordinate[e]=boundary_zeta[0];   
         Boundary_id[e]=b;
        }
       else
        {
         std::ostringstream error_message;
         error_message << "never get here!";
         throw OomphLibError(error_message.str(),
                             OOMPH_CURRENT_FUNCTION,
                             OOMPH_EXCEPTION_LOCATION);
        }
       Facet_is_on_boundary[e]=true;
      }
     Facet_pt[e]=new TetMeshFacet(n_vertex_on_facet);
     unsigned one_based_boundary_id=first_one_based_boundary_id_for_disk+e;
     Facet_pt[e]->set_one_based_boundary_id(one_based_boundary_id);
     last_one_based_boundary_id_for_disk=one_based_boundary_id;
     for (unsigned j=0;j<3;j++)
      {
       Facet_pt[e]->set_vertex_pt(j,Equivalent_vertex_pt[el_pt->node_pt(j)]);
      }
    }
  }
 
 
 ~DiskTetMeshFacetedSurface()
  {
   delete Tri_mesh_pt;
   Tri_mesh_pt=0;
  }
 
 void boundary_zeta01(const unsigned& facet_id,
                      const double& zeta_boundary,
                      Vector<double>& zeta)
 {
  if (Facet_is_on_boundary[facet_id])
   {
    double phi0=Left_boundary_coordinate[facet_id];
    double phi1=Right_boundary_coordinate[facet_id];
    double phi=phi0+(phi1-phi0)*zeta_boundary;
    
    unsigned b=Boundary_id[facet_id];
    double zeta_bound=phi;
    Geom_object_with_boundaries_pt->zeta_on_boundary(b,
                                                     zeta_bound,
                                                     zeta);
   }
  else
   {
    Vector<double> zeta0=Facet_pt[facet_id]->
     vertex_pt(0)->zeta_in_geom_object();
    Vector<double> zeta1=Facet_pt[facet_id]->
     vertex_pt(1)->zeta_in_geom_object();
    zeta[0]=zeta0[0]+(zeta1[0]-zeta0[0])*zeta_boundary;
    zeta[1]=zeta0[1]+(zeta1[1]-zeta0[1])*zeta_boundary;
   }
 }
 
 
protected:
 
 unsigned Nfacet;
 
 std::vector<bool> Facet_is_on_boundary;
 
 Vector<double> Left_boundary_coordinate;
 
 Vector<double> Right_boundary_coordinate;
 
 Vector<unsigned> Boundary_id;
 
 unsigned Nelement_on_disk_boundary;
 
 
 TriangleMesh<TPoissonElement<2,2> >* Tri_mesh_pt;
 
 
 std::map<Node*,TetMeshVertex*> Equivalent_vertex_pt;
 
};
 
 
 
 
 
 
 
 
//================================================================
//================================================================
class DiskWithTorusAroundEdgeTetMeshFacetedSurface : public virtual
DiskTetMeshFacetedSurface, public virtual TetMeshFacetedClosedSurface
{
 
 
public:
 
 DiskWithTorusAroundEdgeTetMeshFacetedSurface(
  DiskLikeGeomObjectWithBoundaries* 
  disk_parametrised_by_nonsingular_coordinates_pt,
  const unsigned& half_nsegment,
  const double& r_torus,
  const unsigned& nvertex_torus,
  const unsigned& first_one_based_boundary_id_for_disk,
  const unsigned& one_based_torus_region_id,
  unsigned& last_one_based_boundary_id_for_disk,
  unsigned& first_one_based_boundary_id_for_torus,
  unsigned& last_one_based_boundary_id_for_torus,
  Vector<unsigned>& one_based_boundary_id_for_disk_within_torus,
  Vector<unsigned>& one_based_boundary_id_for_disk_outside_torus) :
 DiskTetMeshFacetedSurface(
  disk_parametrised_by_nonsingular_coordinates_pt,
  half_nsegment,
  first_one_based_boundary_id_for_disk,
  last_one_based_boundary_id_for_disk)
  {
   // Wipe
   one_based_boundary_id_for_disk_within_torus.clear();
   one_based_boundary_id_for_disk_outside_torus.clear();
 
 
   // Is element in torus region or not?
   std::map<FiniteElement*,bool> is_in_torus_region;
   {
    unsigned r=1;
    /* ofstream some_file; */
    /* some_file.open("tri_mesh_in_region1.dat");  */
    unsigned nel=Tri_mesh_pt->nregion_element(r);
    for(unsigned e=0;e<nel;e++)
     {
      FiniteElement* el_pt=Tri_mesh_pt->region_element_pt(r,e);
      is_in_torus_region[el_pt]=true;
      //unsigned npts=3;
      //el_pt->output(some_file,npts);
     }
    //some_file.close();
   }
    
   
   // Now loop over all elements   
   unsigned nel=Tri_mesh_pt->nelement();
   for (unsigned e=0;e<nel;e++)  
    {
     FiniteElement* el_pt=Tri_mesh_pt->finite_element_pt(e); 
     unsigned one_based_boundary_id=Facet_pt[e]->one_based_boundary_id();
     if (is_in_torus_region[el_pt])
      {
       Facet_pt[e]->set_one_based_region_that_facet_is_embedded_in
        (one_based_torus_region_id);
       one_based_boundary_id_for_disk_within_torus.push_back
        (one_based_boundary_id);
      }
     else
      {
       one_based_boundary_id_for_disk_outside_torus.push_back
        (one_based_boundary_id);
      }
    }
 
 
   // Now add torus:
   //---------------
   unsigned one_based_boundary_id=last_one_based_boundary_id_for_disk+1;
   first_one_based_boundary_id_for_torus=one_based_boundary_id;
 
   // Find sorted sequence of vertices going around the
   // inner boundary (innermost part of torus)
   Vector<TetMeshVertex*> sorted_vertex_pt;
   Vector<std::pair<TetMeshVertex*,double> > vertex_pt_and_boundary_coord;
   std::map<Node*,bool> done;
   Vector<double> boundary_zeta(1);
   for (unsigned b=2;b<=3;b++)
    {
     unsigned nnod=Tri_mesh_pt->nboundary_node(b);
     for (unsigned j=0;j<nnod;j++)
      {
       Node* nod_pt=Tri_mesh_pt->boundary_node_pt(b,j);
       if (!done[nod_pt])
        {
         nod_pt->get_coordinates_on_boundary(b,boundary_zeta);
         vertex_pt_and_boundary_coord.push_back(
          std::make_pair(Equivalent_vertex_pt[nod_pt],
                         boundary_zeta[0]));
         done[nod_pt]=true;
        }
      }
    }
 
   
 
   oomph_info << "Number of vertices on inner ring: "
              <<  vertex_pt_and_boundary_coord.size() << " "
              << " sum: "
              << Tri_mesh_pt->nboundary_node(2)+Tri_mesh_pt->nboundary_node(3)
              << std::endl;
   
   // Identify point in torus (for tetgen)
   Vector<double> point_in_torus(3);
   unsigned b=0;
   unsigned j=0;
   Node* nod_pt=Tri_mesh_pt->boundary_node_pt(b,j);
   point_in_torus[0]=Equivalent_vertex_pt[nod_pt]->x(0);
   point_in_torus[1]=Equivalent_vertex_pt[nod_pt]->x(1);
   point_in_torus[2]=Equivalent_vertex_pt[nod_pt]->x(2)+0.5*r_torus;
 
   oomph_info << "point in torus: " 
              << point_in_torus[0] << " " 
              << point_in_torus[1] << " " 
              << point_in_torus[2] << " " 
              << std::endl;
 
  // Sort based on boundary coordinate
  std::sort (vertex_pt_and_boundary_coord.begin(),
             vertex_pt_and_boundary_coord.end(),
             HelperForSortingVectorOfPairsOfVertexPtAndDouble::
             less_than_based_on_double);
 
  // Loop over annular rings
  unsigned n_vertices_on_annular_ring=nvertex_torus;
  unsigned n_new_vertices_on_annular_ring=
   n_vertices_on_annular_ring-2;
 
 
  // Storage for vertices and facets
  unsigned nv=Vertex_pt.size();
  unsigned jv=nv;
  nv+=n_new_vertices_on_annular_ring*
   vertex_pt_and_boundary_coord.size();
  Vertex_pt.resize(nv);
 
  unsigned jf=Nfacet;
  Nfacet+=(n_new_vertices_on_annular_ring+1)*
   (vertex_pt_and_boundary_coord.size());
  Facet_pt.resize(Nfacet);
 
 
  unsigned new_vertex_count=0;
  unsigned new_facet_count=0;
  unsigned n=vertex_pt_and_boundary_coord.size();
  Vector<Vector<TetMeshVertex*> > annulus_vertex_pt(n);
  for (unsigned j=0;j<n;j++)
   {
 
    // Wrap around for simplicity
    annulus_vertex_pt[j].resize(n_vertices_on_annular_ring);
    annulus_vertex_pt[j][0]=vertex_pt_and_boundary_coord[j].first;
    annulus_vertex_pt[j][n_vertices_on_annular_ring-1]=
     vertex_pt_and_boundary_coord[j].first;
 
    // Position of vertex
    Vector<double> vertex_position(3);
    vertex_position[0]=vertex_pt_and_boundary_coord[j].first->x(0);
    vertex_position[1]=vertex_pt_and_boundary_coord[j].first->x(1);
    vertex_position[2]=vertex_pt_and_boundary_coord[j].first->x(2);
 
    // Polar angle
    double theta=vertex_pt_and_boundary_coord[j].second;
 
    // Which boundary is this point located on?
    unsigned b=0;
    if (theta>Geom_object_with_boundaries_pt->zeta_boundary_end(0))
     {
      b=1;
     }
    double zeta_bound=theta;
 
    // Get coordinate on edge and normal and tangent vectors
    Vector<double> x(3);
    Vector<double> r_edge(3);
    Vector<double> normal(3);
    Vector<double> tangent(3);
    Vector<double> normal_normal(3);
    Geom_object_with_boundaries_pt->boundary_triad(b,
                                                   zeta_bound,
                                                   r_edge,
                                                   tangent,
                                                   normal,
                                                   normal_normal);
    
    // What rho do we need to create a circle, centred at the edge
    // going through current vertex?
    Vector<double> distance_to_edge(3);
    distance_to_edge[0]=r_edge[0]-vertex_position[0];
    distance_to_edge[1]=r_edge[1]-vertex_position[1];
    distance_to_edge[2]=r_edge[2]-vertex_position[2];
    double rho=sqrt(distance_to_edge[0]*distance_to_edge[0]+
                    distance_to_edge[1]*distance_to_edge[1]+
                    distance_to_edge[2]*distance_to_edge[2]);
 
    // What is the starting angle?
    double cos_phi0=(distance_to_edge[0]*normal[0]+
                     distance_to_edge[1]*normal[1]+
                     distance_to_edge[2]*normal[2])/rho;
    if (cos_phi0>1.0) cos_phi0=1.0;
    if (cos_phi0<-1.0) cos_phi0=-1.0;
 
    double phi_0=acos(cos_phi0);
    
    for (unsigned i=1;i<n_vertices_on_annular_ring-1;i++)
     {
      double phi=phi_0+double(i)/double(n_vertices_on_annular_ring-1)*
       2.0*MathematicalConstants::Pi-MathematicalConstants::Pi;
      x[0]=r_edge[0]+rho*cos(phi)*normal[0]+rho*sin(phi)*normal_normal[0];
      x[1]=r_edge[1]+rho*cos(phi)*normal[1]+rho*sin(phi)*normal_normal[1];
      x[2]=r_edge[2]+rho*cos(phi)*normal[2]+rho*sin(phi)*normal_normal[2];
      
      TetMeshVertex* new_vertex_pt=new TetMeshVertex(x);
      annulus_vertex_pt[j][i]=new_vertex_pt;
      Vertex_pt[jv]=new_vertex_pt;
      jv++;
      new_vertex_count++;
     }
   }
 
  for (unsigned i=0;i<n-1;i++)
   {
    unsigned m=annulus_vertex_pt[i].size();
    for (unsigned j=0;j<m-1;j++)
     {
      unsigned n_vertex_on_facet=4;
      TetMeshFacet* new_facet_pt=new TetMeshFacet(n_vertex_on_facet);
      new_facet_count++;
      new_facet_pt->set_vertex_pt(0,annulus_vertex_pt[i][j]);
      new_facet_pt->set_vertex_pt(1,annulus_vertex_pt[i][j+1]);
      new_facet_pt->set_vertex_pt(2,annulus_vertex_pt[i+1][j+1]);
      new_facet_pt->set_vertex_pt(3,annulus_vertex_pt[i+1][j]);
      Facet_pt[jf]=new_facet_pt;
      unsigned region_id=one_based_torus_region_id;
      Facet_pt[jf]->set_one_based_adjacent_region_id(region_id);
      Facet_pt[jf]->set_one_based_boundary_id(one_based_boundary_id);
      one_based_boundary_id++;
      jf++;
     }
   }
  
  unsigned m=annulus_vertex_pt[n-1].size();
  for (unsigned j=0;j<m-1;j++)
   {
    unsigned n_vertex_on_facet=4;
    TetMeshFacet* new_facet_pt=new TetMeshFacet(n_vertex_on_facet);
    new_facet_count++;
    new_facet_pt->set_vertex_pt(0,annulus_vertex_pt[n-1][j]);
    new_facet_pt->set_vertex_pt(1,annulus_vertex_pt[n-1][j+1]);
    new_facet_pt->set_vertex_pt(2,annulus_vertex_pt[0][j+1]);
    new_facet_pt->set_vertex_pt(3,annulus_vertex_pt[0][j]);
    Facet_pt[jf]=new_facet_pt;
    unsigned region_id=one_based_torus_region_id;
    Facet_pt[jf]->set_one_based_adjacent_region_id(region_id);
    Facet_pt[jf]->set_one_based_boundary_id(one_based_boundary_id);
    last_one_based_boundary_id_for_torus=one_based_boundary_id;
    one_based_boundary_id++;
    jf++;
   }
 
  oomph_info << "First/last one-based disk boundary ID: "
             << first_one_based_boundary_id_for_disk << " "
             << last_one_based_boundary_id_for_disk << std::endl;
 
 
  oomph_info << "First/last one-based torus boundary ID: "
             << first_one_based_boundary_id_for_torus << " "
             << last_one_based_boundary_id_for_torus << std::endl;
 
 
   // Is it a genuine region?
   if (one_based_torus_region_id>0)
    {
     set_region_for_tetgen(one_based_torus_region_id-1,point_in_torus);
    }
   else
    {
     std::ostringstream error_message;
     error_message << "one_based_torus_region_id must be strictly positive";
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
 
 
 
  }
 
 
};
 
 
 
 
 
 
 
 
 
 
 
//================================================================
//================================================================
class DiskWithTwoLayersTetMeshFacetedSurface : public virtual
DiskTetMeshFacetedSurface, public virtual TetMeshFacetedClosedSurface
{
 
 
public:
 
 DiskWithTwoLayersTetMeshFacetedSurface(
  DiskLikeGeomObjectWithBoundaries* 
  disk_parametrised_by_nonsingular_coordinates_pt,
  const unsigned& half_nsegment,
  const unsigned& first_one_based_boundary_id_for_disk,
  const unsigned& one_based_region_id_above_disk,
  const unsigned& one_based_region_id_below_disk,
  unsigned& last_one_based_boundary_id_for_disk,
  unsigned& last_one_based_boundary_id) :
 DiskTetMeshFacetedSurface(disk_parametrised_by_nonsingular_coordinates_pt,
                           half_nsegment,
                           first_one_based_boundary_id_for_disk,
                           last_one_based_boundary_id_for_disk)
  {
 
   // Currently: Number of facets = number of elements on disk.
   // Set bounding region numbers for existing facets:
   unsigned nel=Nfacet;
   for (unsigned e=0;e<Nfacet;e++)
    {
     Facet_pt[e]->set_one_based_adjacent_region_id(
      one_based_region_id_above_disk);
     Facet_pt[e]->set_one_based_adjacent_region_id(
      one_based_region_id_below_disk);
    }
 
   // Additional facets: one square facet moving up/down from each
   // side and two lids.
   Nfacet+=2*Nelement_on_disk_boundary+2;
   Facet_pt.resize(Nfacet);
   Facet_is_on_boundary.resize(Nfacet,false);
 
   // Find sorted sequence of vertices going around the
   // perimeter of the disk
   Vector<TetMeshVertex*> sorted_perimeter_vertex_pt;
   std::set<TetMeshFacet*> boundary_facet_pt;
   bool first=true;
   for (unsigned e=0;e<nel;e++)
    {
     if (Facet_is_on_boundary[e])
      {
       boundary_facet_pt.insert(Facet_pt[e]);
       if (first)
        {
         TetMeshVertex* vertex0_pt=Facet_pt[e]->vertex_pt(0);
         sorted_perimeter_vertex_pt.push_back(vertex0_pt);
         first=false;
        }
      }
    }
   while (boundary_facet_pt.size()!=1)
    {
     for (std::set<TetMeshFacet*>::iterator it=boundary_facet_pt.begin();
          it!=boundary_facet_pt.end();it++)
      {
       unsigned n_last=sorted_perimeter_vertex_pt.size();
       if (sorted_perimeter_vertex_pt[n_last-1]==
           (*it)->vertex_pt(0))
        {
         sorted_perimeter_vertex_pt.push_back((*it)->vertex_pt(1));
         boundary_facet_pt.erase(*it);
         break;
        }
      }
    }
 
   //Find mean position
   double z_mean=0.0;
   double z_min= DBL_MAX;
   double z_max=-DBL_MAX;
   unsigned n=sorted_perimeter_vertex_pt.size();
   for (unsigned j=0;j<n;j++)
    {
     double z=sorted_perimeter_vertex_pt[j]->x(2);
     z_mean+=z;
     if (z>z_max) z_max=z;
     if (z<z_min) z_min=z;
    }
   z_mean/=double(n);
   double layer_thickness=0.6*(z_max-z_min);
   
   // Loop over upper/lower layer
   unsigned facet_count=nel;
   int offset_sign=-1;
   unsigned one_based_boundary_id=last_one_based_boundary_id_for_disk+1;
   unsigned n_vertex_on_disk=Vertex_pt.size();
   unsigned n_vertex_on_disk_perimeter=sorted_perimeter_vertex_pt.size();
   for (unsigned i_offset=0;i_offset<2;i_offset++)
    {
     unsigned one_based_region_id=one_based_region_id_below_disk;
     if (i_offset==1) one_based_region_id=one_based_region_id_above_disk;
 
     for (unsigned j=0;j<n_vertex_on_disk_perimeter;j++)
      {
       Vector<double> posn(3);
       posn[0]=sorted_perimeter_vertex_pt[j]->x(0);
       posn[1]=sorted_perimeter_vertex_pt[j]->x(1);
       posn[2]=z_mean+double(offset_sign)*layer_thickness;
       Vertex_pt.push_back(new TetMeshVertex(posn));
      }
     unsigned count=0;
     for (unsigned j=0;j<n_vertex_on_disk_perimeter-1;j++)
      {
       Facet_pt[facet_count]=new TetMeshFacet(4);
       Facet_pt[facet_count]->set_one_based_boundary_id(one_based_boundary_id);
       Facet_pt[facet_count]->set_vertex_pt(0,sorted_perimeter_vertex_pt[j]);
       Facet_pt[facet_count]->set_vertex_pt(1,sorted_perimeter_vertex_pt[j+1]);
       Facet_pt[facet_count]->
        set_vertex_pt(2,Vertex_pt[n_vertex_on_disk+
                                  i_offset*n_vertex_on_disk_perimeter+j+1]);
       Facet_pt[facet_count]->
        set_vertex_pt(3,Vertex_pt[n_vertex_on_disk+
                                  i_offset*n_vertex_on_disk_perimeter+j  ]);
       Facet_pt[facet_count]->
        set_one_based_adjacent_region_id(one_based_region_id);
       count++;
       facet_count++;
       one_based_boundary_id++;
      }
     Facet_pt[facet_count]=new TetMeshFacet(4);
     Facet_pt[facet_count]->set_one_based_boundary_id(one_based_boundary_id);
     Facet_pt[facet_count]->set_vertex_pt(0,sorted_perimeter_vertex_pt[count]);
     Facet_pt[facet_count]->set_vertex_pt(1,sorted_perimeter_vertex_pt[0]);
     Facet_pt[facet_count]->
      set_vertex_pt(2,Vertex_pt[n_vertex_on_disk+
                                i_offset*n_vertex_on_disk_perimeter+0]);
     Facet_pt[facet_count]->
      set_vertex_pt(3,Vertex_pt[n_vertex_on_disk+
                                i_offset*n_vertex_on_disk_perimeter+count]);
     Facet_pt[facet_count]->
      set_one_based_adjacent_region_id(one_based_region_id);
 
 
     // Closing lid
     facet_count++;
     one_based_boundary_id++;
     Facet_pt[facet_count]=new TetMeshFacet(n_vertex_on_disk_perimeter);
     Facet_pt[facet_count]->set_one_based_boundary_id(one_based_boundary_id);
     for (unsigned j=0;j<n_vertex_on_disk_perimeter;j++)
      {
       Facet_pt[facet_count]->
        set_vertex_pt(j,Vertex_pt[n_vertex_on_disk+
                                  i_offset*n_vertex_on_disk_perimeter+j]);
      }
     Facet_pt[facet_count]->
      set_one_based_adjacent_region_id(one_based_region_id);
 
     // Last one used (gets assigned twice but last one sticks!)
     last_one_based_boundary_id=one_based_boundary_id;
 
     // Bump
     facet_count++;
     one_based_boundary_id++;
 
     Vector<double> region_point(3,0.0);
     region_point[2]=z_mean+0.5*double(offset_sign)*layer_thickness;
     set_region_for_tetgen(one_based_region_id-1,region_point);
 
     offset_sign*=-1;
    }
  }
 
 
 
};
 
}
 
 
#endif