splitmesh6.cpp

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

// $Id: splitmesh6.cpp,v 1.3 2008/05/14 12:48:24 hecht Exp $

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

  using namespace  Fem2D;

Mesh * SplitMesh6(Fem2D::Mesh * const & pTh)
{
  assert(pTh);
  const Mesh & Th(*pTh);  // le maillage d'origne a decoupe
  using  Fem2D::Triangle;
  using  Fem2D::Vertex;
  using  Fem2D::R2;
  using  Fem2D::BoundaryEdge;
  using  Fem2D::Mesh;
 // using  Fem2D::R;
  int nbv=Th.nv; // nombre de sommet 
  int nbt=Th.nt; // nombre de triangles
  int neb=Th.neb; // nombre d'aretes fontiere
  // allocation des nouveaux items du maillage  
  int nbe = 0;
  for(int k=0; k <nbt; ++k)
    for (int e = 0; e<3;++e)
    {
	int ee=e;
	int kk = Th.ElementAdj(k,ee);
	if( kk<=k) nbe++;
    }
  Vertex * v= new Vertex[nbv+nbt+nbe];
  Triangle *t= new Triangle[nbt*6];
  BoundaryEdge *b= new BoundaryEdge[neb*2];
  // generation des nouveaus sommets 
  Vertex *vv=v;
  KN<int> mm(3*nbt);
  // copie des anciens sommets (remarque il n'y a pas operateur de copy des sommets)
  for (int i=0;i<nbv;i++)
   {
     Vertex & V=Th(i);
     vv->x=V.x;
     vv->y=V.y;
     vv->lab = V.lab;
     vv++;      
   }
  // generation des points barycentre de trianngles 
  for (int k=0;k<nbt;k++)
    {
      Triangle & K=Th[k];
      R2 G= ( (R2) K[0] + K[1] + K[2] )  / 3.;
      vv->x=G.x;
      vv->y=G.y;
      vv->lab = 0;
      vv++;
    }   
// generation des milieux des cote
  int nn=0;
  for(int k=0; k <nbt; ++k)
      for (int e = 0; e<3;++e)
      {
	  int ee=e;
	  int kk = Th.ElementAdj(k,ee);
	  if( (kk>=k) || (kk <0) ) { 
	      int v0=Th(k,EdgesVertexTriangle[e][0]);
	      int v1=Th(k,EdgesVertexTriangle[e][1]);
	      R2  M = ((R2) Th(v0) + Th(v1)) /2.;
	      int lab=0;
	      BoundaryEdge *be=Th.TheBoundaryEdge(v0,v1);
	      if( be) 
	      {
		  lab= be->lab;
	      }
		
	      vv->x=M.x;
	      vv->y=M.y;
	      vv->lab = lab;
	      
	       mm[k*3+e]= vv-v;// numero du sommet
	       vv++;      
	       nn++;
	  }
	  else   mm[k*3+e] = mm[kk*3+ee];
      }
	  cout << " nb edge = " << nbe << " == " << nn << endl;
  ffassert(nbe==nn);	
  
  //  generation des triangles 
  Triangle *tt= t; 
  int nberr=0;
   
  for (int i=0;i<nbt;i++)
    {
      int i0=Th(i,0), i1=Th(i,1),i2=Th(i,2);
      int j0=mm[i*3], j1=mm[i*3+1],j2=mm[i*3+2];
      int ii = nbv + i; // numero du 
      // les 3 triangles par triangles origines 
      (*tt++).set(v,ii,i1,j0,Th[i].lab);
      (*tt++).set(v,ii,j0,i2,Th[i].lab);
      (*tt++).set(v,i0,ii,j1,Th[i].lab);
      (*tt++).set(v,j1,ii,i2,Th[i].lab);
      (*tt++).set(v,i0,j2,ii,Th[i].lab);
      (*tt++).set(v,j2,i1,ii,Th[i].lab);
    }  

  // les arete frontieres qui n'ont pas change
  BoundaryEdge * bb=b;
  for (int i=0;i<neb;i++)
    {        
      int ki;
      int k=Th.BoundaryElement(i,ki);
      int i1=Th(Th.bedges[i][0]);
      int i2=Th(Th.bedges[i][1]);
      int ii = mm[3*k+ki]; 
      int lab=Th.bedges[i].lab;     
      *bb++ = BoundaryEdge(v,i1,ii,lab);   
      *bb++ = BoundaryEdge(v,ii,i2,lab);   
    }
  //  generation de la class Mesh a partir des 3 tableaux : v,t,b
  {
    Mesh * m = new Mesh(nbv+nbt+nbe,nbt*6,neb*2,v,t,b);
    R2 Pn,Px;
    m->BoundingBox(Pn,Px);
    m->quadtree=new Fem2D::FQuadTree(m,Pn,Px,m->nv);
    m->decrement();
    return m;
  }
}

//  truc pour que la fonction 
// Init::Init() soit appele a moment du chargement dynamique
// du fichier 
//  
class Init { public:
  Init();
};

static Init init;  //  une variable globale qui serat construite  au chargement dynamique 

Init::Init(){  // le constructeur qui ajoute la fonction "splitmesh3"  a freefem++ 
  if (verbosity)
    cout << " lood: Split6  " << endl;
  Global.Add("splitmesh6","(",new OneOperator1_<Mesh *,Mesh *>(SplitMesh6));
  // utilisation 
  // mesh Th,Th3;
  //  ... construction du maillage Th ici 
  //   Th3=splitmesh3(Th);
  /*  example complet : splitmesh3.edp
      load "splitmesh3"
      mesh Th=square(5,5);
      mesh Th3=splitmesh3(Th);
      plot(Th3,wait=1);
  */

}