doperator.hpp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

#include <queue>
// class for 1 linear form
// just copy a array 

//  Warning this class are use at compilating time
//  ----------------------------------------------
using Fem2D::operatortype;
using Fem2D::op_id;
using Fem2D::op_dx;
using Fem2D::op_dy;
using Fem2D::last_operatortype;

template<class T> 
inline  T * NewCopy(const  T * const o,int n)
{ 
  int  m=n;
  T * c = new T [n];
  for (int i=0;i<m;i++)
    c[i]=o[i];
  return c;
}
template<class T> 
inline  T * NewCopy(const  T * const o,int n,int m)
{ 
  throwassert(m<=n);
  T * c = new T [n];
  for (int i=0;i<m;i++) c[i]=o[i];
  return c;
}

template <class T1, class T2, class T3>
struct triplet
{
	typedef T1 first_type;
	typedef T2 second_type;
	typedef T3 third_type;
	T1 first;
	T2 second;
	T3 third;
	triplet() : first(),second(),third()   {}
	triplet(const T1& x, const T2& y, const T3& z)  : first(x),second(y), third(z) {}
	template<class U, class V,class W> inline
	triplet(const triplet<U, V, W>& p) : first(p.first),second(p.second),third(p.third) {}
};

template <class T1, class T2,class T3>
inline triplet<T1, T2, T3> make_triplet(T1 x, T2 y, T3 z)
{
	return triplet<T1, T2, T3>(x, y, z);
}


template<class I,class R>
 class LinearComb : public E_F0mps { public: 
   typedef I TI;
   typedef R TR;
   typedef size_t size_type;
   typedef pair<I,R> K;
   typedef vector<K> array;
   typedef typename array::const_iterator const_iterator;
   typedef typename array::iterator iterator;
   array v;
   vector<size_type> where_in_stack_opt;
   vector<bool> mesh_indep_stack_opt;
   const E_F0_Optimize * optiexp0,*optiexpK;
   bool isoptimize;
   
   LinearComb(): v(),
		 where_in_stack_opt(),mesh_indep_stack_opt(),
		 optiexp0(0),optiexpK(0),isoptimize(false) {}
   
   LinearComb(const I& i,const R& r) :v(1),
				      where_in_stack_opt(),mesh_indep_stack_opt(),
				      optiexp0(),optiexpK(),isoptimize(false) 
   {
     v[0]=make_pair<I,R>(i,r);
   }
   
   LinearComb(const LinearComb &l) 
      :v(l.v),
        where_in_stack_opt(l.where_in_stack_opt),
        mesh_indep_stack_opt(l.mesh_indep_stack_opt),
        optiexp0(l.optiexp0),optiexpK(l.optiexpK),isoptimize(false){}  
        
   int nbtrue(bool *ok) const 
     { 
        int k=0; 
        for (size_t i=0;i<v.size();i++) 
          if (ok[i]) k++;
         return k;
     }
     
   //  copy just les parties ok   
   LinearComb(const LinearComb &l,bool * ok) 
     :v(l.nbtrue(ok) ),
      where_in_stack_opt(v.size()),
      mesh_indep_stack_opt(v.size()),
      optiexp0(l.optiexp0),optiexpK(l.optiexpK),
      isoptimize(l.isoptimize) 
      {
	// int k=0;
        const_iterator lll= l.v.begin();
        iterator ll=v.begin(); 
        for (int i=0,k=0;lll!= l.v.end();++lll,++i) 
          if (ok[i]) 
            {
            *ll++=*lll;
            where_in_stack_opt[k]=l.where_in_stack_opt[i];
            mesh_indep_stack_opt[k]=l.mesh_indep_stack_opt[i];
            ++k;            
            }
         
      }  
        
       
   void operator=(const LinearComb<I,R> &l) {
     v=l.v;
     where_in_stack_opt=l.where_in_stack_opt;
     mesh_indep_stack_opt=l.mesh_indep_stack_opt;
     optiexp0=l.optiexp0;
     optiexpK=l.optiexpK;
     isoptimize=l.isoptimize; 
     }
   const I * simple() const { if (v.size()==1) return & v.begin()->first;else return  0;}     
   void  add(const I& i,const R &r)  { 
     for (iterator k=v.begin();k!=v.end();k++)
       if (k->first == i) {k->second += r;return ;}
     v.push_back(make_pair<I,R>(i,r));     
   }
   
   size_type size() const { return v.size();}   
   
   const K & operator[](size_type i) const { return v[i];}
   
   void operator+=(const LinearComb & l) {
     for (const_iterator k=l.v.begin();k!=l.v.end();k++)
       { const K & kk(*k);
       add(kk.first,kk.second);}
   }
   
   void operator*=(const R & r) {
     for (iterator k=v.begin();k!=v.end();k++)
       {K & kk(*k);
       kk.second = kk.second*r;}
   }     
   
   void operator/=(const R & r) {
     for (iterator k=v.begin();k!=v.end();k++)
       {K & kk(*k);
       kk.second = kk.second/r; }
   }     
   
  AnyType operator()(Stack )  const {
    return SetAny<const LinearComb<I,R> * >(this);}
  operator aType () const { return atype<const LinearComb >();}         
    
    
  bool mappable(bool (*f)(const  R &)) const {
     for (const_iterator k=v.begin();k!=v.end();k++)
       if (!(*f)(k->second)) return false;
      return true;}
      
  void mapping(R (*f)(const R &)) {
     for (iterator k=v.begin();k!=v.end();k++)
       k->second=(*f)(k->second) ;}
  
  int MaxOp() const {
    int m=0;
     for (const_iterator k=v.begin();k!=v.end();k++)
       m= maxop(m,(k->first)) ;
     return m;}
    
 void DiffOp(KN_<bool> &d) const {
     assert(d.N() >= last_operatortype);
     d=false;
     for (const_iterator k=v.begin();k!=v.end();k++)
       SetOp(d,k->first);
       
   
 }  
  unsigned int DiffOp(int & lastop) const
 {
     unsigned int d=0;
    lastop=0;
    for (const_iterator k=v.begin();k!=v.end();k++)
       d |=GetDiffOp(k->first,lastop);
    assert(lastop<last_operatortype);
    lastop++;
    return d;
  }  
  
  ostream & dump(ostream &f) const {
    int n = size();
    for (int i=0; i<n; i++)
      {
        const K & ri=v[i];
        Expression ee= ri.second.LeftValue();
        f  << "\n\t\t"<< i << " " << ri.first << ": "  ;
       f << " :  type exp: " << typeid(*ee).name() << " "<<endl;
      }

    return f;
  } 
 LinearComb * Optimize(Block * b)  
  {
    const bool kdump=(verbosity/1000)%10==1;
    if (kdump)
    cout << "\n\n Optimize " << endl;
    LinearComb * r=new LinearComb(*this);
    LinearComb  &rr= *r;
    int n = rr.size();
    
    deque<pair<Expression,int> > ll;
    MapOfE_F0 m;
    rr.where_in_stack_opt.resize(n);
    rr.mesh_indep_stack_opt.resize(n);
    size_type top = b->OffSet(0), topbb=top; // FH. bofbof ??? 
    for (int i=0; i<n; i++)
      {
      const K & ri=rr.v[i];
       Expression ee= ri.second.LeftValue();
       if (kdump)
       cout << "Optimize :  type exp: " << typeid(*ee).name() << " "<<endl;
        rr.mesh_indep_stack_opt[i]=ee->MeshIndependent(); 
       rr.where_in_stack_opt[i]=ee->Optimize(ll, m, top);
       if (kdump)
       cout  << "\n\t\t"<< i << " " << ri.first << ": " << rr.where_in_stack_opt[i] << endl;
      }
      
    b->OffSet(top-topbb);
    //  
    int k=ll.size(),k0=0,k1=0;
    for (int i=0;i<k;i++)
        if (ll[i].first->MeshIndependent()) k0++;
    deque<pair<Expression,int> > l0(k0),l1(k-k0);
    k0=0,k1=0;
    for (int i=0;i<k;i++)
       if (ll[i].first->MeshIndependent()) 
         {
          if (kdump)
          cout << " mi " << ll[i].second << " " << *(ll[i].first) << endl;
          l0[k0++]=ll[i];
         }
        else 
         {
          if (kdump)
          cout << " md " << ll[i].second << " " << *(ll[i].first) << endl;
          l1[k1++]=ll[i];
         }
    if (k0)      
      rr.optiexp0 = new E_F0_Optimize(l0,m,0);  
    if (k1) 
      rr.optiexpK = new E_F0_Optimize(l1,m,0);
    rr.isoptimize=true;
    if (kdump) cout << "LinearCom Optimize k0(mi) = " << k0 << " k1 = " << k1 << "\n\n"<<endl;
    return r;
  }    
   
};



template<class I,class R>
LinearComb<I,R> operator+(const LinearComb<I,R> & a,const LinearComb<I,R> & b)
 {LinearComb<I,R> r(a);r+=b;return r;}

template<class I,class R>
LinearComb<I,R> operator*(const LinearComb<I,R> & a,const R & b)