problem.cpp

来自「FreeFem++可以生成高质量的有限元网格。可以用于流体力学」· C++ 代码 · 共 1,861 行 · 第 1/5 页

// -*- 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>
using namespace std;

#include "rgraph.hpp"
#include "error.hpp"
#include "AFunction.hpp"

//#include "lex.hpp"
#include "MatriceCreuse_tpl.hpp"
#include "Mesh3dn.hpp"
#include "MeshPoint.hpp"
#include "lgfem.hpp"
#include "lgmesh3.hpp"
#include "lgsolver.hpp"
#include "problem.hpp"
#include <set>



basicAC_F0::name_and_type  CDomainOfIntegration::name_param[]= {
    { "qft", &typeid(const Fem2D::QuadratureFormular *)},
    { "qfe", &typeid(const Fem2D::QuadratureFormular1d *)},
    { "qforder",&typeid(long)},
    { "qfnbpT",&typeid(long)},
    { "qfnbpE",&typeid(long)},
    { "optimize",&typeid(bool)},
    { "binside",&typeid(double)},
    { "mortar",&typeid(bool)}
};


basicAC_F0::name_and_type  Problem::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 )},
  {  "strategy",&typeid(long )},
  {  "save",&typeid(string* )},
  {  "cadna",&typeid(KN<double>*)},
  {  "tolpivot", &typeid(double)},
  {  "tolpivotsym", &typeid(double)},
  {  "nbiter", &typeid(long)} // 12 
  
};



namespace Fem2D {




void Check(const Opera &Op,int N,int  M)
 {
   int err=0;
   for (BilinearOperator::const_iterator l=Op.v.begin();l!=Op.v.end();l++)
     {  // attention la fonction test donne la ligne 
       //  et la fonction test est en second      
       BilinearOperator::K ll(*l);
       pair<int,int> jj(ll.first.first),ii(ll.first.second);
       if (ii.first <0 || ii.first >= M) err++;
       if (jj.first <0 || jj.first >= N) err++;
       
     }
   if (err) {
     cout << "Check Bilinear Operator" << N << " " << M << endl;
     for (BilinearOperator::const_iterator l=Op.v.begin();l!=Op.v.end();l++)
       {  // attention la fonction test donne la ligne 
         //  et la fonction test est en second      
         BilinearOperator::K ll(*l);
         pair<int,int> jj(ll.first.first),ii(ll.first.second);
         cout << " +  " << jj.first << " " << jj.second << "*" << ii.first << " " << ii.second << endl;
          }
     ExecError("Check BilinearOperator N M");
   }
 }
 void Check(const  BC_set * bc,int N)
 {
   int err=0;
   int kk=bc->bc.size();
   for (int k=0;k<kk;k++)
     {
       pair<int,Expression> xx=bc->bc[k];
       if (xx.first >= N) { 
         err++;
         cerr << " Sorry : just " << N << " componant in FE space \n"
              << "   and Boundary condition refere to " << xx.first << "+1 componant " << endl;
       }
     }
   if (err) 
     ExecError("Incompatibility beetwen  boundary condition  and FE space");
 }
 
  void Check(const  Ftest * fl,int N)
 {
   assert(fl);
   int err=0;
   Ftest::const_iterator kk= fl->v.end(),k;
   int ii=0;
   for (k=fl->v.begin();k<kk;k++)
     { 
       ii++;
       int j=k->first.first;
       if (  j >= N) { 
         err++;
         cerr << " Sorry : just " << N << " componant in FE space \n"
              << " and linear var form  refere to " << j << "+1 componant (part " << ii << ")" << endl;
       }
     }
   if (err) 
     ExecError("Incompatibility beetwen linear varf  and FE space");
 }
 //---------------------------------------------------------------------------------------
  template<class R>
  void  Element_OpVF(MatriceElementairePleine<R,FESpace3> & mat,
		     const FElement3 & Ku,const FElement3 & KKu,
		     const FElement3 & Kv,const FElement3 & KKv,
		     double * p,int ie,int iie, int label,void *bstack)
  {
    ffassert(0); 
  }
  
  template<class R>
  void  Element_OpVF(MatriceElementairePleine<R,FESpace> & mat,
			const FElement & Ku,const FElement & KKu,
			const FElement & Kv,const FElement & KKv,
			double * p,int ie,int iie, int label,void *bstack)
  {
    
    pair<void *,double *> * bs=static_cast<pair<void *,double *> *>(bstack);   
    void * stack= bs->first;
    double binside = *bs->second; // truc FH pour fluide de grad2 (decentrage bizard)
    assert(mat.onFace); //   Finite Volume or discontinous Galerkine 
    assert(ie>=0 && ie < 3); //  int on edge 
    MeshPoint mp= *MeshPointStack(stack);
    R ** copt = Stack_Ptr<R*>(stack,ElemMatPtrOffset);
    
    bool same = &Ku == & Kv;
    assert(same); 
    const Triangle & T  = Ku.T;
    int nTonEdge =  &Ku == &KKu ? 1 : 2; 
    double cmean = 1./nTonEdge;
    
    throwassert(&T == &Kv.T);  
  // const QuadratureFormular & FI = mat.FIT;
    const QuadratureFormular1d & FIb = mat.FIE;
    long npi;
    R *a=mat.a;
    R *pa=a;
    long i,j;
    long n= mat.n,m=mat.m,nx=n*m;
    assert(nx<=mat.lga);
    long N= Kv.N;
    long M= Ku.N;
    
    long mu=Ku.NbDoF();
    long mmu=KKu.NbDoF();
    long nv=Kv.NbDoF();
    long nnv=Kv.NbDoF();
    assert(mu==mmu && nv == nnv) ; 
    
    
    
    const Opera &Op(*mat.bilinearform);
    bool classoptm = copt && Op.optiexpK;
    //  if (Ku.number<1 && verbosity/100 && verbosity % 10 == 2) 
    if (Ku.number<1 && ( verbosity > 1 ) )
      cout << "Element_OpVF P: copt = " << copt << " " << classoptm << " binside (For FH) =" << binside <<endl;
  
    
    KN<bool> Dop(last_operatortype); //  sinon ca plate bizarre 
    Op.DiffOp(Dop);  
    int lastop=1+Dop.last(binder1st<equal_to<bool> >(equal_to<bool>(),true));
    //assert(lastop<=3);
    int lffv = nv*N*last_operatortype;  
    int lffu = mu*M*last_operatortype;
    int loffset =  same ? 0 :  (nv+nnv)*N*last_operatortype;
    
    RNMK_ fv(p,nv,N,lastop); //  the value for basic fonction in K
    RNMK_ ffv(p + lffv ,nnv,N,lastop); //  the value for basic fonction in KK
    RNMK_ fu(  (double*) fv   + loffset  ,mu,M,lastop); //  the value for basic fonction
    RNMK_ ffu( (double*) fu  + lffu  ,mmu,M,lastop); //  the value for basic fonction
    
    R2 E=T.Edge(ie);
    double le = sqrt((E,E));
    R2 PA(TriangleHat[VerticesOfTriangularEdge[ie][0]]),
      PB(TriangleHat[VerticesOfTriangularEdge[ie][1]]),
      PC(TriangleHat[OppositeVertex[ie]]);
    // warning the to edge are in opposite sens         
    R2 PP_A(TriangleHat[VerticesOfTriangularEdge[iie][1]]),
      PP_B(TriangleHat[VerticesOfTriangularEdge[iie][0]]),
      PP_C(TriangleHat[OppositeVertex[ie]]);
    R2 Normal(E.perp()/-le); 
    for (npi=0;npi<FIb.n;npi++) // loop on the integration point
      {
        pa =a;
        QuadratureFormular1dPoint pi( FIb[npi]);        
        double coef = le*pi.a;
        double sa=pi.x,sb=1-sa;
        R2 Pt(PA*sa+PB*sb ); //
        R2 PP_t(PP_A*sa+PP_B*sb ); //  
        if (binside) {
	  Pt   = (1-binside)*Pt + binside*PC; 
	  PP_t  = (1-binside)*PP_t + binside*PP_C; }
        Ku.BF(Dop,Pt,fu);
        KKu.BF(Dop,PP_t,ffu);
        if (!same) { Kv.BF(Dop,Pt,fv); KKv.BF(Dop,PP_t,ffv); }     
        // int label=-999999; // a passer en argument 
        MeshPointStack(stack)->set(T(Pt),Pt,Kv,label, Normal,ie);
        if (classoptm) (*Op.optiexpK)(stack); // call optim version 
        
        
        for ( i=0;  i<n;   i++ )  
            { 
              int ik= mat.nik[i];
              int ikk=mat.nikk[i]; 
               
              RNM_ wi(fv(Max(ik,0),'.','.'));   
              RNM_ wwi(ffv(Max(ikk,0),'.','.'));   
                  
              for ( j=0;  j<m;   j++,pa++ ) 
                { 
                  int jk= mat.njk[j];
                  int jkk=mat.njkk[j];
                    
                  RNM_ wj(fu(Max(jk,0),'.','.'));
                  RNM_ wwj(ffu(Max(jkk,0),'.','.'));
                  
                  int il=0;
                  for (BilinearOperator::const_iterator l=Op.v.begin();l!=Op.v.end();l++,il++)
                      {       
                        BilinearOperator::K ll(*l);
                        pair<int,int> jj(ll.first.first),ii(ll.first.second);
                        int iis = ii.second, jjs=jj.second;
                        
                        int iicase  = iis / last_operatortype;
                        int jjcase  = jjs / last_operatortype;
                        
                         iis %= last_operatortype;
                         jjs %= last_operatortype;
                        double w_i=0,w_j=0,ww_i=0,ww_j=0;
                        
                        if(ik>=0) w_i =   wi(ii.first,iis ); 
                        if(jk>=0) w_j =   wj(jj.first,jjs );
                        
                        if( iicase>0 && ikk>=0) ww_i =  wwi(ii.first,iis ); 
                        if( jjcase>0 && jkk>=0) ww_j =  wwj(jj.first,jjs );
                       
                        
                        if       (iicase==Code_Jump) w_i = ww_i-w_i; // jump
                        else  if (iicase==Code_Mean) { 
                            
                            w_i = cmean*  (w_i + ww_i );} // average
                        else  if (iicase==Code_OtherSide) w_i = ww_i;  // valeur de autre cote
                        
                        if      (jjcase==Code_Jump) w_j = ww_j-w_j; // jump
                        else if (jjcase==Code_Mean) w_j = cmean*  (w_j +ww_j ); // average
                        else if (jjcase==Code_OtherSide) w_j = ww_j;  //  valeur de l'autre cote    
                        
                       // R ccc = GetAny<R>(ll.second.eval(stack));
                       
                        R ccc = copt ? *(copt[il]) : GetAny<R>(ll.second.eval(stack));
                       if ( copt && Kv.number <1)
                        {
                         R cc  =  GetAny<R>(ll.second.eval(stack));
                         if ( ccc != cc) { 
                          cerr << cc << " != " << ccc << " => ";
                         cerr << "Sorry error in Optimization (b) add:  int2d(Th,optimize=0)(...)" << endl;
                         ExecError("In Optimized version "); }
                 }
                         *pa += coef * ccc * w_i*w_j;
                      }
                }
            } 
         // else pa += m;
      }
  
  
     pa=a;
     if (verbosity > 55 && (Ku.number <=0 || KKu.number <=0 )) { 
       cout <<endl  << " edge between " << Ku.number << " , " <<  KKu.number   << " =  "<<  T[0] << ", " << T[1] << ", " << T[2] << " " << nx << endl;
       cout << " K u, uu =  " << Ku.number << " " << KKu.number << " " <<  " K v, vv =  " << Kv.number << " " << KKv.number << " " <<endl; 
       for (int i=0;i<n;i++)
	 {
	   cout << setw(2) << i << setw(4) << mat.ni[i] <<  setw(4) << mat.nik[i] << setw(4) << mat.nikk[i]  <<  " :";
	   for (int j=0;j<m;j++)
	     cout << setw(5)  << (*pa++) << " ";
	   cout << endl;
	 } } 
     
     *MeshPointStack(stack) = mp;
  }  

 //--------------------------------------------------------------------------------------
 
// --------- FH 120105
 template<class R>
 void AssembleBilinearForm(Stack stack,const Mesh & Th,const FESpace & Uh,const FESpace & Vh,bool sym,
			   MatriceCreuse<R>  & A, const  FormBilinear * b  )
   
 {
   StackOfPtr2Free * sptr = WhereStackOfPtr2Free(stack);
   bool sptrclean=true;
   const CDomainOfIntegration & di= *b->di;
   const Mesh * pThdi = GetAny<pmesh>( (* di.Th)(stack));
   if ( pThdi != &Th) { 
     ExecError("No way to compute bilinear form with integrale of on mesh \n"
	       "  test  or unkown function  defined on an other mesh! sorry to hard.   ");
   }
   SHOWVERB(cout << " FormBilinear " << endl);
    MatriceElementaireSymetrique<R,FESpace> *mates =0;
    MatriceElementairePleine<R,FESpace> *matep =0;
    const bool useopt=di.UseOpt(stack);    
    double binside=di.binside(stack);
    
    const vector<Expression>  & what(di.what);             
    CDomainOfIntegration::typeofkind  kind = di.kind;
    set<int> setoflab;
    bool all=true; 
    const QuadratureFormular1d & FIE = di.FIE(stack);
    const QuadratureFormular & FIT = di.FIT(stack);
    bool VF=b->VF();  // finite Volume or discontinous Galerkin
    if (verbosity>2) cout << "  -- discontinous Galerkin  =" << VF << " size of Mat =" << A.size()<< " Bytes\n";
    if (verbosity>3) 
      if (CDomainOfIntegration::int1d==kind) cout << "  -- boundary int border ( nQP: "<< FIE.n << ") ,"  ;
      else  if (CDomainOfIntegration::intalledges==kind) cout << "  -- boundary int all edges ( nQP: "<< FIE.n << "),"  ;
      else  if (CDomainOfIntegration::intallVFedges==kind) cout << "  -- boundary int all VF edges nQP: ("<< FIE.n << ")," ;
      else cout << "  --  int    (nQP: "<< FIT.n << " ) in "  ;
    for (size_t i=0;i<what.size();i++)
      {
	long  lab  = GetAny<long>( (*what[i])(stack));
	setoflab.insert(lab);
	if ( verbosity>3) cout << lab << " ";
	all=false;
      }
    if (verbosity>3) cout <<" Optimized = "<< useopt << ", ";
    const E_F0 & optiexp0=*b->b->optiexp0;
    
    int n_where_in_stack_opt=b->b->where_in_stack_opt.size();
    R** where_in_stack =0;
    if (n_where_in_stack_opt && useopt)
      where_in_stack = new R * [n_where_in_stack_opt];
    if (where_in_stack)