mesh_tetgen_support.c

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

// $Id: mesh_tetgen_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 "libmesh_config.h"
#ifdef HAVE_TETGEN


// C++ includes
#include <sstream>
#include <set>

// Local includes
#include "cell_tet4.h"
#include "face_tri3.h"
#include "unstructured_mesh.h"
#include "mesh_tetgen_support.h"


//----------------------------------------------------------------------
// TetGenMeshInterface functions
TetGenWrapper::TetGenWrapper()
{
  tetgen_output = new tetgenio;

  this->tetgen_data.mesh_dim                = 3;
  this->tetgen_data.numberofpointattributes = 0;
  this->tetgen_data.firstnumber             = 0;
}



TetGenWrapper::~TetGenWrapper()
{
  delete tetgen_output;
}



void TetGenWrapper::set_node(const int i, const REAL x, const REAL y, const REAL z)
{
  int index = i*3;
  tetgen_data.pointlist[index++] = x;
  tetgen_data.pointlist[index++] = y;
  tetgen_data.pointlist[index++] = z;
}



void TetGenWrapper::set_hole(const int i, const REAL x, const REAL y, const REAL z)
{
  int index = i*3;
  tetgen_data.holelist[index++] = x;
  tetgen_data.holelist[index++] = y;
  tetgen_data.holelist[index++] = z;
}



void TetGenWrapper::set_numberofpoints(const int i)
{
  tetgen_data.numberofpoints = i;
}



void TetGenWrapper::get_output_node(const int i, REAL& x, REAL& y, REAL& z)
{
  x = tetgen_output->pointlist[3*i];
  y = tetgen_output->pointlist[3*i+1];
  z = tetgen_output->pointlist[3*i+2];
}



int TetGenWrapper::get_numberoftetrahedra()
{
  return tetgen_output->numberoftetrahedra;
}



int TetGenWrapper::get_numberoftrifaces()
{
  return tetgen_output->numberoftrifaces;
}



int TetGenWrapper::get_numberofpoints()
{
  return tetgen_output->numberofpoints;
}



int TetGenWrapper::get_element_node(const int i, const int j)
{
  return tetgen_output->tetrahedronlist[i*4+j];
}



int TetGenWrapper::get_triface_node(const int i, const int j)
{
  return tetgen_output->trifacelist[i*3+j];
}



REAL TetGenWrapper::get_element_attribute(const int i)
{
  libmesh_assert(tetgen_output->numberoftetrahedronattributes>0);
  return tetgen_output->tetrahedronattributelist[tetgen_output->numberoftetrahedronattributes*i];
}



void TetGenWrapper::set_pointlist(const int numofpoints)
{
  this->set_numberofpoints(numofpoints);
  this->tetgen_data.pointlist = new REAL[tetgen_data.numberofpoints * 3];
}



void TetGenWrapper::set_switches(const std::string& s)
{
  // Copy the string to a temporary buffer for passing to the C API
  char buffer[256];
  libmesh_assert (s.size() < sizeof(buffer)-1);
  buffer[ s.copy( buffer , sizeof( buffer ) - 1 ) ] = '\0' ;
  
  if (!tetgen_be.parse_commandline(buffer)) 
    std::cout << "TetGen replies: Wrong switches!" << std::endl;
}



void TetGenWrapper::run_tetgen()
{
  // Call tetrahedralize from the TetGen library.
  tetrahedralize(&tetgen_be, &tetgen_data, tetgen_output);
}



void TetGenWrapper::set_numberoffacets(const int i)
{ 
  this->tetgen_data.numberoffacets = i;
}



void TetGenWrapper::set_numberofholes(const int i)
{ 
  this->tetgen_data.numberofholes = i;
}



void TetGenWrapper::set_numberofregions(const int i)
{
  this->tetgen_data.numberofregions = i;
}



void TetGenWrapper::set_facetlist(const int numoffacets, const int numofholes)
{
  set_numberoffacets(numoffacets);
  set_numberofholes(numofholes);
  this->tetgen_data.facetlist = new tetgenio::facet[this->tetgen_data.numberoffacets];
  for (int i=0; i<numoffacets; i++)
    this->tetgen_data.init(&(this->tetgen_data.facetlist[i]));
  this->tetgen_data.holelist = new REAL[this->tetgen_data.numberofholes * 3];
}



void TetGenWrapper::set_regionlist(const int numofregions)
{
  set_numberofregions(numofregions);
  this->tetgen_data.regionlist = new REAL[this->tetgen_data.numberofregions * 5];
}



void TetGenWrapper::set_facet_numberofpolygons(const int i, const int num)
{
  this->tetgen_data.facetlist[i].numberofpolygons = num;
}



void TetGenWrapper::set_facet_numberofholes(const int i, const int num)
{
  this->tetgen_data.facetlist[i].numberofholes = num;
}




void TetGenWrapper::set_facet_polygonlist(const int i, const int numofpolygons)
{
  set_facet_numberofpolygons(i, numofpolygons);
  set_facet_numberofholes(i, 0);
  this->tetgen_data.facetlist[i].polygonlist = new tetgenio::polygon[numofpolygons];
  for (int j=0; j<this->tetgen_data.facetlist[i].numberofpolygons; j++)
    this->tetgen_data.init(&(this->tetgen_data.facetlist[i].polygonlist[j]));
}



void TetGenWrapper::set_polygon_numberofvertices(const int i, const int j, const int num)
{
  this->tetgen_data.facetlist[i].polygonlist[j].numberofvertices = num;
}



void TetGenWrapper::set_polygon_vertexlist(const int i, const int j, const int numofvertices)
{
  set_polygon_numberofvertices(i, j, numofvertices);
  this->tetgen_data.facetlist[i].polygonlist[j].vertexlist = new int[numofvertices];
}




void TetGenWrapper::set_vertex(const int i, const int j, const int k, const int nodeindex)
{
  this->tetgen_data.facetlist[i].polygonlist[j].vertexlist[k] = nodeindex;
}



void TetGenWrapper::set_region(const int i, const REAL x, const REAL y, const REAL z,
			       const REAL attribute, const REAL vol_constraint)
{
  int index = i*5;
  tetgen_data.regionlist[index++] = x;
  tetgen_data.regionlist[index++] = y;
  tetgen_data.regionlist[index++] = z;
  tetgen_data.regionlist[index++] = attribute;
  tetgen_data.regionlist[index++] = vol_constraint;
}


// class type TetGen_access is cast to TetGenWrapper.
typedef TetGenWrapper TetGen_access;






//----------------------------------------------------------------------
// TetGenMeshInterface class members
TetGenMeshInterface::TetGenMeshInterface (UnstructuredMesh& mesh) :
  _mesh         (mesh)
{
}



void TetGenMeshInterface::triangulate_pointset ()   
{
  // class tetgen_wrapper allows library access on a basic level:
  TetGen_access tetgen_wrapper; 

  // fill input structure with point set data:
  tetgen_wrapper.set_pointlist( this->_mesh.n_nodes() );

  {
    int index = 0;
    MeshBase::node_iterator it  = this->_mesh.nodes_begin();