TetMeshFacetedSurface that defines disk. More...

#include <tetmesh_faceted_surfaces.h>

Inheritance diagram for oomph::DiskTetMeshFacetedSurface:

Public Member Functions
	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)

	~DiskTetMeshFacetedSurface ()
	Destructor. More...

void	boundary_zeta01 (const unsigned &facet_id, const double &zeta_boundary, Vector< double > &zeta)

Public Member Functions inherited from oomph::TetMeshFacetedSurface
	TetMeshFacetedSurface ()
	Constructor: More...

virtual	~TetMeshFacetedSurface ()
	Empty destructor. More...

unsigned	nvertex () const
	Number of vertices. More...

unsigned	nfacet () const
	Number of facets. More...

unsigned	one_based_facet_boundary_id (const unsigned &j) const
	One-based boundary id of j-th facet. More...

unsigned	one_based_vertex_boundary_id (const unsigned &j) const
	First (of possibly multiple) one-based boundary id of j-th vertex. More...

double	vertex_coordinate (const unsigned &j, const unsigned &i) const
	i-th coordinate of j-th vertex More...

unsigned	nvertex_on_facet (const unsigned &j) const
	Number of vertices defining the j-th facet. More...

bool	boundaries_can_be_split_in_tetgen ()
	Test whether boundary can be split in tetgen. More...

void	enable_boundaries_can_be_split_in_tetgen ()
	Test whether boundaries can be split in tetgen. More...

void	disable_boundaries_can_be_split_in_tetgen ()
	Test whether boundaries can be split in tetgen. More...

TetMeshFacet *	facet_pt (const unsigned &j) const
	Pointer to j-th facet. More...

TetMeshVertex *	vertex_pt (const unsigned &j) const
	Pointer to j-th vertex. More...

DiskLikeGeomObjectWithBoundaries *	geom_object_with_boundaries_pt ()

void	output (std::ostream &outfile) const
	Output. More...

void	output (const std::string &filename) const
	Output. More...

virtual void	boundary_zeta12 (const unsigned &facet_id, const double &zeta_boundary, Vector< double > &zeta)

virtual void	boundary_zeta20 (const unsigned &facet_id, const double &zeta_boundary, Vector< double > &zeta)

Vector< unsigned >	vertex_index_in_tetgen (const unsigned &f)

Protected Attributes
unsigned	Nfacet
	Number of facets. More...

std::vector< bool >	Facet_is_on_boundary
	Is facet on boundary? More...

Vector< double >	Left_boundary_coordinate
	Left boundary coordinate of i-th facet. More...

Vector< double >	Right_boundary_coordinate
	Right boundary coordinate of i-th facet. More...

Vector< unsigned >	Boundary_id
	ID of boundary the i-th facet is located on. More...

unsigned	Nelement_on_disk_boundary
	Number of elements on disk boundary. More...

TriangleMesh< TPoissonElement< 2, 2 > > *	Tri_mesh_pt
	Mesh used to facet-ise (discretise) disk. More...

std::map< Node , TetMeshVertex >	Equivalent_vertex_pt
	Mapping between nodes and vertices. More...

Protected Attributes inherited from oomph::TetMeshFacetedSurface
Vector< TetMeshVertex * >	Vertex_pt
	Vector pointers to vertices. More...

Vector< TetMeshFacet * >	Facet_pt
	Vector of pointers to facets. More...

bool	Boundaries_can_be_split_in_tetgen

Vector< Vector< unsigned > >	Facet_vertex_index_in_tetgen

DiskLikeGeomObjectWithBoundaries *	Geom_object_with_boundaries_pt
	GeomObject with boundaries associated with this surface. More...

Detailed Description

TetMeshFacetedSurface that defines disk.

////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////

Constructor & Destructor Documentation

◆ DiskTetMeshFacetedSurface()

oomph::DiskTetMeshFacetedSurface::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
	)

inline

Constructor: Pass pointer to GeomObject (with boundaries, and parametrised by coordinates without coordinate singularities, i.e. not polars, say) that defines the shape of the disk. Other args specify half the number of segments on perimeter of disk and first one-based boundary ID to be used to enumerate the boundaries on the disk. Returns last one-based boundary id used to enumerate the 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)]);
       }
     }
   }

References oomph::TriangleMeshParameters::add_region_coordinates(), b, Boundary_id, e(), oomph::TriangleMeshParameters::element_area(), Equivalent_vertex_pt, Facet_is_on_boundary, oomph::TetMeshFacetedSurface::Facet_pt, oomph::TetMeshFacetedSurface::Geom_object_with_boundaries_pt, oomph::Node::get_coordinates_on_boundary(), oomph::TriangleMeshParameters::internal_closed_curve_pt(), oomph::Node::is_on_boundary(), j, Left_boundary_coordinate, oomph::DiskLikeGeomObjectWithBoundaries::nboundary(), Nelement_on_disk_boundary, Nfacet, oomph::FiniteElement::node_pt(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION, oomph::MathematicalConstants::Pi, oomph::GeomObject::position(), oomph::Node::position(), Eigen::bfloat16_impl::pow(), UniformPSDSelfTest::r, Right_boundary_coordinate, sqrt(), Tri_mesh_pt, oomph::TetMeshFacetedSurface::Vertex_pt, Eigen::zeta(), and oomph::DiskLikeGeomObjectWithBoundaries::zeta_in_region().

◆ ~DiskTetMeshFacetedSurface()

oomph::DiskTetMeshFacetedSurface::~DiskTetMeshFacetedSurface ( )

inline

Destructor.

   {
    delete Tri_mesh_pt;
    Tri_mesh_pt=0;
   }

References Tri_mesh_pt.

Member Function Documentation

◆ boundary_zeta01()

void oomph::DiskTetMeshFacetedSurface::boundary_zeta01	(	const unsigned &	facet_id,
		const double &	zeta_boundary,
		Vector< double > &	zeta
	)

inlinevirtual

Function that specifies the variation of the zeta coordinates in the GeomObject along the boundary connecting vertices 0 and 1 in the facet. NOTE: input zeta_boundary ranges between 0 and 1; gets mapped to actual boundary coordinate inside!

Reimplemented from oomph::TetMeshFacetedSurface.

  {
   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;
    }
  }