Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
       3             : 
       4             : // This library is free software; you can redistribute it and/or
       5             : // modify it under the terms of the GNU Lesser General Public
       6             : // License as published by the Free Software Foundation; either
       7             : // version 2.1 of the License, or (at your option) any later version.
       8             : 
       9             : // This library is distributed in the hope that it will be useful,
      10             : // but WITHOUT ANY WARRANTY; without even the implied warranty of
      11             : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      12             : // Lesser General Public License for more details.
      13             : 
      14             : // You should have received a copy of the GNU Lesser General Public
      15             : // License along with this library; if not, write to the Free Software
      16             : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
      17             : 
      18             : 
      19             : #include "libmesh/libmesh_config.h"
      20             : 
      21             : #ifdef LIBMESH_HAVE_TRIANGLE
      22             : 
      23             : // libmesh includes
      24             : #include "libmesh/mesh_triangle_wrapper.h"
      25             : #include "libmesh/boundary_info.h"
      26             : #include "libmesh/enum_elem_type.h"
      27             : #include "libmesh/face_tri3.h"
      28             : #include "libmesh/face_tri6.h"
      29             : #include "libmesh/point.h"
      30             : #include "libmesh/unstructured_mesh.h"
      31             : #include "libmesh/enum_to_string.h"
      32             : 
      33             : namespace libMesh
      34             : {
      35             : 
      36        1500 : void TriangleWrapper::init(TriangleWrapper::triangulateio & t)
      37             : {
      38        1500 :   t.pointlist                    = static_cast<REAL*>(nullptr);
      39        1500 :   t.pointattributelist           = static_cast<REAL*>(nullptr);
      40        1500 :   t.pointmarkerlist              = static_cast<int *>(nullptr);
      41        1500 :   t.numberofpoints               = 0 ;
      42        1500 :   t.numberofpointattributes      = 0 ;
      43             : 
      44        1500 :   t.trianglelist                 = static_cast<int *>(nullptr);
      45        1500 :   t.triangleattributelist        = static_cast<REAL*>(nullptr);
      46        1500 :   t.trianglearealist             = static_cast<REAL*>(nullptr);
      47        1500 :   t.neighborlist                 = static_cast<int *>(nullptr);
      48        1500 :   t.numberoftriangles            = 0;
      49        1500 :   t.numberofcorners              = 0;
      50        1500 :   t.numberoftriangleattributes   = 0;
      51             : 
      52        1500 :   t.segmentlist                  = static_cast<int *>(nullptr);
      53        1500 :   t.segmentmarkerlist            = static_cast<int *>(nullptr);
      54        1500 :   t.numberofsegments             = 0;
      55             : 
      56        1500 :   t.holelist                     = static_cast<REAL*>(nullptr);
      57        1500 :   t.numberofholes                = 0;
      58             : 
      59        1500 :   t.regionlist                   = static_cast<REAL*>(nullptr);
      60        1500 :   t.numberofregions              = 0;
      61             : 
      62        1500 :   t.edgelist                     = static_cast<int *>(nullptr);
      63        1500 :   t.edgemarkerlist               = static_cast<int *>(nullptr);
      64        1500 :   t.normlist                     = static_cast<REAL*>(nullptr);
      65        1500 :   t.numberofedges                = 0;
      66        1500 : }
      67             : 
      68             : 
      69             : 
      70             : 
      71             : 
      72             : 
      73        1500 : void TriangleWrapper::destroy(TriangleWrapper::triangulateio & t, TriangleWrapper::IO_Type io_type)
      74             : {
      75        1500 :   std::free (t.pointlist            );
      76        1500 :   std::free (t.pointattributelist   );
      77        1500 :   std::free (t.pointmarkerlist      );
      78        1500 :   std::free (t.trianglelist         );
      79        1500 :   std::free (t.triangleattributelist);
      80        1500 :   std::free (t.trianglearealist     );
      81        1500 :   std::free (t.neighborlist         );
      82        1500 :   std::free (t.segmentlist          );
      83        1500 :   std::free (t.segmentmarkerlist    );
      84             : 
      85             :   // Only attempt to free these when t was used as an input struct!
      86        1500 :   if (io_type==INPUT)
      87             :     {
      88         500 :       std::free (t.holelist  );
      89         500 :       std::free (t.regionlist);
      90             :     }
      91             : 
      92        1500 :   std::free (t.edgelist      );
      93        1500 :   std::free (t.edgemarkerlist);
      94        1500 :   std::free (t.normlist      );
      95             : 
      96             :   // Reset
      97             :   // TriangleWrapper::init(t);
      98        1500 : }
      99             : 
     100             : 
     101             : 
     102             : 
     103             : 
     104             : 
     105         500 : void TriangleWrapper::copy_tri_to_mesh(const triangulateio & triangle_data_input,
     106             :                                        UnstructuredMesh & mesh_output,
     107             :                                        const ElemType type,
     108             :                                        const triangulateio * voronoi)
     109             : {
     110             :   // Transfer the information into the LibMesh mesh.
     111         500 :   mesh_output.clear();
     112             : 
     113             :   // Make sure the new Mesh will be 2D
     114         500 :   mesh_output.set_mesh_dimension(2);
     115             : 
     116             :   // Node information
     117        4472 :   for (int i=0, c=0; c<triangle_data_input.numberofpoints; i+=2, ++c)
     118             :     {
     119             :       // Specify ID when adding point, otherwise, if this is DistributedMesh,
     120             :       // it might add points with a non-sequential numbering...
     121        3972 :       mesh_output.add_point( Point(triangle_data_input.pointlist[i],
     122        3972 :                                    triangle_data_input.pointlist[i+1]),
     123        3972 :                              /*id=*/c);
     124             :     }
     125             : 
     126             :   // Element information
     127        4512 :   for (int i=0; i<triangle_data_input.numberoftriangles; ++i)
     128             :     {
     129        4012 :       switch (type)
     130             :         {
     131        3964 :         case TRI3:
     132             :           {
     133        5946 :             Elem * elem = mesh_output.add_elem(Elem::build(TRI3));
     134             : 
     135       15856 :             for (unsigned int n=0; n<3; ++n)
     136       11892 :               elem->set_node(n, mesh_output.node_ptr(triangle_data_input.trianglelist[i*3 + n]));
     137             : 
     138             :             // use the first attribute to set the subdomain ID
     139        3964 :             if (triangle_data_input.triangleattributelist)
     140           0 :               elem->subdomain_id() =
     141           0 :                 std::round(triangle_data_input.
     142           0 :                            triangleattributelist[i * triangle_data_input.numberoftriangleattributes]);
     143        1982 :             break;
     144             :           }
     145             : 
     146          48 :         case TRI6:
     147             :           {
     148          48 :             Elem * elem = mesh_output.add_elem(Elem::build(TRI6));
     149             : 
     150             :             // Triangle number TRI6 nodes in a different way to libMesh
     151          48 :             elem->set_node(0, mesh_output.node_ptr(triangle_data_input.trianglelist[i*6 + 0]));
     152          48 :             elem->set_node(1, mesh_output.node_ptr(triangle_data_input.trianglelist[i*6 + 1]));
     153          48 :             elem->set_node(2, mesh_output.node_ptr(triangle_data_input.trianglelist[i*6 + 2]));
     154          48 :             elem->set_node(3, mesh_output.node_ptr(triangle_data_input.trianglelist[i*6 + 5]));
     155          48 :             elem->set_node(4, mesh_output.node_ptr(triangle_data_input.trianglelist[i*6 + 3]));
     156          48 :             elem->set_node(5, mesh_output.node_ptr(triangle_data_input.trianglelist[i*6 + 4]));
     157             : 
     158             :             // use the first attribute to set the subdomain ID
     159          48 :             if (triangle_data_input.triangleattributelist)
     160           0 :               elem->subdomain_id() =
     161           0 :                 std::round(triangle_data_input.
     162           0 :                            triangleattributelist[i * triangle_data_input.numberoftriangleattributes]);
     163          24 :             break;
     164             :           }
     165             : 
     166           0 :         default:
     167           0 :           libmesh_error_msg("ERROR: Unrecognized triangular element type == " << Utility::enum_to_string(type));
     168             :         }
     169             :     }
     170             : 
     171             :   // Note: If the input mesh was a parallel one, calling
     172             :   // prepare_for_use() now will re-parallelize it by a call to
     173             :   // delete_remote_elements()... We do not actually want to
     174             :   // reparallelize it here though: the triangulate() function may
     175             :   // still do some Mesh smoothing.  The main thing needed (for
     176             :   // smoothing) is the neighbor information, so let's just find
     177             :   // neighbors...
     178             :   //mesh_output.prepare_for_use(/*skip_renumber =*/false);
     179         500 :   mesh_output.find_neighbors();
     180             : 
     181             :   // set boundary info
     182         500 :   if (voronoi && triangle_data_input.edgemarkerlist)
     183             :   {
     184           0 :     BoundaryInfo & boundary_info = mesh_output.get_boundary_info();
     185           0 :     for (int e=0; e<triangle_data_input.numberofedges; ++e)
     186             :     {
     187           0 :       if (triangle_data_input.edgemarkerlist[e] != 0)
     188             :       {
     189           0 :         int p1 = triangle_data_input.edgelist[e + e];
     190           0 :         int p2 = triangle_data_input.edgelist[e + e + 1];
     191           0 :         int elem_id = voronoi->edgelist[e + e];
     192             :         unsigned short int s;
     193           0 :         if (p1 == triangle_data_input.trianglelist[elem_id*3] &&
     194           0 :             p2 == triangle_data_input.trianglelist[elem_id*3 + 1])
     195           0 :           s = 0;
     196           0 :         else if (p1 == triangle_data_input.trianglelist[elem_id*3 + 1] &&
     197           0 :                  p2 == triangle_data_input.trianglelist[elem_id*3 + 2])
     198           0 :           s = 1;
     199           0 :         else if (p1 == triangle_data_input.trianglelist[elem_id*3 + 2] &&
     200           0 :                  p2 == triangle_data_input.trianglelist[elem_id*3])
     201           0 :           s = 2;
     202             :         else
     203           0 :           libmesh_error_msg("ERROR: finding element errors for boundary edges.");
     204             : 
     205           0 :         boundary_info.add_side(elem_id, s, triangle_data_input.edgemarkerlist[e]);
     206             :       }
     207             :     }
     208             :   }
     209         500 : }
     210             : 
     211             : 
     212             : }
     213             : 
     214             : #endif // LIBMESH_HAVE_TRIANGLE