lgmesh.cpp

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

// -*- 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
 */
/*
#include  <iostream>
#include  <cfloat>
using namespace std;
#include "error.hpp"
#include "AFunction.hpp"
#include "rgraph.hpp"

#include "RNM.hpp"
#include "fem.hpp"


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


//#include "fem3.hpp"
#include "MeshPoint.hpp"
#include <complex>
#include "Operator.hpp" 

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

#include "lex.hpp"

#include "Mesh2.h"

#include "BamgFreeFem.hpp"
#include "lgfem.hpp"
*/
#include "ff++.hpp"

using Fem2D::Mesh;
using Fem2D::MeshPoint;

extern bool NoWait; 

typedef Mesh * pmesh;




class classBuildMesh :  public E_F0mps { public:  
 
   typedef pmesh  Result;

   static basicAC_F0::name_and_type name_param[] ;
   static const int n_name_param =2;
   
    Expression nargs[n_name_param];
   
   Expression getborders;
   
   long arg(int i,Stack stack,long a) const{ return nargs[i] ? GetAny<long>( (*nargs[i])(stack) ): a;}
   bool arg(int i,Stack stack,bool a) const{ return nargs[i] ? GetAny<bool>( (*nargs[i])(stack) ): a;}
    
    classBuildMesh(const basicAC_F0 & args) 
    {   
      args.SetNameParam(n_name_param,name_param,nargs);
      getborders=to<const E_BorderN *>(args[0]); 
     }   
     
    static ArrayOfaType  typeargs() { return  ArrayOfaType(atype<const E_BorderN *>());}
    static  E_F0 * f(const basicAC_F0 & args){ return new classBuildMesh(args);} 
    AnyType operator()(Stack s) const ;
    operator aType () const { return atype<Result>();} 

};

basicAC_F0::name_and_type  classBuildMesh::name_param[]= {
    {  "nbvx", &typeid(long)} ,
    {"fixeborder", &typeid(bool)}  
};
// modif aout 2007
class BuildMeshFile :  public E_F0mps { public:  
    
    typedef pmesh  Result;
    
    static basicAC_F0::name_and_type name_param[] ;
    static const int n_name_param =1;
    
    Expression nargs[n_name_param];
    
    Expression getfilename;
    
    long arg(int i,Stack stack,long a) const{ return nargs[i] ? GetAny<long>( (*nargs[i])(stack) ): a;}
    
    BuildMeshFile(const basicAC_F0 & args) 
    {   
	args.SetNameParam(n_name_param,name_param,nargs);
	getfilename=to<string* >(args[0]); 
    }   
    
    static ArrayOfaType  typeargs() { return  ArrayOfaType(atype<string *>());}
    static  E_F0 * f(const basicAC_F0 & args){ return new BuildMeshFile(args);} 
    AnyType operator()(Stack s) const ;
    operator aType () const { return atype<Result>();} 
    
};

basicAC_F0::name_and_type  BuildMeshFile::name_param[]= {
    {  "nbvx", &typeid(long) }
};
// fin modif 2007 

class MoveMesh :  public E_F0mps { public:  
 
   typedef pmesh  Result;
   Expression getmesh;
   Expression U,V;
   int nbsol;    
    vector<Expression> sol;
   
    MoveMesh(const basicAC_F0 & args) :nbsol(args.size()-2),sol(args.size()-2)
    {   
  using  Fem2D::Triangle;
  using  Fem2D::Vertex;
  using  Fem2D::R2;
  using  Fem2D::BoundaryEdge;
  using  Fem2D::Mesh;
 // using  Fem2D::R;
  using  Fem2D::MeshPointStack;
      args.SetNameParam();
      getmesh=to<pmesh>(args[0]); 
      const E_Array * a = dynamic_cast<const E_Array *>(args[1].LeftValue());
      
      ffassert(a);
      if (a->size() !=2) CompileError("movemesh(Th,[u,v],...) need 2 componate in array ",atype<pmesh>());
      U=to<double>( (*a)[0]);
      V=to<double>( (*a)[1]);
      
      for (int i=2;i<args.size();i++)
        sol[i-2]=to<double>(args[i]);      
    }   
    static ArrayOfaType  typeargs() { return  ArrayOfaType(atype<pmesh>(),atype<E_Array>(),true);}
    static  E_F0 * f(const basicAC_F0 & args){ return new MoveMesh(args);} 
    AnyType operator()(Stack s) const ;
      operator aType () const { return atype<Result>();} 

};

class SplitMesh :  public E_F0mps { public:  
 
   typedef pmesh  Result;
   Expression getmesh;
   Expression U;
   
    SplitMesh(const basicAC_F0 & args) 
    {   
      args.SetNameParam();
      getmesh=to<pmesh>(args[0]); 
      U = to<long>(args[1]);      
     }   
    static ArrayOfaType  typeargs() { return  ArrayOfaType(atype<pmesh>(),atype<long>());}
    static  E_F0 * f(const basicAC_F0 & args){ return new SplitMesh(args);} 
    AnyType operator()(Stack s) const ;
    operator aType () const { return atype<Result>();} 

};

class SaveMesh :  public E_F0 { public:  
 
   typedef pmesh  Result;
   Expression getmesh;
   Expression filename; 
   Expression xx,yy,zz;  
   SaveMesh(const basicAC_F0 & args) 
    {   
      xx=0;
      yy=0;
      zz=0;
      args.SetNameParam();
      getmesh=to<pmesh>(args[0]); 
      filename=to<string*>(args[1]); 
      if (args.size() >2) 
        {
          const E_Array * a = dynamic_cast<const E_Array *>(args[2].LeftValue());
          if (!a) CompileError("savemesh(Th,\"filename\",[u,v,w],...");
          int k=a->size() ;
         // cout << k << endl;
          if ( k!=2 && k !=3) CompileError("savemesh(Th,\"filename\",[u,v,w]) need 2 or 3  componate in array ",atype<pmesh>());
          xx=to<double>( (*a)[0]);
          yy=to<double>( (*a)[1]);
          if(k==3)
           zz=to<double>( (*a)[2]);
         }
      
   }   
    static ArrayOfaType  typeargs() { return  ArrayOfaType(atype<pmesh>(),atype<string*>(),true);}
    static  E_F0 * f(const basicAC_F0 & args){ return new SaveMesh(args);} 
    AnyType operator()(Stack s) const ;
  
};


class Adaptation :   public E_F0mps { public:
  typedef pmesh  Result;
  
  static basicAC_F0::name_and_type name_param[] ;
  static const int n_name_param =27;
  
  int nbsol;    
  Expression nargs[n_name_param];
  Expression getmesh;
  Expression em11,em22,em12;
  int  typesol[100];
  vector<Expression> sol;
  int nbcperiodic;
  Expression *periodic;
  
  
  double arg(int i,Stack stack,double a) const { return nargs[i] ? GetAny<double>( (*nargs[i])(stack) ): a;}
  long arg(int i,Stack stack,long a) const{ return nargs[i] ? GetAny<long>( (*nargs[i])(stack) ): a;}
  bool arg(int i,Stack stack,bool a) const{ return nargs[i] ? GetAny<bool>( (*nargs[i])(stack) ): a;}
  int arg(int i,Stack stack,int a) const{ return nargs[i] ? GetAny<int>( (*nargs[i])(stack) ): a;}
  
  Adaptation(const basicAC_F0 & args) :nbsol(args.size()-1),sol(args.size()-1)
  {   
      em11=0;
      em22=0;
      em12=0;
      
      args.SetNameParam(n_name_param,name_param,nargs);
      getmesh=to<pmesh>(args[0]); 
      int ksol=0; 
      ffassert(nbsol<100);
      for (int i=1;i<nbsol+1;i++)       
         if (args[i].left()==atype<E_Array>())
          {
            const E_Array * a = dynamic_cast<const E_Array *>(args[i].LeftValue());
            ffassert(a);
            ksol+=a->size(); 
          }
         else
           ksol++;
      sol.resize(ksol); 
      ksol=0; 
      for (int i=1;i<nbsol+1;i++)       
         if (args[i].left()==atype<E_Array>())
          {
            const E_Array * a = dynamic_cast<const E_Array *>(args[i].LeftValue());
            ffassert(a);
             int N=a->size();
            typesol[i-1]=N-1; // ok en 2D
            if (N<=4) {
             for (int j=0;j<N;j++)             
              sol[ksol++]=to<double>((*a)[j]); }
            else {
             lgerror(" Adaptation vecteur a plus de 4 compossantes inconnue");
            }
              
          }
         else {
          typesol[i-1]=0;
          sol[ksol++]=to<double>(args[i]);
          }
    const E_Array * expmetrix = dynamic_cast<const E_Array *>(nargs[24]);
  
  if(expmetrix)
   {
      if(expmetrix->nbitem()!=3)
        ExecError("\nSorry we wait an array with 3 componants in: metrix=[m11,m12,m22]");
        
      em11=(*expmetrix)[0];
      em12=(*expmetrix)[1];
      em22=(*expmetrix)[2];
      
      if(  (*expmetrix)[0].left()!= atype<KN<double> *>() )
          CompileError("Sorry the fist  array componant in metrix=[m11,m12,m22] must be vector");
      if( (*expmetrix)[1].left()!= atype<KN<double> *>() )
          CompileError("Sorry the second  array componant in metrix=[m11,m12,m22] must be vector");
      if( (*expmetrix)[2].left()!= atype<KN<double> *>() ) 
          CompileError("Sorry the third  array componant in metrix=[m11,m12,m22] must be vector");
          
          
     }
     nbcperiodic=0;
     periodic=0;
     GetPeriodic(nargs[25],nbcperiodic,periodic);
       
   }  
    
    static ArrayOfaType  typeargs() { return  ArrayOfaType(atype<pmesh>(),true);}
    static  E_F0 * f(const basicAC_F0 & args){ return new Adaptation(args);} 
    AnyType operator()(Stack s) const ;    
    operator aType () const { return atype<pmesh>();}         
};


 basicAC_F0::name_and_type Adaptation::name_param[Adaptation::n_name_param] = {
       {   "hmin",             &typeid(double)},  //