problem.cpp

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

	  R3 Pt(pi);
	  pa =a;
	  Ku.BF(Dop,Pt,fu);
	  MeshPointStack(stack)->set(T(Pt),Pt,Kv);
	  if (classoptm) {
	    if( oldopt) (*Op.optiexpK)(stack); // call old optim version 
	    else Op.optiexpK->eval(stack,iloop++,unvarexp); // new optim version 
	  }
	  if (!same) Kv.BF(Dop,Pt,fv);
	  int il=0;
	  for (BilinearOperator::const_iterator l=Op.v.begin();l!=Op.v.end();l++,il++)
	    {  // 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);
	      long jcomp= ll.first.first.first,jop=ll.first.first.second;
	      long icomp= ll.first.second.first,iop=ll.first.second.second;
	      
	      R ccc = copt ? *(copt[il]) : GetAny<R>(ll.second.eval(stack));
	      if ( copt && Kv.number <1)
		{
		  R cc  =  GetAny<R>(ll.second.eval(stack));
		  //cout << *(copt[il]) << " == " <<  cc << endl;
		  if ( ccc != cc) { 
		    cerr << cc << " != " << ccc << " => ";
		    cerr << "Sorry error in Optimization (a) add:  int2d(Th,optimize=0)(...)" << endl;
		    ExecError("In Optimized version "); }
		}
	      int fi=Kv.dfcbegin(icomp);
	      int li=Kv.dfcend(icomp);
	      int fj=Ku.dfcbegin(jcomp);
	      int lj=Ku.dfcend(jcomp);
	      ccc *= coef;
	      
	      // attention la fonction test donne la ligne 
	      //  et la fonction test est en second      
	      
	      for ( i=fi;  i<li;   i++ )  
		{ 
		  for ( j=fj;  j<lj;   j++ ) 
		    { 		      
		      R w_i =  fv(i,icomp,iop); 
		      R w_j =  fu(j,jcomp,jop);		      
		      mat(i,j) += ccc * w_i*w_j;
		    }
		}
	    }
	}
    
    else // int on edge ie 
     for (npi=0;npi<FIb.n;npi++) // loop on the integration point
      {
        pa =a;
	GQuadraturePoint<R2> pi( FIb[npi]);
	R3 NN= T.N(ie);
	  double mes=NN.norme();
	  NN/=mes;
	  double coef = 0.5*mes*pi.a; // correction 0.5 050109 FH
	  R3 Pt(T.PBord(ie,pi));
	  Ku.BF(Dop,Pt,fu);
	  if (!same) Kv.BF(Dop,Pt,fv);      
          MeshPointStack(stack)->set(T(Pt),Pt,Ku,label,NN,ie);
         if (classoptm) (*Op.optiexpK)(stack); // call optim version 
        
        
	  for ( i=0;  i<n;   i++ )  
	    { 
		RNM_ wi(fv(i,'.','.'));       
		for ( j=0;  j<m;   j++,pa++ ) 
		  { 
		      RNM_ wj(fu(j,'.','.'));
		      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);
			    
			    double w_i =  wi(ii.first,ii.second); 
			    double w_j =  wj(jj.first,jj.second);
			    
			    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;
			}
		  }
	    } 
      }
   
    
    if (Ku.Vh.Th(T) <1 && verbosity>100) {
      pa=mat.a;
      cout <<endl  << " Tet " << Ku.Vh.Th(T) << " =  " << T  << " " << nx << endl;
      for (int i=0;i<n;i++)
	{
	  cout << setw(2) << i << setw(4) << mat.ni[i] << " :";
	  for (int j=0;j<m;j++)
	    cout << setw(5)  << (*pa++) << " ";
	  cout << endl;
	} } 
    

  } 
  // xxxxxxxxxxxxxxxxx  modif a faire 
  template<class R> 
  void  Element_Op(MatriceElementairePleine<R,FESpace> & mat,const FElement & Ku,const FElement & Kv,double * p,int ie,int label,void *stack)
  {
    typedef  FElement::Element Element;
    MeshPoint mp= *MeshPointStack(stack);
    R ** copt = Stack_Ptr<R*>(stack,ElemMatPtrOffset);

  bool same = &Ku == & Kv;
  const Element & T  = Ku.T;
  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;
  long N= Kv.N;
  long M= Ku.N;
  
  
  
  
  
  const Opera &Op(*mat.bilinearform);
  bool classoptm = copt && Op.optiexpK;
  bool oldopt=1;  // juin 2007 FH ???? a voir 
  int  iloop=0;
  KN<bool> unvarexp(classoptm ? Op.optiexpK->sizevar() : 1);
  if (Ku.number<1 && verbosity/100 && verbosity % 10 == 2) 
     cout << "Element_Op P: copt = " << copt << " " << classoptm << endl;
    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);
  RNMK_ fv(p,n,N,lastop); //  the value for basic fonction
  RNMK_ fu(p+ (same ?0:n*N*lastop) ,m,M,lastop); //  the value for basic fonction
  for (i=0;i< nx;i++) 
    *pa++ = 0.; 
  if (ie<0)    
    for (npi=0;npi<FI.n;npi++) // loop on the integration point
      {
        QuadraturePoint pi(FI[npi]);
        R coef = T.area*pi.a;
        R2 Pt(pi);
        pa =a;
        Ku.BF(Dop,Pt,fu);
        MeshPointStack(stack)->set(T(Pt),Pt,Kv);
        if (classoptm) {
	    if( oldopt) (*Op.optiexpK)(stack); // call old optim version 
	    else Op.optiexpK->eval(stack,iloop++,unvarexp); // new optim version 
	}
        if (!same) Kv.BF(Dop,Pt,fv);      
	  int il=0;
	  for (BilinearOperator::const_iterator l=Op.v.begin();l!=Op.v.end();l++,il++)
	    {  // 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);
	      long jcomp= ll.first.first.first,jop=ll.first.first.second;
	      long icomp= ll.first.second.first,iop=ll.first.second.second;
	      
	      R ccc = copt ? *(copt[il]) : GetAny<R>(ll.second.eval(stack));
	      if ( copt && Kv.number <1)
		{
		  R cc  =  GetAny<R>(ll.second.eval(stack));
		  //cout << *(copt[il]) << " == " <<  cc << endl;
		  if ( ccc != cc) { 
		    cerr << cc << " != " << ccc << " => ";
		    cerr << "Sorry error in Optimization (a) add:  int2d(Th,optimize=0)(...)" << endl;
		    ExecError("In Optimized version "); }
		}
	      int fi=Kv.dfcbegin(icomp);
	      int li=Kv.dfcend(icomp);
	      int fj=Ku.dfcbegin(jcomp);
	      int lj=Ku.dfcend(jcomp);
	      fi=0,fj=0;
	      li=n,lj=m;
	      ccc *= coef;
	      
	      // attention la fonction test donne la ligne 
	      //  et la fonction test est en second      
	      for ( i=fi;  i<li;   i++ )  
		{ 
		  for ( j=fj;  j<lj;   j++ ) 
		    { 		      
		      R w_i =  fv(i,icomp,iop); 
		      R w_j =  fu(j,jcomp,jop);		      
		      mat(i,j) += ccc * w_i*w_j;
		    }
		}
	    }

	/*
        for ( i=0;  i<n;   i++ )  
          { 
            
            // attention la fonction test donne la ligne 
            //  et la fonction test est en second      
            
            RNM_ wi(fv(i,'.','.'));         
            for ( j=0;  j<m;   j++,pa++ ) 
              { 
                RNM_ wj(fu(j,'.','.'));
                int il=0;
                for (BilinearOperator::const_iterator l=Op.v.begin();l!=Op.v.end();l++,il++)
                  {  // 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);
                    R w_i =  wi(ii.first,ii.second); 
                    R w_j =  wj(jj.first,jj.second);
                    R ccc = copt ? *(copt[il]) : GetAny<R>(ll.second.eval(stack));
                if ( copt && Kv.number <1)
                 {
                     R cc  =  GetAny<R>(ll.second.eval(stack));
                     //cout << *(copt[il]) << " == " <<  cc << endl;
                     if ( ccc != cc) { 
                        cerr << cc << " != " << ccc << " => ";
                       cerr << "Sorry error in Optimization (a) add:  int2d(Th,optimize=0)(...)" << endl;
                       ExecError("In Optimized version "); }
                 }
                    
                    
                    *pa += coef * ccc * w_i*w_j;
                  }
              }
          }
	*/
      }
  else // int on edge ie 
    for (npi=0;npi<FIb.n;npi++) // loop on the integration point
      {
        pa =a;
        QuadratureFormular1dPoint pi( FIb[npi]);
        R2 E=T.Edge(ie);
        double le = sqrt((E,E));
        double coef = le*pi.a;
        double sa=pi.x,sb=1-sa;
        R2 PA(TriangleHat[VerticesOfTriangularEdge[ie][0]]),
          PB(TriangleHat[VerticesOfTriangularEdge[ie][1]]);
        R2 Pt(PA*sa+PB*sb ); //  
        Ku.BF(Dop,Pt,fu);
        if (!same) Kv.BF(Dop,Pt,fv);      
        // int label=-999999; // a passer en argument 
        MeshPointStack(stack)->set(T(Pt),Pt,Kv,label,R2(E.y,-E.x)/le,ie);
        if (classoptm) (*Op.optiexpK)(stack); // call optim version 
	int il=0;
	for (BilinearOperator::const_iterator l=Op.v.begin();l!=Op.v.end();l++,il++)
	  {  // 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);
	    long jcomp= ll.first.first.first,jop=ll.first.first.second;
	    long icomp= ll.first.second.first,iop=ll.first.second.second;

	    
	    R ccc = copt ? *(copt[il]) : GetAny<R>(ll.second.eval(stack));
	    if ( copt && Kv.number <1)
	      {
		R cc  =  GetAny<R>(ll.second.eval(stack));
		//cout << *(copt[il]) << " == " <<  cc << endl;
		if ( ccc != cc) { 
		  cerr << cc << " != " << ccc << " => ";
		  cerr << "Sorry error in Optimization (a) add:  int2d(Th,optimize=0)(...)" << endl;
		  ExecError("In Optimized version "); }
	      }
	    int fi=Kv.dfcbegin(icomp);
	    int li=Kv.dfcend(icomp);
	    int fj=Ku.dfcbegin(jcomp);
	    int lj=Ku.dfcend(jcomp);
	    ccc *= coef;
	    
	    // attention la fonction test donne la ligne 
	    //  et la fonction test est en second      
	    
	    for ( i=fi;  i<li;   i++ )  
	      { 
		  for ( j=fj;  j<lj;   j++ ) 
		    { 		      
		      R w_i =  fv(i,icomp,iop); 
		      R w_j =  fu(j,jcomp,jop);		      
		      mat(i,j) += ccc * w_i*w_j;
		    }
	      }
	  }
	  
	/*        
        for ( i=0;  i<n;   i++ )  
         // if (onWhatIsEdge[ie][Kv.DFOnWhat(i)]) // juste the df on edge bofbof generaly wrong FH dec 2003
            { 
              RNM_ wi(fv(i,'.','.'));       
              for ( j=0;  j<m;   j++,pa++ ) 
                { 
                  RNM_ wj(fu(j,'.','.'));
                  int il=0;
                  for (BilinearOperator::const_iterator l=Op.v.begin();l!=Op.v.end();l++,il++)
                   // if (onWhatIsEdge[ie][Kv.DFOnWhat(j)]) // juste the df on edge bofbof generaly wrong FH dec 2003
                      {       
                        BilinearOperator::K ll(*l);
                        pair<int,int> jj(ll.first.first),ii(ll.first.second);
                        
                        double w_i =  wi(ii.first,ii.second); 
                        double w_j =  wj(jj.first,jj.second);
                       // 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;  FH dec 2003
	 */
      }
  
  
  /*  pa=a;
      if (Ku.Vh.Th(T) >=0 ) {
      cout <<endl  << " Triangle " << Ku.Vh.Th(T) << " =  "<<  T[0] << ", " << T[1] << ", " << T[2] << " " << nx << endl;
      for (int i=0;i<n;i++)
      {
      cout << setw(2) << i << setw(4) << mat.ni[i] << " :";
      for (int j=0;j<m;j++)
      cout << setw(5)  << (*pa++) << " ";
      cout << endl;
      } } 
  */ 
  *MeshPointStack(stack) = mp;
  }  
  
  
  
 template<class R>
 void  Element_Op(MatriceElementaireSymetrique<R,FESpace3> & mat,const FElement3 & Ku,double * p,int ie,int label, void * stack)
  {
   typedef FESpace3 FESpace;
   typedef typename FESpace3::Mesh Mesh;
   typedef Mesh *pmesh ;
   typedef typename Mesh::Element Element;
    MeshPoint mp= *MeshPointStack(stack);
    R ** copt = Stack_Ptr<R*>(stack,ElemMatPtrOffset);
    const Element & T  = Ku.T;

    const GQuadratureFormular<R3> & FI = mat.FIT;
    const GQuadratureFormular<R2> & 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+1)/2;