bamgfreefem.cpp

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

// ORIG-DATE:     Dec 97
// -*- Mode : c++ -*-
//
// SUMMARY  :      
// USAGE    :        
// ORG      : 
// 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
 */

// #define TRACETRIANGLE 3
//#pragma dont_inline on
//#pragma global_optimizer off
//#pragma inline_depth(0)

#undef NDEBUG
extern long verbosity ;
//#define strcasecmp strcmp
#include <cstdio>
#include <string>
#include <cmath>
#include <ctime>
#include <iomanip>
#include <fstream>

using namespace std;

#include "Meshio.h"
#include "Mesh2.h"
#include "QuadTree.h"
#include "SetOfE4.h"

#include "rgraph.hpp"
#include "fem.hpp"
#include "AFunction.hpp"
#include "BamgFreeFem.hpp"
#include "RNM.hpp"
#include "FESpace.hpp"
#include "Mesh3dn.hpp"
#include "MeshPoint.hpp"
#include "PlotStream.hpp"
#include <set>
Fem2D::Mesh *bamg2msh( bamg::Triangles* tTh,bool renumbering)
{ 
  using namespace bamg;
  bamg::Triangles & th (*tTh);
  tTh->ReNumberingTheTriangleBySubDomain(!renumbering);//  just compress 
  //tTh->NbRef++;
  Int4  i,j,k=0;
  int nv  =  tTh->nbv;
  int nt  =   tTh->nbt - tTh->NbOutT;
  int neb =   tTh->nbe;
  
  int nbcrakev = 0;
  tTh->ReMakeTriangleContainingTheVertex();
  Fem2D::Triangle * t =  new Fem2D::Triangle[nt]  ;
  Fem2D::BoundaryEdge * b_e = new Fem2D::BoundaryEdge[neb];
  
  Fem2D::Vertex vbase;        
  Fem2D::Vertex *vb(&vbase);
  if (verbosity>5)
    cout << "  -- Before cracking mesh:  Nb Triangles = " << nt << " Nb of Vertices " << nv << endl;
  for ( i=0;i<nt;i++) // unset all triangles 
    for (j=0;j<3;j++)
      t[i](j)=0;
  //  nv=0;    
  for (int iv=0;iv<th.nbv;iv++) // vertex 
    {
      // cout << iv << " : " ;
      const Vertex & v(th[iv]); 
      int kk=0; // nb cracked
      int kc=0; 
      int kkk =0; // nb triangle  with same number 
      Triangle * tbegin = v.t;
      Fem2D::Vertex * vv = vb+iv;
      int i  = v.vint;       
      throwassert(tbegin && (i >= 0 ) && (i <3));
      // turn around the vertex v
      TriangleAdjacent ta(tbegin,EdgesVertexTriangle[i][0]);// previous edge
      int k=0;
      do {
        int kv = VerticesOfTriangularEdge[ta][1];
        k++; 
        Triangle * tt (ta);
        throwassert( &v == & (*  tt)[kv] );            
        if ( ta.Cracked() ) 
          {   // cout << " || "    ;                    
            if ( kk == 0) tbegin=ta,kkk=0;  //  begin by a cracked edge  => restart                
            if (  kkk ) { kc =1;vv = vb +  nv++;  kkk = 0; } // new vertex if use 
            kk++;
            // number of cracked edge view                 
          }
        if ( tt->link ) { // if good triangles store the value 
          int it = th.Number(tt);
          throwassert(it < nt);
          //int iiv=vv-vb;
          t[it](kv) = vv;
          /*
          cout << it << " " << kv << " "<< iiv  << endl;
          if (&th(it)[kv] != &th[iiv])
             cout << it << " " << kv << " "<< iiv << " != " << th.Number(th(it)[kv]) << endl ;
          */
          kkk++;
        } else if (kk) { // crack + boundary 
          if (  kkk ) { kc =1;vv = vb +  nv++;  kkk = 0; } // new vertex if use 
        }
        
        ta = Next(ta).Adj(); 
      } while ( (tbegin != ta)); 
      throwassert(k);
      if (kc)  nbcrakev++;
    }
  Fem2D::Vertex * v = new Fem2D::Vertex[nv];
  //  set the vertices --
  for (i=0;i<nt;i++)
    { 
      for (j=0;j<3;j++)
        {
          throwassert( t[i](j) );             
          int k = t[i](j) - vb;
          t[i](j) = v+ k;
          throwassert(k>=0 && k < nv);
          Vertex & thv(th(i)[j]);
          v[k].x    =  thv.r.x;
          v[k].y    =  thv.r.y;
          v[k].lab  =  thv.ref();        
        }  
    }
  // warning in cracked edges 
  // construction of the edges --
  
  if (nbcrakev && verbosity>2)
    cout << " -- Nb of craked vertices = " << nbcrakev << " Nb of created vertices " << nv - th.nbv << endl;
  
  
  for (i=0;i<tTh->nbe;i++)
    {
      int i0=tTh->Number(tTh->edges[i][0]),i1=tTh->Number(tTh->edges[i][1]);
      throwassert(i0>=0 && i0 <nv);
      throwassert(i1>=0 && i1 <nv);      
      b_e[i]=Fem2D::BoundaryEdge(v,i0,i1,tTh->edges[i].ref);
    }      
  Int4 *reft = new Int4[tTh->nbt];
  //Int4 nbref =
  tTh->ConsRefTriangle(reft);
  for( i=0,k=0;i<tTh->nbt;i++)
    if(tTh->triangles[i].link)
      { 
        
        Fem2D::R2 A(t[k][0]),B(t[k][1]),C(t[k][2]);
        t[k].area = (( B-A)^(C-A))*0.5 ;
        t[k].lab = tTh->subdomains[reft[i]].ref;  // a faire
        throwassert(k == i);
        k++;
      }
  delete [] reft;
  throwassert ( nt == k);
  tTh->ReMakeTriangleContainingTheVertex();
  
  if (verbosity)
    cout << "  --  mesh:  Nb of Triangles = "  << setw(6) <<  nt << ", Nb of Vertices " << nv << endl;
  
  {  
    Fem2D::Mesh *m = new Fem2D::Mesh(nv,nt,neb,v,t,b_e);
    if (renumbering) m->renum();
    m->MakeQuadTree();
    return m;
  }
}

Fem2D::Mesh *bamg2msh(const bamg::Geometry &Gh)
{ 
  // ------------------
  int nv=  Gh.nbv;
  int neb=Gh.nbe;
  Fem2D::Triangle * t = 0 ;
  Fem2D::BoundaryEdge * b_e = new Fem2D::BoundaryEdge[neb];
  Fem2D::Vertex *v = new Fem2D::Vertex[nv]  ;
  for (int i=0;i<nv;i++)
    {
      const bamg::GeometricalVertex & vg( Gh[i]);
      v[i].x=vg.r.x;
      v[i].y=vg.r.y;
      v[i].lab=vg.ref();
      
    }     
  for (int ie=0;ie<neb;ie++)
    {
      const bamg::GeometricalEdge & e= Gh(ie);
      int i0=Gh.Number(e[0]),i1=Gh.Number(e[0]);
      b_e[ie]= Fem2D::BoundaryEdge(v,i0,i1,e.ref);
    }   
  
  {  
    Fem2D::Mesh *m = new Fem2D::Mesh(nv,0,neb,v,t,b_e);
    m->MakeQuadTree();
    return m;
  }
  // ------------------
}





bamg::Triangles * msh2bamg(const Fem2D::Mesh & Th,double cutoffradian,long * reqedgeslab,int nreqedgeslab)
  
{
  using namespace bamg;
  Triangles *Tn=new Triangles(Th.nv);
  Tn->nbv = Th.nv;
  Tn->nbt = Th.nt;
  Tn->nbe = Th.neb;
  Tn->name= new char[strlen("msh2bamg")+1];
  strcpy(Tn->name,"msh2bamg");
  //  Tn->triangles = new Triangle [Tn->nbtx];
  throwassert(Tn->triangles);
  //  Tn->vertices = new Vertex [Tn->nbvx];
  //  Tn->ordre = new (Vertex* [Tn->nbvx]);
  Tn->edges = new Edge [Th.neb];
 
  Int4 i;
  Metric Mid(1.);
  for (i = 0; i < Th.nv; i++)
    {
      Tn->vertices[i].r.x = Th(i).x;
      Tn->vertices[i].r.y = Th(i).y;
	Tn->vertices[i].m=Mid;
      Tn->vertices[i].ReferenceNumber = Th(i).lab;
    }
  
  //  Int4 i1 [nbt],i2 [nbt],i3 [nbt];
  for (i = 0; i < Th.nt; i++)
    {
      int i1 = Th(Th[i][0]);
      int i2 = Th(Th[i][1]);
      int i3 = Th(Th[i][2]);
      Tn->triangles[i]= Triangle( Tn,i1 ,i2 ,i3 );
      Tn->triangles[i].color = Th[i].lab;
    }
    //  Real8 cutoffradian = -1;    
    // add code   un change boundary part ...  frev 2009 JYU FH
    set<int> labreq;
    if(nreqedgeslab && verbosity) cout << " label of required edges " ;
    for (int i=0; i <nreqedgeslab;++i)
      {
	  if(verbosity)
	      cout << " " << reqedgeslab[i];
	  labreq.insert(reqedgeslab[i]);
      }
    bamg::GeometricalEdge paszero;  // add JYU    fevr 2009   for  required edge ....
    if(nreqedgeslab && verbosity) cout << endl;
    int k=0;  
    for (i = 0; i < Th.neb; i++)
    {
      Tn->edges[i].v[0] = Tn->vertices + Th(Th.bedges[i][0]);
      Tn->edges[i].v[1] = Tn->vertices + Th(Th.bedges[i][1]);
      Tn->edges[i].ref = Th.bedges[i].lab;
      Tn->edges[i].on = 0; 
      if( labreq.find( Tn->edges[i].ref) != labreq.end())
	{
	    k++; 
	    Tn->edges[i].on = &paszero; 
	}

    }
  if(verbosity)cout << "  number of required edges : "<< k << endl;
   
    
  Tn->ConsGeometry(cutoffradian);   
  Tn->Gh.AfterRead();    
  Tn->SetIntCoor();
  Tn->FillHoleInMesh();
  return Tn;
}


bamg::Triangles * msh2bamg(const Fem2D::Mesh & Th,double cutoffradian, 
                           int  nbdfv, int * ndfv,int  nbdfe, int * ndfe,
			   long * reqedgeslab,int nreqedgeslab)
{
  using namespace bamg;
  Triangles *Tn=new Triangles(Th.nv);
  KN<int> equiedges(Th.neb);
  for(int i=0;i<Th.neb;i++)
    equiedges[i]=2*i;
  if(nbdfe !=0 )
  {
  KN<int>  kk(Th.neb),kn(Th.neb);
   kk=0;
    for(int i=0;i<Th.neb;i++)
      {
         int df=ndfe[i];
         kk[df]++;
         if(kk[df]==1) kn[df]=i;
         else { 
           int k=kn[df],sens=0;
           int di0=ndfv[Th(Th.bedges[i][0])];
           int di1=ndfv[Th(Th.bedges[i][1])];
           int dk0=ndfv[Th(Th.bedges[k][0])];
           int dk1=ndfv[Th(Th.bedges[k][1])];
           if ((di0==dk0) &&(di1==dk1) ) sens=0;
           else if ((di1==dk0) &&(di0==dk1) ) sens=1;
           else  {
              cout << "Error in periodic mesh " << di0 << " " << di1 << " <=> " << dk0 << " " << dk1 << endl;
              ExecError("bug periodic mesh in ??? ");
           }
           equiedges[i]=2*k+sens;

         }
      }
    
  }; // a faire pour les maillages periodique 
  
  Tn->nbv = Th.nv;
  Tn->nbt = Th.nt;
  Tn->nbe = Th.neb;
  Tn->name= new char[strlen("msh2bamg")+1];
  strcpy(Tn->name,"msh2bamg");
  //  Tn->triangles = new Triangle [Tn->nbtx];
  throwassert(Tn->triangles);
  //  Tn->vertices = new Vertex [Tn->nbvx];
  //  Tn->ordre = new (Vertex* [Tn->nbvx]);
  Tn->edges = new Edge [Th.neb];
  
  Int4 i;
    Metric Mid(1.);  
  for (i = 0; i < Th.nv; i++)
    {
      Tn->vertices[i].r.x = Th(i).x;
      Tn->vertices[i].r.y = Th(i).y;
      Tn->vertices[i].ReferenceNumber = Th(i).lab;
      Tn->vertices[i].m=Mid;
    }
  
  //  Int4 i1 [nbt],i2 [nbt],i3 [nbt];
  for (i = 0; i < Th.nt; i++)
    {
      int i1 = Th(Th[i][0]);
      int i2 = Th(Th[i][1]);
      int i3 = Th(Th[i][2]);
      Tn->triangles[i]= Triangle( Tn,i1 ,i2 ,i3 );
      Tn->triangles[i].color = Th[i].lab;
    }
    
    // add code   un change boundary part ...  frev 2009 JYU FH
    set<int> labreq;
    if(nreqedgeslab && verbosity) cout << " label of required edges " ;
    for (int i=0; i <nreqedgeslab;++i)
      {
	  if(verbosity)
	      cout << " " << reqedgeslab[i];
	  labreq.insert(reqedgeslab[i]);
      }
    bamg::GeometricalEdge paszero;  // add JYU    fevr 2009   for  required edge ....
    if(nreqedgeslab && verbosity) cout << endl;
    int k=0;  
    
  for (i = 0; i < Th.neb; i++)
    {
      Tn->edges[i].v[0] = Tn->vertices + Th(Th.bedges[i][0]);
      Tn->edges[i].v[1] = Tn->vertices + Th(Th.bedges[i][1]);
      Tn->edges[i].ref = Th.bedges[i].lab;
      Tn->edges[i].on = 0; 
      if( labreq.find( Tn->edges[i].ref) != labreq.end())
	  {
	      k++; 
	      Tn->edges[i].on = &paszero; 
	  }
    }
  //  Real8 cutoffradian = -1;
  Tn->ConsGeometry(cutoffradian,equiedges);
  Tn->Gh.AfterRead();    
  Tn->SetIntCoor();
  Tn->FillHoleInMesh();
  return Tn;
}



Fem2D::Mesh *  BuildMesh(Stack stack, E_BorderN const * const & b,bool justboundary,int nbvmax,bool Requiredboundary) 
{
  

  using namespace bamg;
  using bamg::Abs;
  using bamg::Max;
  using bamg::Min;
  using bamg::Pi;
  Fem2D::MeshPoint & mp (*Fem2D::MeshPointStack(stack)), mps = mp;
  
  int nbvx=0,nbe=0,nbsd=0;
  for (E_BorderN const * k=b;k;k=k->next)
    {long n=  Max(1L,Abs(k->Nbseg(stack))); ;
    nbvx += n+1;
    nbe += n;
    nbsd++;
    }
  
  Geometry * Gh =  new Geometry;
  
  if(verbosity>2)
    cout <<"\t\t"  << "  Begin: ConstGeometry from nb Border  "  << nbsd <<endl;
  //throwassert(empty());
  const char * filename = "FREEFEM.gh";
  Gh->name=new char [strlen(filename)+1];
  strcpy(Gh->name,filename);
  Real8 Hmin = HUGE_VAL;// the infinie value 
  Int4 hvertices =0;
  Int4 i,nn,n;
  //Int4 dim=0;
  Gh->MaximalAngleOfCorner =30.00*Pi/180.0;
  Gh->nbv = 0;
  Gh->nbvx = nbvx;
  
  Gh->nbe = nbe;
  Gh->edges = new GeometricalEdge[Gh->nbe];
  bamg::Vertex *vertices =  new GeometricalVertex[Gh->nbvx];
  double lmin= HUGE_VAL;
  //  generation des points et des lignes 
  i=0;
  for (E_BorderN const * k=b;k;k=k->next)
    {
      assert(k->b->xfrom); // a faire 
      double & t = *  k->var(stack);
      double a(k->from(stack)),b(k->to(stack));
      n=Max(Abs(k->Nbseg(stack)),1L);