lgfem.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
 */
#ifdef __MWERKS__
#pragma optimization_level 0
//#pragma inline_depth(0)
#endif

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


#include "MatriceCreuse_tpl.hpp"


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

#include <set>
#include <map>
#include <vector>

#include "lex.hpp"
#include "lgfem.hpp"
#include "lgmesh3.hpp"
#include "lgsolver.hpp"
#include "problem.hpp"
#include "CGNL.hpp"
#include "AddNewFE.h"
#include "array_resize.hpp"
#include "PlotStream.hpp"

// add for the gestion of the endianness of the file.
//PlotStream::fBytes PlotStream::zott; //0123;
//PlotStream::hBytes PlotStream::zottffss; //012345678;
// ---- FH
namespace bamg { class Triangles; }
namespace Fem2D { void DrawIsoT(const R2 Pt[3],const R ff[3],const RN_ & Viso);
   extern GTypeOfFE<Mesh3> &P1bLagrange3d;
}

#include "BamgFreeFem.hpp"

static bool TheWait=false;
bool  NoWait=false;

extern long verbosity;
extern FILE *ThePlotStream; //  Add for new plot. FH oct 2008
void init_lgmesh() ;

namespace FreeFempp {
template<class R>
 TypeVarForm<R> * TypeVarForm<R>::Global;
}

const int nTypeSolveMat=10;
int kTypeSolveMat;
TypeSolveMat *dTypeSolveMat[nTypeSolveMat];


basicAC_F0::name_and_type  OpCall_FormBilinear_np::name_param[]= {
  {   "init", &typeid(bool)},
  {   "solver", &typeid(TypeSolveMat*)},
  {   "eps", &typeid(double) } ,
  {   "precon",&typeid(Polymorphic*)}, 
  {   "dimKrylov",&typeid(long)},
  {   "bmat",&typeid(Matrice_Creuse<R>* )},
  {   "tgv",&typeid(double )},
  {   "factorize",&typeid(bool)},
  {   "strategy",&typeid(long )},
  {   "tolpivot", &typeid(double)},
  {   "tolpivotsym", &typeid(double) },
  {   "nbiter", &typeid(long)} // 12 

     
          
};


basicAC_F0::name_and_type  OpCall_FormLinear_np::name_param[]= {
  "tgv",&typeid(double )
};


/*
template<class R,class TA0,class TA1,class TA2>
 class E_F_F0F0F0 :public  E_F0 { public:
   template <class T> struct remove_reference     {typedef T type;};
   template <class T> struct remove_reference<T&> {typedef T type;};
   typedef typename remove_reference<TA0>::type A0;
   typedef typename remove_reference<TA1>::type A1;
   typedef typename remove_reference<TA2>::type A2;
   typedef  R (*func)( A0 , A1, A2 ) ;
    
  func f;
  Expression a0,a1,a2;
  E_F_F0F0F0(func ff,Expression aa0,Expression aa1,Expression aa2) 
    : f(ff),a0(aa0),a1(aa1),a2(aa2) {}
  AnyType operator()(Stack s)  const 
    {return SetAny<R>( f( GetAny<A0>((*a0)(s)) , GetAny<A1>((*a1)(s)),  GetAny<A2>((*a2)(s)) ) );}  
   bool EvaluableWithOutStack() const 
      {return a0->EvaluableWithOutStack() && a1->EvaluableWithOutStack() && a2->EvaluableWithOutStack() ;} // 
   bool MeshIndependent() const 
      {return a0->MeshIndependent() && a1->MeshIndependent() && a2->MeshIndependent();} // 
  int compare (const E_F0 *t) const { 
     int rr;
    // cout << "cmp " << typeid(*this).name() << " and " << typeid(t).name() << endl;
     const  E_F_F0F0F0* tt=dynamic_cast<const E_F_F0F0F0 *>(t);
     if (tt && f == tt->f) rr= clexico(a0->compare(tt->a0),a1->compare(tt->a1),a2->compare(tt->a2));
     else rr = E_F0::compare(t);
     return rr;
     } // to give a order in instuction 

      
   int Optimize(deque<pair<Expression,int> > &l,MapOfE_F0 & m, size_t & n) const
    {

       int rr = find(m);
       if (rr) return rr;

       return insert(new Opt(*this,a0->Optimize(l,m,n),a1->Optimize(l,m,n),a2->Optimize(l,m,n)),l,m,n);
    }
     // build optimisation

     class Opt: public E_F_F0F0F0  { public :
       size_t ia,ib,ic;  
       Opt(const  E_F_F0F0F0 &t,size_t iaa,size_t ibb,size_t icc) 
         : E_F_F0F0F0(t) ,
         ia(iaa),ib(ibb),ic(icc) {}
      AnyType operator()(Stack s)  const 
       {
         //A0 aa =*static_cast<A0 *>(static_cast<void*>(static_cast<char *>(s)+ia));
         //A1 bb=*static_cast<A1 *>(static_cast<void*>(static_cast<char *>(s)+ib)) ;
         //cout << ia << " " << ib <<  "f( " << aa << "," << bb  << " )   = "<< f(aa,bb) << endl;
         return SetAny<R>( f( *static_cast<A0 *>(static_cast<void*>(static_cast<char *>(s)+ia)) , 
                             *static_cast<A1 *>(static_cast<void*>(static_cast<char *>(s)+ib)) ,
                             *static_cast<A2 *>(static_cast<void*>(static_cast<char *>(s)+ic)) 
                             ) );}  
         
      };     
       
    
};
template<class R,class A=R,class B=A,class C=A,class CODE=E_F_F0F0F0<R,A,B,C> >
class  OneOperator3 : public OneOperator {
    aType r; //  return type 
    typedef typename CODE::func func;
    //typedef  R (*func)(A,B) ;
    func f;
    public: 
    E_F0 * code(const basicAC_F0 & args) const 
     { return  new CODE(f,t[0]->CastTo(args[0]),t[1]->CastTo(args[1]),t[2]->CastTo(args[3]));} 
    OneOperator3(func  ff): 
      OneOperator(map_type[typeid(R).name()],map_type[typeid(A).name()],map_type[typeid(B).name()],map_type[typeid(C).name()]),
      f(ff){}
};
*/
const E_Array * Array(const C_F0 & a) { 
  if (a.left() == atype<E_Array>() )
    return dynamic_cast<const E_Array *>(a.LeftValue());
  else 
    return 0;
}
bool Box2(const C_F0 & bb, Expression * box)
{
  const E_Array * a= Array(bb);
  if(a && a->size() == 2)
   {
    box[0] = to<double>((*a)[0]);
    box[1] = to<double>((*a)[1]);
    return true;
    }
  else 
    return false;
  
}
bool Box2x2(Expression  bb, Expression * box)
{
  const E_Array * a= dynamic_cast<const E_Array *>(bb);
  if(a && a->size() == 2)
     return Box2((*a)[0],box)  &&  Box2((*a)[1],box+2) ;
  else 
    return false;
}
void dump_table()
{
   cout << " dump the table of the language " << endl;
   cout << " ------------------------------ " <<endl <<endl;
	map<const string,basicForEachType *>::const_iterator i; ;
	       
	for (i= map_type.begin();i !=map_type.end();i++)
	  {
	   cout << " type : " << i->first << endl;
	   if( i->second )
	     i->second->ShowTable(cout);
	   else cout << " Null \n";
	   cout << "\n\n";     
	  }

	for (i= map_type.begin();i !=map_type.end();i++)
	  {
	   cout << " type : " << i->first << endl;
	   if( i->second )
	     i->second->ShowTable(cout);
	   else cout << " Null \n";
	   cout << "\n\n";    
	   
	  }
	  
	 cout << "--------------------- " << endl;
	 cout << *TheOperators << endl;
	 cout << "--------------------- " << endl;

}
/*
 class LocalArrayVariable:public E_F0
 { 
  size_t offset;
  aType t; //  type of the variable just for check  
  Expression  n; // expression of the size 
  public:
  AnyType operator()(Stack s) const { 
    SHOWVERB( cout << "\n\tget var " << offset << " " <<  t->name() << endl);  
   return PtrtoAny(static_cast<void *>(static_cast<char *>(s)+offset),t);}
  LocalArrayVariable(size_t o,aType tt,Expression nn):offset(o),t(tt),n(nn) {ffassert(tt);
    SHOWVERB(cout << "\n--------new var " << offset << " " <<  t->name() << endl);
    }
};
*/






bool In(long *viso,int n,long v) 
{
  int i=0,j=n,k;
  if  (v <viso[0] || v >viso[j-1]) 
    return false;  
  while (i<j-1)    
   if ( viso[k=(i+j)/2]> v) j=k;
   else i=k;
  return (viso[i]=v);
}



class  LinkToInterpreter { public:
 Type_Expr   P,N,x,y,z,label,region,nu_triangle,nu_edge,lenEdge,hTriangle,area,inside;
  LinkToInterpreter() ; 
};

LinkToInterpreter * l2interpreter;


  using namespace Fem2D;
  using namespace EF23;
/*
inline pmesh  ReadMesh( string * const & s) {
  Mesh * m=new Mesh(*s);
  R2 Pn,Px;
  m->BoundingBox(Pn,Px);
  m->quadtree=new FQuadTree(m,Pn,Px,m->nv);
 //  delete s; modif FH 2006 (stack ptr)
  return m;
 }

inline pmesh3  ReadMesh3( string * const & s) {
  Mesh3 * m=new Mesh3(s->c_str());
  R3 Pn,Px;
  // m->BoundingBox(Pn,Px);
  m->gtree=new Mesh3::GTree(m->vertices,m->Pmin,m->Pmax,m->nv);
 //  delete s; modif FH 2006 (stack ptr)
  return m;
 }
*/
 

template<class Result,class A>
class E_F_A_Ptr_o_R :public  E_F0 { public:
  typedef Result A::* ptr;
  Expression a0;
  ptr p; 
  E_F_A_Ptr_o_R(Expression aa0,ptr pp) 
    : a0(aa0),p(pp) {}
  AnyType operator()(Stack s)  const {
    return SetAny<Result*>(&(GetAny<A*>((*a0)(s))->*p));}
  bool MeshIndependent() const {return a0->MeshIndependent();} // 
  
};
//  ----
//  remarque pas de template, cela ne marche pas encore ......
 class E_P_Stack_P   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const { throwassert(* (long *) s);
    return SetAny<R3*>(&MeshPointStack(s)->P);} 
    operator aType () const { return atype<R3*>();}         
    
}; 
 class E_P_Stack_Px   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const { throwassert(* (long *) s);
    return SetAny<R*>(&MeshPointStack(s)->P.x);} 
    operator aType () const { return atype<R*>();}         
    
}; 
 class E_P_Stack_Py   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const {throwassert(* (long *) s);
    return SetAny<R*>(&MeshPointStack(s)->P.y);} 
    operator aType () const { return atype<R*>();}         
    
}; 
 class E_P_Stack_Pz   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const { throwassert(* (long *) s);
    return SetAny<R*>(&MeshPointStack(s)->P.z);} 
    operator aType () const { return atype<R*>();}         
    
}; 

 class E_P_Stack_N   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const { throwassert(* (long *) s);
    return SetAny<R3*>(&MeshPointStack(s)->N);} 
        operator aType () const { return atype<R3*>();}         

}; 
 class E_P_Stack_Nx   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const { throwassert(* (long *) s);
    return SetAny<R*>(&MeshPointStack(s)->N.x);} 
    operator aType () const { return atype<R*>();}         
    
}; 
 class E_P_Stack_Ny   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const { throwassert(* (long *) s);
    return SetAny<R*>(&MeshPointStack(s)->N.y);} 
    operator aType () const { return atype<R*>();}             
}; 
 class E_P_Stack_Nz   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const { throwassert(* (long *) s);
    return SetAny<R*>(&MeshPointStack(s)->N.z);} 
    operator aType () const { return atype<R*>();}         
    
}; 

 class E_P_Stack_Region   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const { throwassert(* (long *) s);
    return SetAny<long*>(&MeshPointStack(s)->region);} 
    operator aType () const { return atype<long*>();}         
    
}; 
 class E_P_Stack_Label   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const { throwassert(* (long *) s);
    return SetAny<long*>(&MeshPointStack(s)->label);} 
    operator aType () const { return atype<long *>();}         
    
}; 
 class E_P_Stack_Mesh   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const { throwassert(* (long *) s);
    return SetAny<pmesh >(const_cast<pmesh>(MeshPointStack(s)->Th));} 
  operator aType () const { return atype<pmesh>();}         

}; 
 class E_P_Stack_Nu_Triangle   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const { throwassert(* (long *) s);
    return SetAny<long>(MeshPointStack(s)->t);} 
    operator aType () const { return atype<long>();}         
    
}; 
 class E_P_Stack_Nu_Vertex   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const { throwassert(* (long *) s);
    return SetAny<long>(MeshPointStack(s)->v);} 
    operator aType () const { return atype<long>();}         
    
}; 
 class E_P_Stack_Nu_Face   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const { throwassert(* (long *) s);
    return SetAny<long>(MeshPointStack(s)->f);} 
    operator aType () const { return atype<long>();}         
    
}; 
 class E_P_Stack_Nu_Edge   :public  E_F0mps { public: 
  AnyType operator()(Stack s)  const { throwassert(* (long *) s);
    return SetAny<long>(MeshPointStack(s)->e);} 
    operator aType () const { return atype<long>();}