problem.cpp

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

      {
	assert(b->b->v.size()==(size_t) n_where_in_stack_opt);
	for (int i=0;i<n_where_in_stack_opt;i++)
	  {
	    int offset=b->b->where_in_stack_opt[i];
	    assert(offset>10);
	    where_in_stack[i]= static_cast<R *>(static_cast<void *>((char*)stack+offset));
	    *(where_in_stack[i])=0;
	  }
	
	
	if(&optiexp0) 
	  optiexp0(stack); 
	KN<bool> ok(b->b->v.size());
	{  //   remove the zero coef in the liste 
	  // R zero=R();  
	  int il=0;
	  for (BilinearOperator::const_iterator l=b->b->v.begin();l!=b->b->v.end();l++,il++)
	    ok[il] =  ! (b->b->mesh_indep_stack_opt[il] && ( norm(*(where_in_stack[il])) < 1e-100 ) );
	}
	BilinearOperator b_nozer(*b->b,ok); 
	if (verbosity % 10 > 3 ) 
	  cout << "   -- nb term in bilinear form  (!0) : " << b_nozer.v.size() 
	       << "  total " << n_where_in_stack_opt << endl;
	
	if ( (verbosity/100) % 10 >= 2)   
	  { 
	    int il=0;
	    
	    for (BilinearOperator::const_iterator l=b->b->v.begin();l!=b->b->v.end();l++,il++)
	      cout << il << " coef (" << l->first << ") = " << *(where_in_stack[il]) 
		   << " offset=" << b->b->where_in_stack_opt[il] 
		   << " dep mesh " << l->second.MeshIndependent() << b->b->mesh_indep_stack_opt[il] << endl;
	  }
      }
    Stack_Ptr<R*>(stack,ElemMatPtrOffset) =where_in_stack;
    void *paramate=stack;
    pair<void *,double *> parammatElement_OpVF;  
    parammatElement_OpVF.first = stack;
    parammatElement_OpVF.second= & binside;
    
    if (verbosity >3) 
      if (all) cout << " all " << endl ;
      else cout << endl;
    if(VF) {
      if(&Uh != &Vh || sym)
	ExecError("To Day in bilinear form with discontinous Galerkin:   \n"
		  "  test or unkown function must be  defined on the same FEspace, \n"
		  "  and the matrix is not symmetric. \n" 
		  " To do other case in a future (F. Hecht) dec. 2003 ");
      
      matep= new MatriceElementairePleine<R,FESpace>(Uh,VF,FIT,FIE);
      matep->faceelement = Element_OpVF;   
      paramate= &parammatElement_OpVF;            
    }
    else if (sym) {
      mates= new MatriceElementaireSymetrique<R,FESpace>(Uh,FIT,FIE);
      mates->element = Element_Op<R>;               
    }
    else {
      matep= new MatriceElementairePleine<R,FESpace>(Uh,Vh,FIT,FIE);
      matep->element = Element_Op<R>;               
    }
    MatriceElementaireFES<R,FESpace> & mate(*( sym? (MatriceElementaireFES<R,FESpace> *)mates : (MatriceElementaireFES<R,FESpace> *) matep));
    
    
    mate.bilinearform=b->b;
    
    Check(*mate.bilinearform,mate.Uh.N,mate.Vh.N);
    
    if (di.kind == CDomainOfIntegration::int1d )
      {
        for( int e=0;e<Th.neb;e++)
          {
            if (all || setoflab.find(Th.bedges[e].lab) != setoflab.end())   
              {                  
                int ie,i =Th.BoundaryElement(e,ie);
                A += mate(i,ie,Th.bedges[e].lab,stack);  
                if(sptrclean) sptrclean=sptr->clean(); // modif FH mars 2006  clean Ptr
		
              }
          }
      }
    else if (di.kind == CDomainOfIntegration::intalledges)
      {
        for (int i=0;i< Th.nt; i++) 
          {
            if ( all || setoflab.find(Th[i].lab) != setoflab.end())
	      for (int ie=0;ie<3;ie++)   
              A += mate(i,ie,Th[i].lab,paramate); 
	    if(sptrclean) sptrclean=sptr->clean(); // modif FH mars 2006  clean Ptr
	    
          }
	
      }      
    else if (di.kind == CDomainOfIntegration::intallVFedges)
      {
	cerr << " a faire intallVFedges " << endl;
	ffassert(0);
        for (int i=0;i< Th.nt; i++) 
          {
            if ( all || setoflab.find(Th[i].lab) != setoflab.end())
             for (int ie=0;ie<3;ie++)   
	       A += mate(i,ie,Th[i].lab,paramate);  
	    if(sptrclean) sptrclean=sptr->clean(); // modif FH mars 2006  clean Ptr
	    
          }
	
      }      
    else if (di.kind == CDomainOfIntegration::int2d )
      {
        for (int i=0;i< Th.nt; i++) 
          {
            if ( all || setoflab.find(Th[i].lab) != setoflab.end())  
              A += mate(i,-1,Th[i].lab,stack);   
            if(sptrclean) sptrclean=sptr->clean(); // modif FH mars 2006  clean Ptr

            // AA += mate;
          }
      } 
    else 
      InternalError(" kind of CDomainOfIntegration unkown");
    
    if (where_in_stack) delete [] where_in_stack;
    delete &mate;
 }
// --------- FH 120105
// 3d
 template<class R>
 void AssembleBilinearForm(Stack stack,const FESpace3::Mesh & Th,const FESpace3 & Uh,const FESpace3 & Vh,bool sym,
			   MatriceCreuse<R>  & A, const  FormBilinear * b  )
   
 {
   typedef FESpace3 FESpace;
   typedef FESpace3::Mesh Mesh;
   typedef Mesh *pmesh ;
   StackOfPtr2Free * sptr = WhereStackOfPtr2Free(stack);
   bool sptrclean=true;
   const CDomainOfIntegration & di= *b->di;
   ffassert(di.d==3);
   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);
    const GQuadratureFormular<R3> & FIV = di.FIV(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::int2d==kind) cout << "  -- boundary int border ( nQP: "<< FIT.n << ") ,"  ;
      else  if (CDomainOfIntegration::intalledges==kind) cout << "  -- boundary int all edges ( nQP: "<< FIT.n << "),"  ;
      else  if (CDomainOfIntegration::intallVFedges==kind) cout << "  -- boundary int all VF edges nQP: ("<< FIE.n << ")," ;
      else cout << "  --  int3d   (nQP: "<< FIV.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)
      {
	assert(b->b->v.size()==(size_t) n_where_in_stack_opt);
	for (int i=0;i<n_where_in_stack_opt;i++)
	  {
	    int offset=b->b->where_in_stack_opt[i];
	    assert(offset>10);
	    where_in_stack[i]= static_cast<R *>(static_cast<void *>((char*)stack+offset));
	    *(where_in_stack[i])=0;
	  }
	
	
	if(&optiexp0) 
	  optiexp0(stack); 
	KN<bool> ok(b->b->v.size());
	{  //   remove the zero coef in the liste 
	  // R zero=R();  
	  int il=0;
	  for (BilinearOperator::const_iterator l=b->b->v.begin();l!=b->b->v.end();l++,il++)
	    ok[il] =  ! (b->b->mesh_indep_stack_opt[il] && ( norm(*(where_in_stack[il])) < 1e-100 ) );
	}
	BilinearOperator b_nozer(*b->b,ok); 
	if (verbosity % 10 > 3 ) 
	  cout << "   -- nb term in bilinear form  (!0) : " << b_nozer.v.size() 
	       << "  total " << n_where_in_stack_opt << endl;
	
	if ( (verbosity/100) % 10 >= 2)   
	  { 
	    int il=0;
	    
	    for (BilinearOperator::const_iterator l=b->b->v.begin();l!=b->b->v.end();l++,il++)
	      cout << il << " coef (" << l->first << ") = " << *(where_in_stack[il]) 
		   << " offset=" << b->b->where_in_stack_opt[il] 
		   << " dep mesh " << l->second.MeshIndependent() << b->b->mesh_indep_stack_opt[il] << endl;
	  }
      }
    Stack_Ptr<R*>(stack,ElemMatPtrOffset) =where_in_stack;
    void *paramate=stack;
    pair<void *,double *> parammatElement_OpVF;  
    parammatElement_OpVF.first = stack;
    parammatElement_OpVF.second= & binside;
    
    if (verbosity >3) 
      if (all) cout << " all " << endl ;
      else cout << endl;
    if(VF) {
      if(&Uh != &Vh || sym)
	ExecError("To Day in bilinear form with discontinous Galerkin:   \n"
		  "  test or unkown function must be  defined on the same FEspace, \n"
		  "  and the matrix is not symmetric. \n" 
		  " To do other case in a future (F. Hecht) dec. 2003 ");
      
      matep= new MatriceElementairePleine<R,FESpace>(Uh,VF,FIV,FIT);
      matep->faceelement = Element_OpVF;   
      paramate= &parammatElement_OpVF;            
    }
    else if (sym) {
      mates= new MatriceElementaireSymetrique<R,FESpace>(Uh,FIV,FIT);
      mates->element = Element_Op<R>;               
    }
    else {
      matep= new MatriceElementairePleine<R,FESpace>(Uh,Vh,FIV,FIT);
      matep->element = Element_Op<R>;               
    }
    MatriceElementaireFES<R,FESpace> & mate(*( sym? (MatriceElementaireFES<R,FESpace> *)mates : (MatriceElementaireFES<R,FESpace> *) matep));
    
    
    mate.bilinearform=b->b;
    
    Check(*mate.bilinearform,mate.Uh.N,mate.Vh.N);
    
    if (di.kind == CDomainOfIntegration::int2d )
      {
        for( int e=0;e<Th.nbe;e++)
          {
            if (all || setoflab.find(Th.be(e).lab) != setoflab.end())   
              {                  
                int ie,i =Th.BoundaryElement(e,ie);
                A += mate(i,ie,Th.be(e).lab,stack);  
                if(sptrclean) sptrclean=sptr->clean(); // modif FH mars 2006  clean Ptr
		
              }
          }
      }
    else if (di.kind == CDomainOfIntegration::intalledges)
      {
        for (int i=0;i< Th.nt; i++) 
          {
            if ( all || setoflab.find(Th[i].lab) != setoflab.end())
	      for (int ie=0;ie<3;ie++)   
              A += mate(i,ie,Th[i].lab,paramate); 
	    if(sptrclean) sptrclean=sptr->clean(); // modif FH mars 2006  clean Ptr
	    
          }
	
      }      
    else if (di.kind == CDomainOfIntegration::intallVFedges)
      {
	cerr << " a faire intallVFedges " << endl;
	ffassert(0);
        for (int i=0;i< Th.nt; i++) 
          {
            if ( all || setoflab.find(Th[i].lab) != setoflab.end())
             for (int ie=0;ie<3;ie++)   
	       A += mate(i,ie,Th[i].lab,paramate);  
	    if(sptrclean) sptrclean=sptr->clean(); // modif FH mars 2006  clean Ptr
	    
          }
	
      }      
    else if (di.kind == CDomainOfIntegration::int3d )
      {
        for (int i=0;i< Th.nt; i++) 
          {
            if ( all || setoflab.find(Th[i].lab) != setoflab.end())  
              A += mate(i,-1,Th[i].lab,stack);   
            if(sptrclean) sptrclean=sptr->clean(); // modif FH mars 2006  clean Ptr
	    
            // AA += mate;
          }
      } 
    else 
      InternalError(" kind of CDomainOfIntegration unkown");
    
    if (where_in_stack) delete [] where_in_stack;
    delete &mate;
 }

// end 3d

// --------- FH 170605  

    template<class R> 
    void  AddMatElem(map<pair<int,int>, R > & A,const Mesh & Th,const BilinearOperator & Op,bool sym,int it,  int ie,int label,
		     const FESpace & Uh,const FESpace & Vh,
		     const QuadratureFormular & FI,
		     const QuadratureFormular1d & FIb,
		     double *p,   void *stack, bool intmortar=false)
    {
	MeshPoint mp= *MeshPointStack(stack);
	R ** copt = Stack_Ptr<R*>(stack,ElemMatPtrOffset);
	const Mesh & Thu(Uh.Th);
	const Mesh & Thv(Vh.Th);
	
	bool same = &Uh == & Vh;
	const Triangle & T  = Th[it];
	long npi;
	long i,j;
	bool classoptm = copt && Op.optiexpK;
	assert(Op.MaxOp() <last_operatortype);
	//
	
	
	KN<bool> Dop(last_operatortype);
	Op.DiffOp(Dop);  
	int lastop=1+Dop.last(binder1st<equal_to<bool> >(equal_to<bool>(),true));
	//assert(lastop<=3);
	
	if (ie<0)    
	    for (npi=0;npi<FI.n;npi++) // loop on the integration point
	    {
		QuadraturePoint pi(FI[npi]);
		double coef = T.area*pi.a;
		R2 Pt(pi),Ptu,Ptv;
		R2 P(T(Pt));
		bool outsideu,outsidev;
		// ici trouve le T
		int iut=0,ivt=0;
		const Triangle * tu,*tv;
		if(&Th == & Thu )
		{
		    tu =&T;
		    Ptu=Pt;
		}
		else
		{
		    tu= Thu.Find(P,Ptu,outsideu);
		    if( !tu ||  outsideu) { 
			if(verbosity>100) cout << " On a pas trouver (u) " << P << " " << endl; 
			continue;}}
		if(same)
		{
		    tv=tu;
		    outsidev=outsideu;
		    Ptv=Ptu;
		}
		else
		{
		    if(&Th == & Thv )
		    {
			tv =&T;
			Ptv=Pt;
		    }
		    else