tetgen_mesh.template.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//====================================================================
 
#ifndef OOMPH_TETGEN_MESH_HEADER
#define OOMPH_TETGEN_MESH_HEADER
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
 
#ifdef OOMPH_HAS_MPI
// mpi headers
#include "mpi.h"
#endif
 
#include "../generic/tetgen_scaffold_mesh.h"
#include "../generic/tet_mesh.h"
 
namespace oomph
{
  //=========start of TetgenMesh class======================================
  //========================================================================
  template<class ELEMENT>
  class TetgenMesh : public virtual TetMeshBase
  {
  public:
    TetgenMesh()
    {
      // Mesh can only be built with 3D Telements.
      MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
    }
 
    TetgenMesh(const std::string& node_file_name,
               const std::string& element_file_name,
               const std::string& face_file_name,
               TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
               const bool& use_attributes = false)
      : Tetgenio_exists(false), Tetgenio_pt(nullptr)
    {
      // Mesh can only be built with 3D Telements.
      MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
 
      // Store the attributes
      Use_attributes = use_attributes;
 
      // Store timestepper used to build elements
      Time_stepper_pt = time_stepper_pt;
 
      // Build scaffold
      Tmp_mesh_pt = new TetgenScaffoldMesh(
        node_file_name, element_file_name, face_file_name);
 
      // Convert mesh from scaffold to actual mesh
      build_from_scaffold(time_stepper_pt, use_attributes);
 
      // Kill the scaffold
      delete Tmp_mesh_pt;
      Tmp_mesh_pt = 0;
 
      // Setup boundary coordinates
      unsigned nb = nboundary();
      for (unsigned b = 0; b < nb; b++)
      {
        bool switch_normal = false;
        setup_boundary_coordinates<ELEMENT>(b, switch_normal);
      }
    }
 
 
    TetgenMesh(tetgenio& tetgen_data,
               TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
               const bool& use_attributes = false)
 
    {
      // Mesh can only be built with 3D Telements.
      MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
 
      // Store the attributes
      Use_attributes = use_attributes;
 
      // Store timestepper used to build elements
      Time_stepper_pt = time_stepper_pt;
 
      // We do not have a tetgenio representation
      Tetgenio_exists = false;
      Tetgenio_pt = 0;
 
      // Build scaffold
      Tmp_mesh_pt = new TetgenScaffoldMesh(tetgen_data);
 
      // Convert mesh from scaffold to actual mesh
      build_from_scaffold(time_stepper_pt, use_attributes);
 
      // Kill the scaffold
      delete Tmp_mesh_pt;
      Tmp_mesh_pt = 0;
 
      // Setup boundary coordinates
      unsigned nb = nboundary();
      for (unsigned b = 0; b < nb; b++)
      {
        bool switch_normal = false;
        setup_boundary_coordinates<ELEMENT>(b, switch_normal);
      }
    }
 
 
    TetgenMesh(const std::string& node_file_name,
               const std::string& element_file_name,
               const std::string& face_file_name,
               const bool& split_corner_elements,
               TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
               const bool& use_attributes = false)
 
    {
      // Mesh can only be built with 3D Telements.
      MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
 
      // Store the attributes
      Use_attributes = use_attributes;
 
      // Store timestepper used to build elements
      Time_stepper_pt = time_stepper_pt;
 
      // We do not have a tetgenio representation
      this->Tetgenio_exists = false;
      this->Tetgenio_pt = 0;
 
      // Build scaffold
      Tmp_mesh_pt = new TetgenScaffoldMesh(
        node_file_name, element_file_name, face_file_name);
 
      // Convert mesh from scaffold to actual mesh
      build_from_scaffold(time_stepper_pt, use_attributes);
 
      // Kill the scaffold
      delete Tmp_mesh_pt;
      Tmp_mesh_pt = 0;
 
      // Split corner elements
      if (split_corner_elements)
      {
        split_elements_in_corners<ELEMENT>();
      }
 
      // Setup boundary coordinates
      unsigned nb = nboundary();
      for (unsigned b = 0; b < nb; b++)
      {
        bool switch_normal = false;
        setup_boundary_coordinates<ELEMENT>(b, switch_normal);
      }
    }
 
    TetgenMesh(tetgenio& tetgen_data,
               const bool& split_corner_elements,
               TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
               const bool& use_attributes = false)
 
    {
      // Mesh can only be built with 3D Telements.
      MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
 
      // Store the attributes
      Use_attributes = use_attributes;
 
      // Store timestepper used to build elements
      Time_stepper_pt = time_stepper_pt;
 
      // We do not have a tetgenio representation
      this->Tetgenio_exists = false;
      this->Tetgenio_pt = nullptr;
 
      // Build scaffold
      Tmp_mesh_pt = new TetgenScaffoldMesh(tetgen_data);
 
      // Convert mesh from scaffold to actual mesh
      build_from_scaffold(time_stepper_pt, use_attributes);
 
      // Kill the scaffold
      delete Tmp_mesh_pt;
      Tmp_mesh_pt = nullptr;
 
      // Split corner elements
      if (split_corner_elements)
      {
        split_elements_in_corners<ELEMENT>();
      }
 
      // Setup boundary coordinates
      unsigned nb = nboundary();
      for (unsigned b = 0; b < nb; b++)
      {
        bool switch_normal = false;
        setup_boundary_coordinates<ELEMENT>(b, switch_normal);
      }
    }
 
 
    TetgenMesh(
      TetMeshFacetedClosedSurface* const& outer_boundary_pt,
      Vector<TetMeshFacetedSurface*>& internal_surface_pt,
      const double& element_volume,
      TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
      const bool& use_attributes = false,
      const bool& split_corner_elements = false,
      Vector<double>* const& target_element_volume_in_region_pt = nullptr)
    {
      // Mesh can only be built with 3D Telements.
      MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
 
      // Store the attributes
      Use_attributes = use_attributes;
 
      // Store timestepper used to build elements
      Time_stepper_pt = time_stepper_pt;
 
      // Copy across
      Outer_boundary_pt = outer_boundary_pt;
      // Setup the reverse lookup scheme
      this->setup_reverse_lookup_schemes_for_faceted_surface(Outer_boundary_pt);
      // Store the internal boundary
      Internal_surface_pt = internal_surface_pt;
      // Setup the reverse lookup schemes
      {
        unsigned n = this->Internal_surface_pt.size();
        for (unsigned i = 0; i < n; i++)
        {
          this->setup_reverse_lookup_schemes_for_faceted_surface(
            Internal_surface_pt[i]);
        }
      }
 
      // Tetgen data structure for the input and output
      tetgenio in;
      this->build_tetgenio(outer_boundary_pt,
                           internal_surface_pt,
                           target_element_volume_in_region_pt,
                           in);
 
 
      // Now tetrahedralise
 
      // The 'p' switch reads a Piecewise Linear Complex, which generates a
      // Delaunay tetrahedralization of the input. The 'q' switch prevents
      // generation of high-aspect ratio tets (slivers), with the trailing float
      // indicating the maximum allowed aspect ratio (default is 2.0). The 'a'
      // switch without  subsequent floating-point number switches on the
      // specific element volume constraints for particular regions (the 5th
      // index in the tetgenio.regionlist array). The 'A' enables region
      // attributes, and the second 'a' switch with the subsequent float sets
      // the global (non-region-specific) element volume constraint.
      std::stringstream input_string;
      input_string << "pq2.0aAa" << element_volume;
 
      // input_string << "Vpq1.414Aa" << element_volume;
      // << "Q"; // Q for quiet!
      // << "V"; // V for verbose incl. quality output!
 
      // If any of the boundaries should not be split add the "Y" flag
      bool can_boundaries_be_split =
        outer_boundary_pt->boundaries_can_be_split_in_tetgen();
      {
        unsigned n_internal = internal_surface_pt.size();
        for (unsigned i = 0; i < n_internal; i++)
        {
          can_boundaries_be_split &=
            internal_surface_pt[i]->boundaries_can_be_split_in_tetgen();
        }
      }
 
      // If we can't split the boundaries add the flag
      if (can_boundaries_be_split == false)
      {
        input_string << "Y";
      }
 
      // Now convert to a C-style string
      char tetswitches[100];
      sprintf(tetswitches, "%s", input_string.str().c_str());
 
      // Make a new tetgen representation
      this->Tetgenio_exists = true;
      Tetgenio_pt = new tetgenio;
      tetrahedralize(tetswitches, &in, this->Tetgenio_pt);
 
      // Build scaffold
      Tmp_mesh_pt = new TetgenScaffoldMesh(*this->Tetgenio_pt);
 
      // If any of the objects are different regions then we need to use
      // the atributes
      bool regions_exist = false;
      {
        unsigned n_internal = internal_surface_pt.size();
        for (unsigned i = 0; i < n_internal; i++)
        {
          TetMeshFacetedClosedSurface* srf_pt =
            dynamic_cast<TetMeshFacetedClosedSurface*>(internal_surface_pt[i]);
          if (srf_pt != 0)
          {
            unsigned n_int_pts = srf_pt->ninternal_point_for_tetgen();
            for (unsigned j = 0; j < n_int_pts; j++)
            {
              regions_exist |=
                srf_pt->internal_point_identifies_region_for_tetgen(j);
            }
          }
        }
      }
 
      // If there are regions, use the attributes
      if (regions_exist)
      {
        Use_attributes = true;
      }
 
      // Convert mesh from scaffold to actual mesh
      build_from_scaffold(time_stepper_pt, Use_attributes);
 
      // Kill the scaffold
      delete Tmp_mesh_pt;
      Tmp_mesh_pt = 0;
 
      // Split corner elements
      if (split_corner_elements)
      {
        split_elements_in_corners<ELEMENT>();
      }
 
      // Setup boundary coordinates
      unsigned nb = nboundary();
      for (unsigned b = 0; b < nb; b++)
      {
        bool switch_normal = false;
        setup_boundary_coordinates<ELEMENT>(b, switch_normal);
      }
 
      // Now snap onto geometric objects associated with triangular facets
      // (if any!)
      snap_nodes_onto_geometric_objects();
    }
 
    void build_tetgenio(
      TetMeshFacetedSurface* const& outer_boundary_pt,
      Vector<TetMeshFacetedSurface*>& internal_surface_pt,
      Vector<double>* const& target_element_volume_in_region_pt,
      tetgenio& tetgen_io)
    {
      // Pointer to Tetgen facet
      tetgenio::facet* f;
      // Pointer to Tetgen polygon
      tetgenio::polygon* p;
 
      // Start all indices from 1 (it's a choice and we've made it
      tetgen_io.firstnumber = 1;
      tetgen_io.useindex = true;
 
      // Find the number of internal surfaces
      const unsigned n_internal = internal_surface_pt.size();
 
      // Find the number of points on the outer boundary
      const unsigned n_outer_vertex = outer_boundary_pt->nvertex();
      tetgen_io.numberofpoints = n_outer_vertex;
 
      // Find the number of points on the inner boundaries and add to the totals
      Vector<unsigned> n_internal_vertex(n_internal);
      Vector<unsigned> internal_vertex_offset(n_internal);
      for (unsigned h = 0; h < n_internal; ++h)
      {
        n_internal_vertex[h] = internal_surface_pt[h]->nvertex();
        internal_vertex_offset[h] = tetgen_io.numberofpoints;
        tetgen_io.numberofpoints += n_internal_vertex[h];
      }
 
      // Read the data into the point list
      tetgen_io.pointlist = new double[tetgen_io.numberofpoints * 3];
      tetgen_io.pointmarkerlist = new int[tetgen_io.numberofpoints];
      unsigned counter = 0;
      for (unsigned n = 0; n < n_outer_vertex; ++n)
      {
        for (unsigned i = 0; i < 3; ++i)
        {
          tetgen_io.pointlist[counter] =
            outer_boundary_pt->vertex_coordinate(n, i);
          ++counter;
        }
      }
      for (unsigned h = 0; h < n_internal; ++h)
      {
        const unsigned n_inner = n_internal_vertex[h];
        for (unsigned n = 0; n < n_inner; ++n)
        {
          for (unsigned i = 0; i < 3; ++i)
          {
            tetgen_io.pointlist[counter] =
              internal_surface_pt[h]->vertex_coordinate(n, i);
            ++counter;
          }
        }
      }
 
 
      // Set up the pointmarkers
      counter = 0;
      for (unsigned n = 0; n < n_outer_vertex; ++n)
      {
        tetgen_io.pointmarkerlist[counter] =
          outer_boundary_pt->one_based_vertex_boundary_id(n);
        ++counter;
      }
      for (unsigned h = 0; h < n_internal; ++h)
      {
        const unsigned n_inner = n_internal_vertex[h];
        for (unsigned n = 0; n < n_inner; ++n)
        {
          tetgen_io.pointmarkerlist[counter] =
            internal_surface_pt[h]->one_based_vertex_boundary_id(n);
          ++counter;
        }
      }
 
 
      // Now the facets
      unsigned n_outer_facet = outer_boundary_pt->nfacet();
      tetgen_io.numberoffacets = n_outer_facet;
      Vector<unsigned> n_inner_facet(n_internal);
      for (unsigned h = 0; h < n_internal; ++h)
      {
        n_inner_facet[h] = internal_surface_pt[h]->nfacet();
        tetgen_io.numberoffacets += n_inner_facet[h];
      }
 
      tetgen_io.facetlist = new tetgenio::facet[tetgen_io.numberoffacets];
      tetgen_io.facetmarkerlist = new int[tetgen_io.numberoffacets];
 
 
      counter = 0;
      for (unsigned n = 0; n < n_outer_facet; ++n)
      {
        // Set pointer to facet
        f = &tetgen_io.facetlist[counter];
        f->numberofpolygons = 1;
        f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
        f->numberofholes = 0;
        f->holelist = NULL;
        p = &f->polygonlist[0];
 
        Vector<unsigned> facet = outer_boundary_pt->vertex_index_in_tetgen(n);
 
        p->numberofvertices = facet.size();
        p->vertexlist = new int[p->numberofvertices];
        for (int i = 0; i < p->numberofvertices; ++i)
        {
          // The offset here is because we have insisted on 1-based indexing
          p->vertexlist[i] = facet[i] + 1;
        }
 
        // Set up the boundary markers
        tetgen_io.facetmarkerlist[counter] =
          outer_boundary_pt->one_based_facet_boundary_id(n);
        // Increase the counter
        ++counter;
      }
 
      // Initialise the number of holes
      tetgen_io.numberofholes = 0;
      // and the number of regions
      tetgen_io.numberofregions = 0;
 
      // Loop over the internal stuff
      for (unsigned h = 0; h < n_internal; ++h)
      {
        for (unsigned n = 0; n < n_inner_facet[h]; ++n)
        {
          // Set pointer to facet
          f = &tetgen_io.facetlist[counter];
          f->numberofpolygons = 1;
          f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
          f->numberofholes = 0;
          f->holelist = NULL;
          p = &f->polygonlist[0];
 
          Vector<unsigned> facet =
            internal_surface_pt[h]->vertex_index_in_tetgen(n);
 
          p->numberofvertices = facet.size();
          p->vertexlist = new int[p->numberofvertices];
          for (int i = 0; i < p->numberofvertices; ++i)
          {
            // Add the 1-based and vertex offsets to get these number correct
            p->vertexlist[i] = facet[i] + internal_vertex_offset[h] + 1;
          }
          // Set up the boundary markers
          tetgen_io.facetmarkerlist[counter] =
            internal_surface_pt[h]->one_based_facet_boundary_id(n);
          ++counter;
        }
 
        // If it's a hole/region add it
        TetMeshFacetedClosedSurface* srf_pt =
          dynamic_cast<TetMeshFacetedClosedSurface*>(internal_surface_pt[h]);
        if (srf_pt != 0)
        {
          unsigned n_int_pts = srf_pt->ninternal_point_for_tetgen();
          for (unsigned j = 0; j < n_int_pts; j++)
          {
            if (srf_pt->internal_point_identifies_hole_for_tetgen(j))
            {
              ++tetgen_io.numberofholes;
            }
            // Otherwise it may be region
            else
            {
              if (srf_pt->internal_point_identifies_region_for_tetgen(j))
              {
                ++tetgen_io.numberofregions;
              }
            }
          }
        }
      }
 
      // Set storage for the holes
      tetgen_io.holelist = new double[3 * tetgen_io.numberofholes];
 
      // Loop over all the internal boundaries again
      counter = 0;
      for (unsigned h = 0; h < n_internal; ++h)
      {
        TetMeshFacetedClosedSurface* srf_pt =
          dynamic_cast<TetMeshFacetedClosedSurface*>(internal_surface_pt[h]);
        if (srf_pt != 0)
        {
          unsigned n_int_pts = srf_pt->ninternal_point_for_tetgen();
          for (unsigned j = 0; j < n_int_pts; j++)
          {
            if (srf_pt->internal_point_identifies_hole_for_tetgen(j))
            {
              for (unsigned i = 0; i < 3; ++i)
              {
                tetgen_io.holelist[counter] =
                  srf_pt->internal_point_for_tetgen(j, i);
                ++counter;
              }
            }
          }
        }
      }
 
      // Set storage for the regions
      tetgen_io.regionlist = new double[5 * tetgen_io.numberofregions];
 
      // Loop over all the internal boundaries again
      counter = 0;
      for (unsigned h = 0; h < n_internal; ++h)
      {
        TetMeshFacetedClosedSurface* srf_pt =
          dynamic_cast<TetMeshFacetedClosedSurface*>(internal_surface_pt[h]);
        if (srf_pt != 0)
        {
          unsigned n_int_pts = srf_pt->ninternal_point_for_tetgen();
          for (unsigned j = 0; j < n_int_pts; j++)
          {
            if (srf_pt->internal_point_identifies_region_for_tetgen(j))
            {
              for (unsigned i = 0; i < 3; ++i)
              {
                tetgen_io.regionlist[counter] =
                  srf_pt->internal_point_for_tetgen(j, i);
                ++counter;
              }
              tetgen_io.regionlist[counter] =
                static_cast<double>(srf_pt->region_id_for_tetgen(j));
              ++counter;
 
              // if there's no target volumes specified, default to zero
              if (target_element_volume_in_region_pt == nullptr)
              {
                tetgen_io.regionlist[counter] = 0;
              }
              else
              {
                // deliberate integer division here to round down to the region
                // number (five doubles per region)
                tetgen_io.regionlist[counter] =
                  (*target_element_volume_in_region_pt)[unsigned(counter / 5)];
              }
 
              ++counter;
            }
          }
        }
      }
    }
 
 
    ~TetgenMesh()
    {
      if (Tetgenio_exists)
      {
        delete Tetgenio_pt;
      }
    }
 
 
    void set_mesh_level_time_stepper(TimeStepper* const& time_stepper_pt,
                                     const bool& preserve_existing_data)
    {
      this->Time_stepper_pt = time_stepper_pt;
    }
 
 
    bool tetgenio_exists() const
    {
      return Tetgenio_exists;
    }
 
 
    tetgenio*& tetgenio_pt()
    {
      return Tetgenio_pt;
    }
 
    void set_deep_copy_tetgenio_pt(tetgenio* const& tetgenio_pt)
    {
      // Delete the existing data
      if (Tetgenio_exists)
      {
        delete Tetgenio_pt;
      }
      this->Tetgenio_pt = new tetgenio;
      // Create a deep copy of tetgenio_pt and store the result in
      // Tetgenio_pt
      this->deep_copy_of_tetgenio(tetgenio_pt, this->Tetgenio_pt);
    }
 
    void deep_copy_of_tetgenio(tetgenio* const& input_pt, tetgenio*& output_pt);
 
 
  protected:
    void build_from_scaffold(TimeStepper* time_stepper_pt,
                             const bool& use_attributes);
 
    void setup_reverse_lookup_schemes_for_faceted_surface(
      TetMeshFacetedSurface* const& faceted_surface_pt);
 
    TetgenScaffoldMesh* Tmp_mesh_pt;
 
    bool Tetgenio_exists;
 
    tetgenio* Tetgenio_pt;
 
    bool Use_attributes;
 
  }; // end class
 
 
 
 
  //==============start_mesh=================================================
  //=========================================================================
  template<class ELEMENT>
  class SolidTetgenMesh : public virtual TetgenMesh<ELEMENT>,
                          public virtual SolidMesh
  {
  public:
    SolidTetgenMesh(const std::string& node_file_name,
                    const std::string& element_file_name,
                    const std::string& face_file_name,
                    const bool& split_corner_elements,
                    TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
                    const bool& use_attributes = false)
      : TetgenMesh<ELEMENT>(node_file_name,
                            element_file_name,
                            face_file_name,
                            split_corner_elements,
                            time_stepper_pt,
                            use_attributes)
    {
      // Assign the Lagrangian coordinates
      set_lagrangian_nodal_coordinates();
    }
 
    SolidTetgenMesh(const std::string& node_file_name,
                    const std::string& element_file_name,
                    const std::string& face_file_name,
                    const bool& split_corner_elements,
                    const bool& switch_normal,
                    TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
                    const bool& use_attributes = false)
      : TetgenMesh<ELEMENT>(node_file_name,
                            element_file_name,
                            face_file_name,
                            split_corner_elements,
                            time_stepper_pt,
                            use_attributes)
    {
      // Assign the Lagrangian coordinates
      set_lagrangian_nodal_coordinates();
 
      // Re-setup boundary coordinates for all boundaries with specified
      // orientation of nnormal
      unsigned nb = this->nboundary();
      for (unsigned b = 0; b < nb; b++)
      {
        this->template setup_boundary_coordinates<ELEMENT>(b, switch_normal);
      }
    }
 
    virtual ~SolidTetgenMesh() {}
  };
 
 
} // namespace oomph
 
#endif