array_tlp.hpp

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
 */
//#pragma dont_inline on
//#pragma inline_depth(1)

#include "config-wrapper.h"

#include <complex>
#include "AFunction.hpp"
#include <cstdarg>
#include <cstring>
#include "error.hpp"
#include "lex.hpp"

#include "RNM.hpp"

#include "Operator.hpp"
// for exec routine 
#include "rgraph.hpp"
#include "InitFunct.hpp"
#include <queue>
#include "array_resize.hpp"

template <class T>
struct affectation: binary_function<T, T, T>
{
	T& operator()(T& x, const T& y) const {return (x=y);}
};

template <class T>
struct affectation_add: binary_function<T, T, T>
{
	T& operator()(T& x, const T& y) const {return (x=+y);}
};

template <class T>
struct affectation_sub: binary_function<T, T, T>
{
	T& operator()(T& x, const T& y) const {return (x=-y);}
};



extern Map_type_of_map map_type_of_map ; //  to store te type 
extern Map_type_of_map map_pair_of_type ; //  to store te type 

extern basicForEachType *  typevarreal,  * typevarcomplex;  //  type of real and complex variable

extern int TheCurrentLine; // unset: by default
extern long mpisize,mpirank;

template<class T> inline T Max (const T &a,const T & b){return a > b ? a : b;}
template<class T> inline T Min (const T &a,const T & b){return a < b ? a : b;}
template<class T> inline T Abs (const T &a){return a <0 ? -a : a;}
template<class T> inline T Max (const T &a,const T & b,const T & c){return Max(Max(a,b),c);}
template<class T> inline T Min (const T &a,const T & b,const T & c){return Min(Min(a,b),c);}
template<class T> inline T Square (const T &a){return a*a;}


 
template<class K> 
struct Op2_dotproduct: public binary_function<Transpose<KN_<K> >,KN<K> *,K> { 
  static K f( Transpose<KN_<K> > const & a, KN<K> * const& b)  
   { return (conj(a.t),*b);} }; 

template<class K> 
struct Op2_dotproduct_: public binary_function<Transpose<KN_<K> >,KN_<K> ,K> { 
  static K f( Transpose<KN_<K> > const & a, KN_<K>  const& b)  
   { return (conj(a.t),b);} }; 
   
template<class A,class B>  A Build(B b) {  return A(b);}

template<class T>
void  HeapSort(T *c,long n,long o)
{ // trie un tableau c de n valeur avec un decalage de o.
   //  le tableau: c[i*o] , pour i = 0 a n-1  
    long l,j,r,i;
    T crit;
    c-=o; // on decale de o pour que le tableau commence a o
    if( n <= 1) return;
    l = (n/2 + 1)*o;
    r = n*o;
    while (1) { // label 2
	if(l <= o ) { // label 20
	    crit = c[r];
	    c[r] = c[o];
	    r-=o;
	    if ( r == o ) { c[o]=crit; return;}
	} else  crit = c[l-=o]; 
	j=l;
	while (1) {// label 4
	    i=j;
	    j=2*j;
	    if  (j>r) {c[i]=crit;break;} // L8 -> G2
	    if ((j<r) && (c[j] < c[j+o])) j+=o; // L5
	    if (crit < c[j]) c[i]=c[j]; // L6+1 G4
	    else {c[i]=crit;break;} //L8 -> G2
	}
    }
}

template<class R,class A> A  SortKn(const A  & ca){ 
    A a(ca);
    HeapSort<R>(&a[0],a.n,a.step);
    return a;}

template<class R> KN<R> *  SortpKn( KN<R> * const & pa){ 
    KN<R> &a(*pa);
    HeapSort<R>(&a[0],a.n,a.step);
    return pa;}

template<class R>
class QuantileKN:  public KN_<R> { public:
    QuantileKN(const KN_<R> &a): KN_<R>(a) {}
    QuantileKN(KN<R>  * p): KN_<R>(*p) {}
    operator R *() const {return this->KN_<R>::operator R *() ;}    
};


template<class R> R   Quantile(QuantileKN<R>  const & a,const double & q){ 
    KN<R> b(a); 
    HeapSort<R>(b,b.n,b.step);
    long m=lrint(b.n*q);
    if( m >= b.n) m=b.n-1;
    if( m < 0) m=0;   
    R qq=b[m];
   // cout <<  "Quantile: m = " << m << " " << b <<endl;
    return qq;}


  
inline void MyAssert(int i,char * ex,char * file,long line)
{if (i) {
    cout << "CompileError assertion :  " << ex << " in file " << file << "  line = " << line << endl; 
     CompileError();}
 }


template<class K>
inline   K * get_element( MyMap<String,K> *  const  &  a,string*  const   & b)
 { K * ret=  &((*a)[*b]); // correction FH feb 2004
  //  cout << "get_element " << *b << " : " << ret << " = "<< * ret << endl;
   // delete b;  modif mars 2006 auto del ptr
    return ret;}
    
template<>
inline   string ** get_element<string*>( MyMap<String,string*> *  const  &  a,string*  const   & b)
 { string** ret=  &((*a)[*b]); // correction FH feb 2004
    if( *ret ==0) *ret = new string(""); //  string vide ???
     // cout << "get_element " << *b << " : " << ret << " = "<< * ret << endl;
    // delete b;  modif mars 2006 auto del ptr
    return ret;}

inline   string ** get_elements( MyMap<String,String> *  const  &  a,string*  const   & b)
 { String* Sret=  &((*a)[*b]); // correction FH feb 2004
   //  delete b;  modif mars 2006 auto del ptr
    return Sret->getap();}

template<class RR,class A,class B>  
RR * get_element_(const A & a,const B & b){ 
  if( b<0 || a.N() <= b) 
   { cerr << " Out of bound  0 <=" << b << " < "  << a.N() << " array type = " << typeid(A).name() << endl;
     ExecError("Out of bound in operator []");}
    return  &((a)[b]);}
    

template<class RR,class A,class B>  
RR * get_elementp_(const A & a,const B & b){ 
  if( b<0 || a->N() <= b) 
   { cerr << " Out of bound  0 <=" << b << " < "  << a->N() << " array type = " << typeid(A).name() << endl;
     ExecError("Out of bound in operator []");}
    return  &((*a)[b]);}

template<class K>  
KN_<K> fSubArray(const KN_<K> & a,const SubArray & b)
 { return a(b);}
template<class K>  
KN_<K> fSubArrayp( KN<K>  * const & a,const SubArray & b)
 { return (*a)(b);}
 
template<class A>  
A fSubArrayc(const A & a,const char & b)
 { return a;}
 
template<class RR,class A,class B,class C>  
RR * get_elementp2_(const A & a,const B & b,const C & c){ 
  if( b<0 || a->N() <= b || c<0 || a->M() <= c  ) 
   { cerr << " Out of bound  0 <=" << b << " < "  << a->N() << " " << c << " < "  << a->M() 
           << " array type = " << typeid(A).name() << endl;
     ExecError("Out of bound in operator (,)");}
    return  &((*a)(b,c));}

template<class RR,class A,class B,class C>  
RR get_element_is(const A &  a,const B & b,const C & c){ 
 //  cout << b << " .... " << ((*a)(SubArray(1,b),c)) << endl;;
    return  ((*a)(b,'.')(c));}

template<class RR,class A,class B,class C>  
RR get_element_si(const A &  a,const B & b,const C & c){ 
 //  cout << c << " .... " << ((*a)(b,SubArray(1,c) )) << endl;;
     return  ((*a)('.',c)(b));}
     
template<class RR,class A,class B,class C>  
RR get_element_lineorcol(const A &  a,const B & b,const C & c){ 
 //  cout << b << " .... " << ((*a)(SubArray(1,b),c)) << endl;;
    return  ((*a)(b,c));}

    

template<class RR,bool isinit>
class  InitArrayfromArray : public OneOperator { 
public:
    typedef KN<RR> * A;
    typedef KN<RR> * R;
    typedef E_Array B;
    
    class CODE : public  E_F0 { public:
       Expression a0;
       int N;
       Expression * tab;
    int * what;//  0  RR, 1 KN<RR>, 
       const  bool mi;

    CODE(Expression a,const E_Array & tt)  
      : a0(a),N(tt.size()),
	tab(new Expression [N]),
	what(new int[N])  ,
	mi(tt.MeshIndependent())

      {
        assert(&tt);
	//      int err=0;
        for (int i=0;i<N;i++)
	if(atype<RR>()->CastingFrom(tt[i].right() ) ) 
	  {
          tab[i]=atype<RR>()->CastTo(tt[i]);
	    what[i]=0;
	  }
	else if(atype<KN_<RR> >()->CastingFrom(tt[i].right() ) ) 
	  {
	    tab[i]=atype<KN_<RR> >()->CastTo(tt[i].RightExp());
	    what[i]=1;
	  }      
	else 
	  CompileError(" we are waiting for scalar or vector of scalar");
    }
    
    AnyType operator()(Stack stack)  const 
    {
      A  a=GetAny<A>((*a0)(stack));
      KN<AnyType> v(N);
      KN<int>  nn(N+1);
      for (int i=0;i<N;i++)
        v[i]= (*(tab[i]))(stack);
      
      int n=0;
      for (int i=0;i<N;i++)
	{
	  if (what[i]==0) nn[i]=1;
	  else if (what[i]==1) nn[i]=GetAny<KN_<RR> >(v[i]).size();
          n += nn[i];
	}
      if (isinit) 
        a->init(n);
      else
	a->resize(n);
      
      for (int i=0,j=0 ;i<N; j += nn[i++])
	
        if (what[i]==0)
          (*a)[j]= GetAny<RR>(v[i]);
        else if (what[i]==1) 
          (*a)(SubArray(nn[i],j)) = GetAny<KN_<RR> >(v[i]);
      return SetAny<R>(a);
    } 
    bool MeshIndependent() const     {return  mi;} // 
    ~CODE() { delete [] tab; delete[] what;}
    operator aType () const { return atype<R>();}    
  }; // end sub class CODE
  
  
    public: 
    E_F0 * code(const basicAC_F0 & args) const 
     { return  new CODE(t[0]->CastTo(args[0]),*dynamic_cast<const E_Array*>( t[1]->CastTo(args[1]).LeftValue()));} 
    InitArrayfromArray():   OneOperator(atype<R>(),atype<A>(),atype<B>())  {}
  
};

template<class RR,bool isinit>
class  InitMatfromAArray : public OneOperator { 
public:
    typedef KNM<RR> * A;
    typedef KNM<RR> * R;
    typedef E_Array B;
    
    class CODE : public  E_F0 { public:
       Expression a0;
       int N;
       int M;
       Expression ** tab;
       const  bool mi;

    CODE(Expression a,const E_Array & tt)  
      : a0(a),N(tt.size()),M(0),
	tab(new Expression* [N]),
	mi(tt.MeshIndependent())

      {
        assert(&tt);
	//       int err=0;
        for (int i=0;i<N;i++)
         {
          const E_Array  *li =  dynamic_cast<const E_Array *>(tt[i].LeftValue());
          if (li)
	  {
	     const E_Array & lli = *li;
	     // -- check ---
	     if( i == 0) {