mesh_triangle_support.c

一个用来实现偏微分方程中网格的计算库

// $Id: mesh_triangle_support.C 2789 2008-04-13 02:24:40Z roystgnr $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



#include "mesh_triangle_support.h"
#include "unstructured_mesh.h"
#include "face_tri3.h"
#include "face_tri6.h"
#include "mesh_generation.h"
#include "mesh_smoother_laplace.h"
#include "boundary_info.h"


#ifdef HAVE_TRIANGLE

// Triangulates a 2D rectangular region with or without holes
void MeshTools::Generation::build_delaunay_square(UnstructuredMesh& mesh,
						  const unsigned int nx, // num. of elements in x-dir
						  const unsigned int ny, // num. of elements in y-dir
						  const Real xmin, const Real xmax,
						  const Real ymin, const Real ymax,
						  const ElemType type,
						  const std::vector<TriangleInterface::Hole*>* holes)
{
  // Check for existing nodes and compatible dimension.
  libmesh_assert (mesh.mesh_dimension() == 2);
  libmesh_assert (nx >= 1); // need at least 1 element in x-direction
  libmesh_assert (ny >= 1); // need at least 1 element in y-direction
  libmesh_assert (xmin < xmax);
  libmesh_assert (ymin < ymax);

  // Clear out any data which may have been in the Mesh
  mesh.clear();
  
  // The x and y spacing between boundary points
  const Real delta_x = (xmax-xmin) / static_cast<Real>(nx);
  const Real delta_y = (ymax-ymin) / static_cast<Real>(ny);  

  // Bottom
  for (unsigned int p=0; p<=nx; ++p)
    mesh.add_point(Point(xmin + p*delta_x, ymin));
  
  // Right side
  for (unsigned int p=1; p<ny; ++p)
    mesh.add_point(Point(xmax, ymin + p*delta_y)); 

  // Top
  for (unsigned int p=0; p<=nx; ++p)
    mesh.add_point(Point(xmax - p*delta_x, ymax));

  // Left side
  for (unsigned int p=1; p<ny; ++p)
    mesh.add_point(Point(xmin,  ymax - p*delta_y)); 
		   
  // Be sure we added as many points as we thought we did
  libmesh_assert (mesh.n_nodes() == 2*(nx+ny));

  // Construct the Triangle Interface object
  TriangleInterface t(mesh);

  // Set custom variables for the triangulation
  t.desired_area()       = 0.5 * (xmax-xmin)*(ymax-ymin) / static_cast<Real>(nx*ny);
  t.triangulation_type() = TriangleInterface::PSLG;
  t.elem_type()          = type;

  if (holes != NULL)
    t.attach_hole_list(holes);

  // Triangulate!
  t.triangulate();

  // The mesh is now generated, but we still need to mark the boundaries
  // to be consistent with the other build_square routines.  Note that only
  // the external boundaries of the square are given boundary ids, the holes
  // are not numbered.
  MeshBase::element_iterator       el     = mesh.elements_begin();
  const MeshBase::element_iterator end_el = mesh.elements_end();
  for ( ; el != end_el; ++el)
    {
      const Elem* elem = *el;
      
      for (unsigned int s=0; s<elem->n_sides(); s++)
	if (elem->neighbor(s) == NULL)
	  {
	    AutoPtr<Elem> side (elem->build_side(s));

	    // Check the location of the side's midpoint.  Since
	    // the square has straight sides, the midpoint is not
	    // on the corner and thus it is uniquely on one of the
	    // sides.
	    Point side_midpoint= 0.5*( (*side->get_node(0)) + (*side->get_node(1)) );

	    // The boundary ids are set following the same convention as Quad4 sides
	    // bottom = 0
	    // right  = 1
	    // top = 2
	    // left = 3
	    // hole = 4
	    short int bc_id=4;

	    // bottom
	    if      (std::fabs(side_midpoint(1) - ymin) < TOLERANCE)
	      bc_id=0;

	    // right
	    else if (std::fabs(side_midpoint(0) - xmax) < TOLERANCE)
	      bc_id=1;

	    // top
	    else if (std::fabs(side_midpoint(1) - ymax) < TOLERANCE)
	      bc_id=2;

	    // left
	    else if (std::fabs(side_midpoint(0) - xmin) < TOLERANCE)
	      bc_id=3;

	    // If the point is not on any of the external boundaries, it
	    // is on one of the holes....
	    
	    // Finally, add this element's information to the boundary info object.
	    mesh.boundary_info->add_side(elem->id(), s, bc_id);
	  }
    }
  
}






// Init helper routine defined in the Triangle namespace
void Triangle::init(Triangle::triangulateio& t)
{
  t.pointlist                    = static_cast<REAL*>(NULL);
  t.pointattributelist           = static_cast<REAL*>(NULL);
  t.pointmarkerlist              = static_cast<int* >(NULL);
  t.numberofpoints               = 0 ;
  t.numberofpointattributes      = 0 ;                                   

  t.trianglelist                 = static_cast<int* >(NULL);
  t.triangleattributelist        = static_cast<REAL*>(NULL);
  t.trianglearealist             = static_cast<REAL*>(NULL);
  t.neighborlist                 = static_cast<int* >(NULL);
  t.numberoftriangles            = 0;
  t.numberofcorners              = 0;
  t.numberoftriangleattributes   = 0;
  
  t.segmentlist                  = static_cast<int* >(NULL);
  t.segmentmarkerlist            = static_cast<int* >(NULL);
  t.numberofsegments             = 0;

  t.holelist                     = static_cast<REAL*>(NULL);
  t.numberofholes                = 0;

  t.regionlist                   = static_cast<REAL*>(NULL);
  t.numberofregions              = 0;
  
  t.edgelist                     = static_cast<int* >(NULL);
  t.edgemarkerlist               = static_cast<int* >(NULL);
  t.normlist                     = static_cast<REAL*>(NULL);
  t.numberofedges                = 0;
}






// Destroy helper routine defined in the Triangle namespace
void Triangle::destroy(Triangle::triangulateio& t, Triangle::IO_Type io_type)
{
  std::free (t.pointlist                    );
  std::free (t.pointattributelist           );
  std::free (t.pointmarkerlist              );
  std::free (t.trianglelist                 );
  std::free (t.triangleattributelist        );
  std::free (t.trianglearealist             );
  std::free (t.neighborlist                 );
  std::free (t.segmentlist                  );
  std::free (t.segmentmarkerlist            );

  // Only attempt to free these when t was used as an input struct!
  if (io_type==INPUT)
    {
      std::free (t.holelist                     );
      std::free (t.regionlist                   );
    }
  
  std::free (t.edgelist                     );
  std::free (t.edgemarkerlist               );
  std::free (t.normlist                     );

  // Reset
  // Triangle::init(t);
}






void Triangle::copy_tri_to_mesh(const triangulateio& triangle_data_input,
				UnstructuredMesh& mesh_output,
				const ElemType type)
{
  // Transfer the information into the LibMesh mesh.
  mesh_output.clear();
  
  // Node information
  for (int i=0, c=0; c<triangle_data_input.numberofpoints; i+=2, ++c)
    {
      mesh_output.add_point( Point(triangle_data_input.pointlist[i],
				   triangle_data_input.pointlist[i+1]) );
    }

  // Element information
  for (int i=0; i<triangle_data_input.numberoftriangles; ++i)
    {
      switch (type)
	{
        case TRI3:
	  {
	    Elem* elem = mesh_output.add_elem (new Tri3);

	    for (unsigned int n=0; n<3; ++n)
	      elem->set_node(n) = mesh_output.node_ptr(triangle_data_input.trianglelist[i*3 + n]);

	    break;
	  }

        case TRI6:
	  {
	    Elem* elem = mesh_output.add_elem (new Tri6);

	    // Triangle number TRI6 nodes in a different way to libMesh
	    elem->set_node(0) = mesh_output.node_ptr(triangle_data_input.trianglelist[i*6 + 0]);
	    elem->set_node(1) = mesh_output.node_ptr(triangle_data_input.trianglelist[i*6 + 1]);
	    elem->set_node(2) = mesh_output.node_ptr(triangle_data_input.trianglelist[i*6 + 2]);
	    elem->set_node(3) = mesh_output.node_ptr(triangle_data_input.trianglelist[i*6 + 5]);
	    elem->set_node(4) = mesh_output.node_ptr(triangle_data_input.trianglelist[i*6 + 3]);
	    elem->set_node(5) = mesh_output.node_ptr(triangle_data_input.trianglelist[i*6 + 4]);

	    break;
	  }

        default:
	  {
	    std::cerr << "ERROR: Unrecognized triangular element type." << std::endl;
	    libmesh_error();
	  }
	}
    }

  // Prepare mesh for usage.
  mesh_output.prepare_for_use();
}






// Function definitions for the TriangleInterface class
void TriangleInterface::triangulate()
{
  // Will the triangulation have holes?
  const bool have_holes = ((_holes != NULL) && (!_holes->empty()));
  
  // Currently, we don't support specifying user-defined segments
  // *and* having holes