mesh3dn.hpp

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

// ORIG-DATE:     Dec 2007
// -*- Mode : c++ -*-
//
// SUMMARY  :  Model  mesh 1d   
// USAGE    : LGPL      
// ORG      : LJLL Universite Pierre et Marie Curi, Paris,  FRANCE 
// AUTHOR   : Frederic Hecht
// 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 
 */


#ifndef MESH3DN_HPP_
#define MESH3DN_HPP_



// definition R
#include <cstdlib>
namespace Fem2D {

#include "R3.hpp"
}
 
 using namespace ::std;
#include "GenericMesh.hpp" 

namespace Fem2D {
  
typedef GenericVertex<R3> Vertex3;

struct DataTriangle3  {
  static const int NbOfVertices =3;
  static const int NbOfEdges =3;
  static const int NbOfFaces =1;
  static const int NT =0;
  static const int NbOfAdjElem =NbOfVertices;
  static const int NbOfVertexOnHyperFace =NbOfVertices-1;
  typedef Vertex3 V;
  typedef  V::Rd Rd ;
  typedef R2 RdHat;
  typedef R1 RdHatBord;
  static RdHat PBord(const int * nvb,const RdHatBord &P)  { 
  return RdHat::KHat[nvb[0]]*(1-P.x)+R2::KHat[nvb[1]]*(P.x) ;}  

  static R mesure(  V *  pv[NbOfVertices]) {    
    return (R3(*pv[0],*pv[1])^R3(*pv[0],*pv[2])).norme()*0.5;
  }
};

struct DataTet  {
  static const int NbOfVertices =4;
  static const int NbOfEdges =6;
  static const int NbOfFaces =4;
  static const int NT =1;
  static const int NbOfAdjElem =NbOfVertices;
  static const int NbOfVertexOnHyperFace =NbOfVertices-1;
  typedef Vertex3 V;
  typedef  V::Rd Rd ;
  static R mesure(  V *  pv[NbOfVertices]) 
  {    
    R3 AB(*pv[0],*pv[1]);
    R3 AC(*pv[0],*pv[2]);
    R3 AD(*pv[0],*pv[3]);
    return det(AB,AC,AD)/6.;
  }
  static const int (* const nvface)[3];// = nvfaceTet;
  static const int (* const nvedge)[2];//  = nvedgeTet;
  typedef R3 RdHat;
  typedef R2 RdHatBord;
  static RdHat PBord(const int * nvb,const RdHatBord& P)  { 
 //     cout << "PBORD : " << nvb[0] << " " <<  nvb[1] <<  nvb[2] << " " << P<< " -> " <<  RdHat::KHat[nvb[0]]*(1-P.x-P.y)+RdHat::KHat[nvb[1]]*(P.x)+RdHat::KHat[nvb[2]]*(P.y) 
//	<< "," <<  RdHat::KHat[nvb[0]] << "," <<  RdHat::KHat[nvb[1]] << "," << RdHat::KHat[nvb[2]] <<endl;
  return RdHat::KHat[nvb[0]]*(1-P.x-P.y)+RdHat::KHat[nvb[1]]*(P.x)+RdHat::KHat[nvb[2]]*(P.y) ;}  

};

class Tet: public GenericElement<DataTet>  {
public: 
  Tet() {}; // constructor empty for array

  
  R3 H(int i) const 
  { ASSERTION(i>=0 && i <4);
    R3 AB(at(this->nvface[i][0]),at(this->nvface[i][1]));
    R3 AC(at(this->nvface[i][0]),at(this->nvface[i][2]));
    return AB^AC/(6.*this->mesure());} // heigth 
 
    R3 n(int i) const 
    { ASSERTION(i>=0 && i <4);
	R3 AB(at(this->nvface[i][0]),at(this->nvface[i][1]));
	R3 AC(at(this->nvface[i][0]),at(this->nvface[i][2]));
	R3 N=AB^AC;
    return N/N.norme();} //  exterior normal  
    
  void Gradlambda(R3 * GradL) const
  {
    R3 V1(at(0),at(1));
    R3 V2(at(0),at(2));
    R3 V3(at(0),at(3));
    R det1=1./(6.*mesure());
    GradL[1]= (V2^V3)*det1;
    GradL[2]= (V3^V1)*det1;
    GradL[3]= (V1^V2)*det1;
    GradL[0]=-GradL[1]-GradL[2]-GradL[3];
  }

};

class Triangle3: public GenericElement<DataTriangle3>  {
public: 
  Triangle3() {}; // constructor empty for array

  Rd Edge(int i) const {ASSERTION(i>=0 && i <3);
    return Rd(this->at((i+1)%3),this->at((i+2)%3));}// opposite edge vertex i
  /*
  Rd H(int i) const { ASSERTION(i>=0 && i <3);
    Rd E=Edge(i);return E.perp()/(2.*this->mesure());} // heigth 
  
  void Gradlambda(Rd * GradL) const
  {
    GradL[1]= H(1);
    GradL[2]= H(2);
    GradL[0]=-GradL[1]-GradL[2];
  }
  */ 

};


class Mesh3 : public GenericMesh<Tet,Triangle3,Vertex3> { 
public:
  Mesh3(){}
  Mesh3(const string); 
  Mesh3(FILE *f);     
  Mesh3(int nnv, int nnt, int nnbe, Vertex3 *vv, Tet *tt, Triangle3 *bb); 
  Mesh3(int nnv, int nnbe, Vertex3 *vv, Triangle3 *bb);  // surface mesh 

  void GSave(FILE * f) const ;
  void GRead(FILE * f);
  int Save(const string & filename) const ;  
  int SaveSurface(const string & filename) const ;  
  int SaveSurface(const string & filename1, const string & filename2) const ;  
  void flipSurfaceMesh3(int surface_orientation);
  void read(istream &);
  //int Popen(const FILE *popenstream);

private:
  int load(const string & filename); 
  Mesh3(const Mesh3 &); // pas de construction par copie
  void operator=(const Mesh3 &);// pas affectation par copy 
};
// for the caracteristic method.
    int  WalkInTet(const Mesh3 & Th,int it, R3 & Phat,const R3 & U, R & dt);

} 
#endif