tetgen.cpp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

// ORIG-DATE:     Juin 2008
// -*- Mode : c++ -*-
//
// SUMMARY  : liaison medit freefem++ : popen  
// USAGE    : LGPL      
// ORG      : LJLL Universite Pierre et Marie Curie, Paris,  FRANCE 
// AUTHOR   : Jacques Morice
// E-MAIL   : [e-mail]
//

/* 
 This file is part of Freefem++
 
 Freefem++ 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.
 
 Freefem++  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 Freefem++; if not, write to the Free Software
 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

 Thank to the ARN ()  FF2A3 grant
 ref:ANR-07-CIS7-002-01 
 */
// $Id: tetgen.cpp,v 1.13 2009/02/03 14:03:02 morice Exp $

#include  <iostream>
#include  <cfloat>
#include <cmath>
#include <complex>
using namespace std;
#include "error.hpp"
#include "AFunction.hpp"
using namespace std;  
#include "rgraph.hpp"
#include "RNM.hpp"
#include "fem.hpp"


#include "FESpacen.hpp" 
#include "FESpace.hpp" 

#include "MatriceCreuse_tpl.hpp"
#include "MeshPoint.hpp"
#include "Operator.hpp" 
#include "lex.hpp"

#include "lgfem.hpp"
#include "lgmesh3.hpp"
#include "lgsolver.hpp"
#include "problem.hpp"
//#include "LayerMesh.hpp"
//#include "TransfoMesh_v2.hpp"
#include "msh3.hpp"
#include "tetgen.h"
//#include "GQuadTree.hpp"

#include <set>
#include <vector>
#include <fstream>


using namespace  Fem2D;

// subroutine use for tetegen call

void mesh3_tetgenio_out(const tetgenio &out, Mesh3 & Th3);

void mesh3_tetgenio_out(const tetgenio &out, const int & label_tet, Mesh3 & Th3);

void mesh3_tetgenio_out(const tetgenio &out, const int & label_tet, const int & label_face, Mesh3 & Th3);

Mesh3 * mesh3_tetgenio_out(const tetgenio &out);

Mesh3 *mesh3_tetgenio_out(const tetgenio &out, const int & label_tet);

Mesh3 * mesh3_tetgenio_out(const tetgenio &out, const int & label_tet, const int & label_face);


Mesh3 * Convexhull_3Dpoints( char* switch_tetgen, const int &nv_t, const double *Xcoord, const double *Ycoord, const double *Zcoord, const int &label_tet);

Mesh3 * RemplissageSurf3D_tetgen( char* switch_tetgen, const Mesh3 & Th3, const int & label_tet);

Mesh3 * RemplissageSurf3D_tetgen_new( char* switch_tetgen, const Mesh3 & Th3, const int & label_tet,
				  const int &nbhole, const double *tabhole, 
				  const int & nbregion, const double *tabregion, 
				  const int &nbfacecl, const double *tabfacecl);

Mesh3 * Transfo_Mesh2_tetgen( const double &precis_mesh, char* switch_tetgen, const Mesh & Th2, const double *tab_XX, const double *tab_YY, const double *tab_ZZ, 
		int &border_only, int &recollement_border, int &point_confondus_ok, 
		const int &label_tet,const map<int, int> &maptri );

Mesh3 *Transfo_Mesh2_tetgen_new(const double &precis_mesh, char *switch_tetgen,const Mesh & Th2, const double *tab_XX, const double *tab_YY, const double *tab_ZZ, 
				 int &border_only, int &recollement_border, int &point_confondus_ok, 
				 const int &label_tet, const map<int, int> &maptri, 
				 const int &nbhole, const double *tabhole, 
				 const int & nbregion, const double *tabregion, 
				 const int &nbfacecl, const double *tabfacecl);

Mesh3 *  ReconstructionRefine_tetgen(char *switch_tetgen,const Mesh3 & Th3, 
				     const int &nbhole, const double *tabhole, 
				     const int & nbregion, const double *tabregion, 
		                     const int &nbfacecl, const double *tabfacecl, const double *tsizevol);


class Build2D3D_Op : public E_F0mps 
{
public:
  Expression eTh;
  Expression xx,yy,zz;
  static const int n_name_param =13; // 
  static basicAC_F0::name_and_type name_param[] ;
  Expression nargs[n_name_param];
  KN_<long>  arg(int i,Stack stack,KN_<long> a ) const
  { return nargs[i] ? GetAny<KN_<long> >( (*nargs[i])(stack) ): a;}
  KN_<double>  arg(int i,Stack stack,KN_<double> a ) const
  { return nargs[i] ? GetAny<KN_<double> >( (*nargs[i])(stack) ): a;}
  double  arg(int i,Stack stack,double a) const{ return nargs[i] ? GetAny< double >( (*nargs[i])(stack) ): a;}
  int  arg(int i,Stack stack, int a) const{ return nargs[i] ? GetAny< int >( (*nargs[i])(stack) ): a;}
  string*  arg(int i,Stack stack, string* a) const{ return nargs[i] ? GetAny< string* >( (*nargs[i])(stack) ): a;}
public:
  Build2D3D_Op(const basicAC_F0 &  args,Expression tth) 
    : eTh(tth),xx(0),yy(0),zz(0)
  {
    if(verbosity) cout << "construction par BuilLayeMesh_Op" << endl;
    args.SetNameParam(n_name_param,name_param,nargs);
    const E_Array * a1=0 ;
    if(nargs[0])  a1  = dynamic_cast<const E_Array *>(nargs[0]);
    int err =0;
   
    if(a1) {
      if(a1->size() !=3) 
      CompileError("Build2D3D (Th,transfo=[X,Y,Z],) ");
      xx=to<double>( (*a1)[0]);
      yy=to<double>( (*a1)[1]);
      zz=to<double>( (*a1)[2]);
    }    
  } 
  
  AnyType operator()(Stack stack)  const ;
};


basicAC_F0::name_and_type Build2D3D_Op::name_param[]= {
  {  "transfo", &typeid(E_Array)},
  {  "switch", &typeid(string*)},
  {  "reftet", &typeid(long)}, 
  {  "refface", &typeid(KN_<long>)},
  {  "facemerge", &typeid(long)},
  {  "ptmerge", &typeid(double)},
  // nouvelle variable
  {  "nbofholes", &typeid(long)},
  {  "holelist", &typeid(KN_<double>)},
  {  "nbofregions", &typeid(long)},
  {  "regionlist", &typeid(KN_<double>)},
  {  "nboffacetcl", &typeid(long)},
  {  "facetcl", &typeid(KN_<double>)},
  // mesure mesh
  {  "mesuremesh", &typeid(long)}
};

class  Build2D3D : public OneOperator { public:  
    Build2D3D() : OneOperator(atype<pmesh3>(),atype<pmesh>()) {}
  
  E_F0 * code(const basicAC_F0 & args) const 
  { 
	return  new Build2D3D_Op( args,t[0]->CastTo(args[0]) ); 
  }
};

AnyType Build2D3D_Op::operator()(Stack stack)  const 
{
  MeshPoint *mp(MeshPointStack(stack)) , mps=*mp;
  Mesh * pTh= GetAny<Mesh *>((*eTh)(stack));
  ffassert( pTh );
  Mesh &Th=*pTh;
  Mesh *m= pTh;   // question a quoi sert *m ??
  int nbv=Th.nv; // nombre de sommet 

  int nbt=Th.nt; // nombre de triangles
  int neb=Th.neb; // nombre d'aretes fontiere
  if(verbosity) cout << " Vertex Triangle Border " << nbv<< "  "<< nbt << " " << neb << endl; 

  if(verbosity >1) cout <<" ======================= " << endl;
  if(verbosity >1) cout <<" == Build2D_3D_Op==" << endl;
 
  KN<long> zzempty;
  string  stringempty = string("pqaAAYCQ");
  string* switch_tet= (arg(1,stack,&stringempty));  
  int label_tet(arg(2,stack,0));  
  KN<long> nrf (arg(3,stack,zzempty));
  int point_confondus_ok(arg(4,stack,0));
  double precis_mesh(arg(5,stack,-1));

  // new parameters
  KN<double> zdzempty;
  int nbhole (arg(6,stack,0));
  KN<double> tabhole (arg(7,stack,zdzempty));
  int nbregion (arg(8,stack,0));
  KN<double> tabregion (arg(9,stack,zdzempty));
  int nbfacecl (arg(10,stack,0));
  KN<double> tabfacecl (arg(11,stack,zdzempty));
  
  // mesuremesh parameters
  int mesureM(arg(13,stack,1));
  int surface_orientation=1; 
  if( mesureM <0 ){
    surface_orientation=-1;
  }

  // assertion au niveau de la taille
  ffassert( tabhole.N()   == 3*nbhole);
  ffassert( tabregion.N() == 5*nbregion);
  ffassert( tabfacecl.N() == 2*nbfacecl);

  //====================================
  //  How to change string* into char* 
  //====================================
  size_t size_switch_tet = switch_tet->size()+1;
  char* switch_tetgen =new char[size_switch_tet];
  strncpy(switch_tetgen, switch_tet->c_str(), size_switch_tet); 

  cout << "switch_tetgen=" << switch_tetgen << endl;

  ffassert( nrf.N() %2 ==0);
  
  map<int,int> mapf;
  for(int i=0;i<nrf.N();i+=2)
    {
      if(nrf[i] != nrf[i+1]){
	mapf[nrf[i]]=nrf[i+1];
      }
    }
  
  map<int, int> mapfme;  
  
  Transfo_Mesh2_map_face( Th, mapfme );  
  
  // Map utilisateur
  map< int, int > :: iterator imap;
  for( int ii=0; ii < nrf.N(); ii+=2){
	imap = mapfme.find(nrf[ii]);
	if( imap != mapfme.end()){
		imap -> second = nrf[ii+1];
	}
  }  
  
  //KN<double> txx(nbv), tyy(nbv), tzz(nbv);
  //KN<int> takemesh(nbv);
  double *txx=new double[nbv];
  double *tyy=new double[nbv];
  double *tzz=new double[nbv];
  int *takemesh=new int[nbv];

  MeshPoint *mp3(MeshPointStack(stack)); 
  
  for(int ii=0; ii<nbv; ii++) 
      takemesh[ii]=0;  

  Mesh &rTh = Th;
  for (int it=0; it<nbt; ++it){
    for( int iv=0; iv<3; ++iv){
      int i=Th(it,iv);
      if(takemesh[i]==0){
	mp3->setP(&Th,it,iv);
	if(xx){ 
	  txx[i]=GetAny<double>((*xx)(stack));
	}
	if(yy){ 
	  tyy[i]=GetAny<double>((*yy)(stack));
	}
	if(zz){ 
	  tzz[i]=GetAny<double>((*zz)(stack));
	}
	takemesh[i] = takemesh[i]+1;
      }
    }
  }

  delete [] takemesh;
  int border_only = 0;
  int recollement_border=1;
  /*
    Mesh3 *Th3=Transfo_Mesh2_tetgen( precis_mesh, switch_tetgen, Th, txx, tyy, tzz, border_only, 
				   recollement_border, point_confondus_ok, label_tet, mapfme);  
  */

  Mesh3 *Th3_tmp = MoveMesh2_func( precis_mesh, Th, txx, tyy, tzz, border_only, recollement_border, point_confondus_ok);

  /* delete array */
  delete [] txx;
  delete [] tyy;
  delete [] tzz;

  /* check orientation of the mesh and flip if necessary*/ 
  Th3_tmp->flipSurfaceMesh3(surface_orientation);
 
  cout << "check :: orientation des surfaces" << endl;
  Th3_tmp->BuildSurfaceAdj();
  cout << "fin check :: orientation des surfaces" << endl;

  /* set label of surface Th3_tmp */
  for(int ii=0; ii< Th3_tmp->nbe; ii++)
    {
      const Triangle3 & K(Th3_tmp->be(ii)); 
      int iv[3];
      int lab;
     
      iv[0] = Th3_tmp->operator()(K[0]);
      iv[1] = Th3_tmp->operator()(K[1]);
      iv[2] = Th3_tmp->operator()(K[2]);
		
      map< int, int>:: const_iterator imap;
      imap = mapfme.find(K.lab); 
		
      if(imap!= mapfme.end()){
	lab=imap->second;	
      } 
      else{
	lab=K.lab;
      }
		
      Th3_tmp->be(ii).set( Th3_tmp->vertices, iv, lab ) ;
    }
  /* mesh domains with tetgen */ 
  Mesh3 *Th3 = RemplissageSurf3D_tetgen_new( switch_tetgen, *Th3_tmp, label_tet,
					     nbhole, tabhole, nbregion, tabregion, 
					     nbfacecl, tabfacecl);


  /*
  Mesh3 *Th3=Transfo_Mesh2_tetgen_new( precis_mesh, switch_tetgen, Th, txx, tyy, tzz, border_only, 
				   recollement_border, point_confondus_ok, label_tet, mapfme, 
				   nbhole, tabhole, nbregion, tabregion, nbfacecl,tabfacecl);

  */

  delete Th3_tmp;
  
  Th3->BuildBound();
  Th3->BuildAdj();
  Th3->Buildbnormalv();  
  Th3->BuildjElementConteningVertex();
  Th3->BuildGTree();
  //Th3->decrement();    
  Add2StackOfPtr2FreeRC(stack,Th3);

  delete [] switch_tetgen;
		
  *mp=mps;
  cout << "FreeFem++: End check mesh given by tetgen" << endl;  
  return Th3;
}

// Fonction pour tetgen
// new parameter
void mesh3_tetgenio_out(const tetgenio &out, Mesh3 & Th3)
{ 
  int i;

// All indices start from 1.
  if(out.firstnumber != 1){
    cout << " probleme ???" << endl;
    exit(1);
  }   
  
  if(out.numberoffacets !=0){
    cout << "tetgen: faces non triangulaire" << endl;
    exit(1);
  }
  
  if(out.numberofcorners !=4){
    cout << "tetgen: element subparametric of order 2" <<endl;
    exit(1);
  }

  cout << "Th3 :: Vertex Element Border :: " << out.numberofpoints << " " <<out.numberoftetrahedra  << " " << out.numberoftrifaces << endl;
  Th3.set(out.numberofpoints, out.numberoftetrahedra, out.numberoftrifaces);

  // new parameter
  if( out.numberoftetrahedronattributes !=  1){
    cout << "out.numberoftetrahedronattributes" << out.numberoftetrahedronattributes << endl;
  }
   
  i=0;
  for(int nnv=0; nnv < Th3.nv; nnv++){
    Th3.vertices[nnv].x=out.pointlist[i];
    Th3.vertices[nnv].y=out.pointlist[i+1];
    Th3.vertices[nnv].z=out.pointlist[i+2];       
    Th3.vertices[nnv].lab=out.pointmarkerlist[nnv];
    i=i+3;    
  }
    
  i=0;
  for(int nnt=0; nnt < Th3.nt; nnt++){
    int iv[4],lab;
    iv[0] = out.tetrahedronlist[i]-1;
    iv[1] = out.tetrahedronlist[i+1]-1;
    iv[2] = out.tetrahedronlist[i+2]-1;
    iv[3] = out.tetrahedronlist[i+3]-1;
    //lab   = label_tet;
    for(int jj=0; jj<4; jj++){
      assert( iv[jj] >=0 && iv[jj] <Th3.nv );
    }
    //cout << "nnt= " <<  nnt << " " << lab  << " " << out.tetrahedronattributelist[nnt] << endl;
    lab =  out.tetrahedronattributelist[nnt];
    //cout << "nnt= " <<  lab  << " " << out.tetrahedronattributelist[nnt] << endl;
      
    Th3.elements[nnt].set( Th3.vertices, iv, lab);
    i=i+4;
  }
 
  for(int ibe=0; ibe < Th3.nbe; ibe++){
    int iv[3];