lgmat.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 "error.hpp"
#include "RNM.hpp"
#include <set>
#include <vector>
#include <cstdio>
#include <cstring>
#include <complex>

#include  <cmath>
#include  <iostream>
using namespace std;
#include "FESpacen.hpp"

#include "AFunction.hpp"
#include "rgraph.hpp"
#include "MatriceCreuse_tpl.hpp"

//#include "fem3.hpp"
#include "MeshPoint.hpp"
#include "Operator.hpp" 

#include "lex.hpp"
#include "lgfem.hpp"
#include "lgmesh3.hpp"
#include "lgsolver.hpp"
#include "problem.hpp"
*/
#include "ff++.hpp" 
#include "array_resize.hpp" 

#include "CGNL.hpp"

namespace bamg { class Triangles; }
namespace Fem2D { void DrawIsoT(const R2 Pt[3],const R ff[3],const RN_ & Viso); }

extern const basicForEachType *aatypeknlongp; //// for  compilation error with g++ 3.2.2
#include "BamgFreeFem.hpp"



// Debut FH Houston -------- avril 2004 ---------
//  class for operator on sparse matrix 
//   ------------------------------------
//  pour le addition de matrice ----matrice creuse in french)
//  MatriceCreuse    real class for matrix sparse
//  Matrice_Creuse   class for matrix sparse +  poiteur on FE space def 
//         to recompute matrice in case of mesh change
//  list<triplet<R,MatriceCreuse<R> *,bool> > * is liste of 
//  \sum_i a_i*A_i  where a_i is a scalare and A_i is a sparse matrix
//
 


template<class R> 
list<triplet<R,MatriceCreuse<R> *, bool > > * to(Matrice_Creuse<R> * M)
{
  list<triplet<R,MatriceCreuse<R> *,bool> >  * l=new list<triplet<R,MatriceCreuse<R> *,bool> >;
   l ->push_back(make_triplet<R,MatriceCreuse<R> *,bool>(1,M->A,false));
  return  l;
}
template<class R> 
list<triplet<R,MatriceCreuse<R> *, bool > > * to(Matrice_Creuse_Transpose<R>  M)
{
  list<triplet<R,MatriceCreuse<R> *,bool> >  * l=new list<triplet<R,MatriceCreuse<R> *,bool> >;
   l ->push_back(make_triplet<R,MatriceCreuse<R> *,bool>(1,M.A->A,true));
  return  l;
}

template<class R>
AnyType M2L3 (Stack , const AnyType & pp)
{
    return to(GetAny<Matrice_Creuse<R> *>(pp));
}


template<class R>
AnyType tM2L3 (Stack , const AnyType & pp)
{
    return to(GetAny<Matrice_Creuse_Transpose<R> >(pp));
}


template<class R> 
struct Op2_ListCM: public binary_function<R,Matrice_Creuse<R> *,list<triplet<R,MatriceCreuse<R> *,bool> > *> 
 { 
   typedef triplet<R,MatriceCreuse<R> *,bool>  P;
   typedef list<P> L;
   typedef L * RR;
   typedef R  AA;
   typedef Matrice_Creuse<R> * BB;
   
 static  RR f(const   AA & a,const   BB & b)  
  {
    RR r=  new list<P> ;
    r ->push_back(make_triplet<R,MatriceCreuse<R> *>(a,b->A,false));
    return r;}
};

template<class R> 
struct Op2_ListMC: public binary_function<Matrice_Creuse<R> *,R,list<triplet<R,MatriceCreuse<R> *,bool> > *> 
 { 
   typedef triplet<R,MatriceCreuse<R> *,bool>  P;
   typedef list<P> L;
   typedef L * RR;
   typedef R  AA;
   typedef Matrice_Creuse<R> * BB;
   
 static  RR f(const   BB & b,const   AA & a)  
  {
    RR r=  new list<P> ;
    r ->push_back(make_triplet<R,MatriceCreuse<R> *>(a,b->A,false));
    return r;}
};
//  ADD FH 16/02/2007

template<class R> 
struct Op2_ListCMt: public binary_function<R,Matrice_Creuse_Transpose<R> ,list<triplet<R,MatriceCreuse<R> *,bool> > *> 
{ 
    typedef triplet<R,MatriceCreuse<R> *,bool>  P;
    typedef list<P> L;
    typedef L * RR;
    typedef R  AA;
    typedef Matrice_Creuse_Transpose<R>  BB;
    
    static  RR f(const   AA & a,const   BB & b)  
    {
	RR r=  new list<P> ;
	r ->push_back(make_triplet<R,MatriceCreuse<R> *>(a,b.A->A,true));
	return r;}
};

template<class R> 
struct Op2_ListMtC: public binary_function<Matrice_Creuse_Transpose<R> ,R,list<triplet<R,MatriceCreuse<R> *,bool> > *> 
{ 
    typedef triplet<R,MatriceCreuse<R> *,bool>  P;
    typedef list<P> L;
    typedef L * RR;
    typedef R  AA;
    typedef Matrice_Creuse_Transpose<R> BB;
    
    static  RR f(const   BB & b,const   AA & a)  
    {
	RR r=  new list<P> ;
	r ->push_back(make_triplet<R,MatriceCreuse<R> *>(a,b.A->A,true));
	return r;}
};
// FIN ADD 16/02/2007



template<class R> 
struct Op1_LCMd: public unary_function<list<triplet<R,MatriceCreuse<R> *,bool> > *,
list<triplet<R,MatriceCreuse<R> *,bool> > *  >
{  //  - ...
    typedef triplet<R,MatriceCreuse<R> *,bool>  P;
    typedef list<P> L;
    typedef L * RR;
    
    static   RR f(const RR & l)  
    { 
	typedef typename list<triplet<R,MatriceCreuse<R> *,bool> >::iterator lci;
	for (lci i= l->begin();i !=l->end();++i)
	    i->first *= R(-1);
	return l;
    }    
    
};

template<class R> 
struct Op2_ListCMCMadd: public binary_function<list<triplet<R,MatriceCreuse<R> *,bool> > *,
                                               list<triplet<R,MatriceCreuse<R> *,bool> > *,
                                               list<triplet<R,MatriceCreuse<R> *,bool> > *  >
{  //  ... + ...
   typedef triplet<R,MatriceCreuse<R> *,bool>  P;
   typedef list<P> L;
   typedef L * RR;

  static   RR f(const RR & a,const RR & b)  
  { 
    a->insert(a->end(),b->begin(),b->end());
      
    delete b;
    return a;
  }    
    
};

template<class R> 
struct Op2_ListMCMadd: public binary_function<Matrice_Creuse<R> *,
                                              list<triplet<R,MatriceCreuse<R> *,bool> > *,                                               
                                               list<triplet<R,MatriceCreuse<R> *,bool> > *  >
{  //  M + ....
   typedef triplet<R,MatriceCreuse<R> *,bool> P;
   typedef list<P> L;
   typedef L * RR;
   typedef Matrice_Creuse<R> * MM;

  static   RR f(const MM & a,const RR & b)  
  { 
    
    b->push_front(make_triplet<R,MatriceCreuse<R> *>(R(1.),a->A,false));
    return b;
  }
    
    
};

template<class R> 
struct Op2_ListCMMadd: public binary_function< list<triplet<R,MatriceCreuse<R> *,bool> > *,                                                                                              
                                               Matrice_Creuse<R> * ,
                                               list<triplet<R,MatriceCreuse<R> *,bool> > *>
{  //   .... + M
   typedef triplet<R,MatriceCreuse<R> *,bool> P;
   typedef list<P> L;
   typedef L * RR;
   typedef Matrice_Creuse<R> * MM;

  static   RR f(const RR & a,const MM & b)  
  { 
    
    a->push_back(make_triplet<R,MatriceCreuse<R> *,bool>(R(1.),b->A,false));
    return a;
  }
    
    
};

template<class R> 
struct Op2_ListMMadd: public binary_function< Matrice_Creuse<R> *,
                                              Matrice_Creuse<R> * ,
                                              list<triplet<R,MatriceCreuse<R> *,bool> > *>
{  //  M + M
   typedef triplet<R,MatriceCreuse<R> *,bool> P;
   typedef list<P> L;
   typedef L * RR;
   typedef Matrice_Creuse<R> * MM;

  static   RR f(const MM & a,const MM & b)  
  { 
    L * l=to(a);
    l->push_back(make_triplet<R,MatriceCreuse<R> *>(R(1.),b->A,false));
    return l;
  }
    
    
};

// Fin Add FH Houston --------

class MatrixInterpolation : public OneOperator { public:  

    class Op : public E_F0info { public:
       typedef pfes * A;
       Expression a,b,c; 
       // if c = 0 => a,b FESpace
       // if c != a FESpace et b,c KN_<double> 
       static const int n_name_param =5;
       static basicAC_F0::name_and_type name_param[] ;
        Expression nargs[n_name_param];
     bool arg(int i,Stack stack,bool a) const{ return nargs[i] ? GetAny<bool>( (*nargs[i])(stack) ): a;}
     long arg(int i,Stack stack,long a) const{ return nargs[i] ? GetAny<long>( (*nargs[i])(stack) ): a;}
     KN_<long>  arg(int i,Stack stack,KN_<long> a ) const{ return nargs[i] ? GetAny<KN_<long> >( (*nargs[i])(stack) ): a;}

       public:
       Op(const basicAC_F0 &  args,Expression aa,Expression bb) : a(aa),b(bb),c(0) {
         args.SetNameParam(n_name_param,name_param,nargs);  }
       Op(const basicAC_F0 &  args,Expression aa,Expression bb,Expression cc) : a(aa),b(bb),c(cc) {
         args.SetNameParam(n_name_param,name_param,nargs); } 
    };
    // interpolation(Vh,Vh)
   MatrixInterpolation() : OneOperator(atype<const MatrixInterpolation::Op *>(),
                                       atype<pfes *>(),
                                       atype<pfes *>()) {}
   // interpolation(Vh,xx,yy)
   MatrixInterpolation(int bidon) : OneOperator(atype<const MatrixInterpolation::Op *>(),
                                                atype<pfes *>(),atype<KN_<double> >(),
                                                atype<KN_<double> >()) {}
   
  
    E_F0 * code(const basicAC_F0 & args) const 
     { 
       if(args.size()==2)
       return  new Op(args,t[0]->CastTo(args[0]),
                           t[1]->CastTo(args[1]));
       else if(args.size()==3)
       return  new Op(args,t[0]->CastTo(args[0]),
                           t[1]->CastTo(args[1]),
                           t[2]->CastTo(args[2]));
       else CompileError("Bug in MatrixInterpolation code nb != 2 or 3 ????? bizarre" );
       return 0;
     }
};

basicAC_F0::name_and_type  MatrixInterpolation::Op::name_param[]= {
   {  "t", &typeid(bool)}, 
   {  "op", &typeid(long)},
   {  "inside",&typeid(bool)},
   {  "composante",&typeid(long)},
   {  "U2Vc",&typeid(KN_<long>)}

};





template<class R>
   class SetMatrix_Op : public E_F0mps { public:
       Expression a; 
       
       static  aType btype;
       static const int n_name_param =12; //  add nbiter FH 30/01/2007 11 -> 12 
       static basicAC_F0::name_and_type name_param[] ;
       Expression nargs[n_name_param];
       const OneOperator * precon;
       bool arg(int i,Stack stack,bool a) const{ return nargs[i] ? GetAny<bool>( (*nargs[i])(stack) ): a;}
       long arg(int i,Stack stack,long a) const{ return nargs[i] ? GetAny<long>( (*nargs[i])(stack) ): a;}

       public:
       SetMatrix_Op(const basicAC_F0 &  args,Expression aa) : a(aa) {
         args.SetNameParam(n_name_param,name_param,nargs);
         precon = 0; //  a changer 
         if ( nargs[3])
          {
           const  Polymorphic * op=  dynamic_cast<const  Polymorphic *>(nargs[3]);
           assert(op);
           precon = op->Find("(",ArrayOfaType(atype<KN<R>* >(),false)); // strange bug in g++ is R become a double
          }
         
       } 
       AnyType operator()(Stack stack)  const ;
    };


template<class R>
class SetMatrix : public OneOperator { public:  

 
  // SetMatrix() : OneOperator(atype<const typename SetMatrix<R>::Op *>(),atype<Matrice_Creuse<R> *>() ) {}
   SetMatrix() : OneOperator(SetMatrix_Op<R>::btype,atype<Matrice_Creuse<R> *>() ) {}
  
    E_F0 * code(const basicAC_F0 & args) const 
     { 
       return  new SetMatrix_Op<R>(args,t[0]->CastTo(args[0])); 
     }
};

template<class R>
       aType  SetMatrix_Op<R>::btype=0;

template <class R> 
basicAC_F0::name_and_type  SetMatrix_Op<R>::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)} // 11


};

template<class R>
AnyType SetMatrix_Op<R>::operator()(Stack stack)  const 
{
   Matrice_Creuse<R> *  A= GetAny<Matrice_Creuse<R> *>((*a)(stack));
   assert(A && A->A);
  long NbSpace = 50; 
  long itmax=0; 
  double eps=1e-6;
//  bool VF=false;
//  VF=isVF(op->largs);
 // assert(!VF); 
  double tgv = 1e30;
  bool VF=false;
  bool factorize=false;
  double tol_pivot=-1;
  double tol_pivot_sym=-1;
// type de matrice par default 
  TypeSolveMat tmat( TypeSolveMat::defaultvalue); 
  if(   tmat !=  TypeSolveMat::SparseSolver  )    
      tmat=TypeSolveMat::GMRES;