1 // ********** DO NOT REMOVE THIS BANNER **********
2 // ORIG-DATE: Jan 2008
3 // -*- Mode : c++ -*-
4 //
5 // SUMMARY : Generic Fiinite Element header 1d, 2d, 3d
6 // USAGE : LGPL
7 // ORG : LJLL Universite Pierre et Marie Curi, Paris, FRANCE
8 // AUTHOR : Frederic Hecht
9 // E-MAIL : frederic.hecht@ann.jussieu.fr
10 //
11
12 /*
13
14 This file is part of Freefem++
15
16 Freefem++ is free software; you can redistribute it and/or modify
17 it under the terms of the GNU Lesser General Public License as published by
18 the Free Software Foundation; either version 2.1 of the License, or
19 (at your option) any later version.
20
21 Freefem++ is distributed in the hope that it will be useful,
22 but WITHOUT ANY WARRANTY; without even the implied warranty of
23 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 GNU Lesser General Public License for more details.
25
26 You should have received a copy of the GNU Lesser General Public License
27 along with Freefem++; if not, write to the Free Software
28 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
29
30 Thank to the ARN () FF2A3 grant
31 ref:ANR-07-CIS7-002-01
32 */
33
34 #ifndef FESpacen_HPP_
35 #define FESpacen_HPP_
36 /*
37 * FEspacen.hpp
38 * EF23n
39 *
40 * Created by Fr�d�ric Hecht on 04/12/07.
41 * Copyright 2007 Universite Pierre et marie Curie All rights reserved.
42 *
43 */
44 #include <cassert>
45 #include <cmath>
46 #include <cstdlib>
47 #include <fstream>
48 #include <iostream>
49 #include <complex>
50 #include <map>
51 using namespace std;
52 #include "error.hpp"
53 #include "ufunction.hpp"
54 #include "MeshLn.hpp"
55 #include "MeshSn.hpp"
56 #include "Mesh3dn.hpp"
57 #include "Mesh2dn.hpp"
58 #include "Mesh1dn.hpp"
59
60 #include "RNM.hpp"
61
62
63 #include "QuadratureFormular.hpp"
64
65 namespace Fem2D {
66
67 template<class Mesh> class GFESpace;
68 template<class Mesh> class GFElement;
69 template<class Mesh> class GbaseFElement;
70 template<class Mesh> class GTypeOfFE;
71
72 // numbering of derivative
73 enum operatortype {
74 op_id=0,
75 op_dx=1,op_dy=2,
76 op_dxx=3,op_dyy=4,
77 op_dyx=5,op_dxy=5,
78 op_dz=6,
79 op_dzz=7,
80 op_dzx=8,op_dxz=8,
81 op_dzy=9,op_dyz=9
82 };
83
84 typedef unsigned int What_d;
85
86 const unsigned int Fop_id= 1<< op_id;
87
88 const unsigned int Fop_dx= 1<< op_dx;
89 const unsigned int Fop_dy= 1<< op_dy;
90 const unsigned int Fop_dz= 1<< op_dz;
91
92 const unsigned int Fop_dxx= 1<< op_dxx;
93 const unsigned int Fop_dxy= 1<< op_dxy;
94 const unsigned int Fop_dxz= 1<< op_dxz;
95
96 const unsigned int Fop_dyx= 1<< op_dyx;
97 const unsigned int Fop_dyy= 1<< op_dyy;
98 const unsigned int Fop_dyz= 1<< op_dyz;
99
100 const unsigned int Fop_dzx= 1<< op_dzx;
101 const unsigned int Fop_dzy= 1<< op_dzy;
102 const unsigned int Fop_dzz= 1<< op_dzz;
103
104 const unsigned int Fop_D0 = Fop_id;
105 const unsigned int Fop_D1 = Fop_dx | Fop_dy | Fop_dz;
106 const unsigned int Fop_D2 = Fop_dxx | Fop_dyy | Fop_dzz | Fop_dxy | Fop_dxz | Fop_dyz;
107 const unsigned int Fop_Dall = Fop_D0| Fop_D1| Fop_D2;
108
Fwhatd(const operatortype op)109 inline What_d Fwhatd(const operatortype op) { return 1<< op;}
110
111
112 const int last_operatortype=10;
113 const bool operatortypeValue[last_operatortype]= {true,false,false,false,false,false,false,false,false,false} ;
114
115
initwhatd(bool * whatd,int k)116 inline void initwhatd(bool *whatd,int k)
117 {
118 for (int i=0;i<last_operatortype;i++)
119 whatd[i]=false;
120 whatd[k]=true;
121 }
122
123 typedef double R;
124 typedef KN_<R> RN_;
125 typedef KN<R> RN;
126 typedef KNM_<R> RNM_;
127 typedef KNMK_<R> RNMK_;
128 typedef KNMK<R> RNMK;
129
130 template<class Mesh> class GFElement;
131 template<class Mesh> class GFESpace;
132 template<class Mesh>class GTypeOfFE ;
133
134 class dataTypeOfFE {
135
136
137 private:
138 const int * data;
139 const int * dataalloc;
140
141 public:
142 const int ndfonVertex;
143 const int ndfonEdge;
144 const int ndfonFace;
145 const int ndfonVolume;
146 const int NbDoF;
147 const int NbNode;
148 int N,nb_sub_fem;
149 const int nbsubdivision; // nb of subdivision for plot
150 const bool discontinue;
151
152 int const * const DFOnWhat;
153 int const * const DFOfNode; // nu du df on Node
154 int const * const NodeOfDF; // nu du node du df
155 int const * const fromFE; // the df come from df of FE
156 int const * const fromDF; // the df come from with FE
157 int const * const fromASubFE; // avril 2006 for CL
158 int const * const fromASubDF; // avril 2006 for CL
159 int const * const dim_which_sub_fem; // from atomic sub FE for CL
ndfOn() const160 const int * ndfOn() const { return & ndfonVertex;}
161
162 dataTypeOfFE(const int *nnitemdim,const int dfon[4],int NN,int nbsubdivisionn,int nb_sub_femm=1,bool discon=true);
163 // pour evite un template
164 // nitemdim : nbitem : si d==2 3,3,1,0 , si d=3: 4,6,4,1 , si d==1 = 2,1,0,0
165 // dfon : nombre de df par item
166 // NN
167 dataTypeOfFE(const int nitemdim[4],const KN< dataTypeOfFE const *> & tef);
168
~dataTypeOfFE()169 virtual ~dataTypeOfFE(){ if(dataalloc) delete [] dataalloc;}
170 };
171
172 template<class RdHat>
173 class InterpolationMatrix {
174 public:
175 const int N,np,ncoef;
176 bool invariant;
177 int k;
178 KN<RdHat> P;
179 KN<R> coef;
180 KN<int> comp;
181 KN<int> p;
182 KN<int> dofe;
183
184 template<class Mesh>
185 InterpolationMatrix(const GFESpace<Mesh> &Vh);
186
187 template<class Mesh>
188 InterpolationMatrix(const GTypeOfFE<Mesh> & tef);
189
190 template<class Mesh>
191 void set(const GFElement<Mesh> & FK);
192
193
194 private: // copie interdit ...
195 InterpolationMatrix(const InterpolationMatrix &);
196 void operator=(const InterpolationMatrix &);
197 };
198
199 template <class RdHat>
operator <<(ostream & f,const InterpolationMatrix<RdHat> & M)200 ostream & operator<<(ostream& f,const InterpolationMatrix<RdHat> &M)
201 { f<< M.N << " "<< M.k << " "<< M.np << " "<< M.ncoef << endl;
202 f<< " = " << M.P ;
203 f << "coef=" <<M.coef ;
204 f << "comp " << M.comp;
205 f << "p = " << M.p;
206 f << "dofe = " << M.dofe;
207 return f;
208 }
209
210
211 template<class Mesh>
212 class GTypeOfFE : public dataTypeOfFE
213 {
214 public:
215 typedef typename Mesh::Element Element;
216 typedef Element E;
217 typedef typename Element::RdHat RdHat;
218 typedef typename Element::Rd Rd;
219 typedef GFElement<Mesh> FElement;
220
221 // generic data build interpolation
222 bool invariantinterpolationMatrix;
223 int NbPtforInterpolation;
224 int NbcoefforInterpolation;
225 KN<RdHat> PtInterpolation;
226 KN<int> pInterpolation,cInterpolation,dofInterpolation;
227 KN<R> coefInterpolation;
228
229 // soit
230 // $(U_pj)_{{j=0,N-1}; {p=0,nbpoint_Pi_h-1}}$ les valeurs de U au point Phat[i];
231 // p est le numero du point d'integration
232 // j est la composante
233 // l'interpole est defini par
234 // Pi_h u = \sum_k u_k w^k , et u_i = \sum_pj alpha_ipj U_pj
235 // la matrice de alpha_ipj est tres creuse
236
237
238
239 KN<GTypeOfFE<Mesh> *> Sub_ToFE;
240 KN<int> begin_dfcomp, end_dfcomp;
241 KN<int> first_comp,last_comp; // for each SubFE
242
init(InterpolationMatrix<RdHat> & M,FElement * pK,int odf,int ocomp,int * pp) const243 virtual void init(InterpolationMatrix<RdHat> & M,FElement * pK,int odf,int ocomp,int *pp) const
244 {
245 assert(pK==0);
246 assert(M.np==NbPtforInterpolation);
247 assert(M.ncoef==NbcoefforInterpolation);
248 M.P=PtInterpolation;
249 M.coef=coefInterpolation;
250 M.comp=cInterpolation;
251 M.p=pInterpolation;
252 M.dofe=dofInterpolation;
253 }
254
255
256 // the full constructor ..
GTypeOfFE(const int dfon[4],const int NN,int nsub,int kPi,int npPi,bool invar,bool discon)257 GTypeOfFE(const int dfon[4],const int NN,int nsub,int kPi,int npPi,bool invar,bool discon)
258 :
259 dataTypeOfFE(Element::nitemdim,dfon, NN, nsub,1,discon),
260 invariantinterpolationMatrix(invar),
261 NbPtforInterpolation(npPi),
262 NbcoefforInterpolation(kPi),
263 PtInterpolation(NbPtforInterpolation),
264 pInterpolation(NbcoefforInterpolation),
265 cInterpolation(NbcoefforInterpolation),
266 dofInterpolation(NbcoefforInterpolation),
267 coefInterpolation(NbcoefforInterpolation),
268
269 Sub_ToFE(this->nb_sub_fem),
270 begin_dfcomp(N,0),
271 end_dfcomp(N,this->NbDoF),
272 first_comp(this->nb_sub_fem,0),
273 last_comp(this->nb_sub_fem,N)
274
275 {
276 Sub_ToFE=this;
277 assert(this->dim_which_sub_fem[this->N-1]>=0 && this->dim_which_sub_fem[this->N-1]< this->nb_sub_fem);
278 }
279
280 // simple constructeur of lagrange type Finite element 1 point par Node for interpolation
GTypeOfFE(const int dfon[4],const int NN,int nsub,bool invar,bool discon)281 GTypeOfFE(const int dfon[4],const int NN,int nsub,bool invar,bool discon)
282 :
283 dataTypeOfFE(Element::nitemdim,dfon, NN, nsub,1,discon),
284
285 invariantinterpolationMatrix(invar),
286 NbPtforInterpolation(this->NbDoF),
287 NbcoefforInterpolation(this->NbDoF),
288 PtInterpolation(NbPtforInterpolation),
289 pInterpolation(NbcoefforInterpolation),
290 cInterpolation(NbcoefforInterpolation),
291 dofInterpolation(NbcoefforInterpolation),
292 coefInterpolation(NbcoefforInterpolation),
293
294 Sub_ToFE(this->nb_sub_fem),
295 begin_dfcomp(N,0),
296 end_dfcomp(N,this->NbDoF),
297 first_comp(this->nb_sub_fem,0),
298 last_comp(this->nb_sub_fem,N)
299 {
300 Sub_ToFE=this;
301 assert(this->dim_which_sub_fem[this->N-1]>=0 && this->dim_which_sub_fem[this->N-1]< this->nb_sub_fem);
302 }
303
304 // simple constructeur pour sum direct d'espace
GTypeOfFE(const KN<GTypeOfFE const * > & tef)305 GTypeOfFE(const KN<GTypeOfFE const *> & tef)
306 :
307 dataTypeOfFE(Element::nitemdim,tef),
308
309 invariantinterpolationMatrix(0),
310 NbPtforInterpolation(0),
311 NbcoefforInterpolation(ncoeftef(tef)),
312 PtInterpolation(NbPtforInterpolation),
313 pInterpolation(NbcoefforInterpolation),
314 cInterpolation(NbcoefforInterpolation),
315 dofInterpolation(NbcoefforInterpolation),
316 coefInterpolation(NbcoefforInterpolation),
317
318
319 Sub_ToFE(this->nb_sub_fem),
320 begin_dfcomp(this->N,0),
321 end_dfcomp(this->N,this->NbDoF),
322 first_comp(this->nb_sub_fem,0),
323 last_comp(this->nb_sub_fem,N)
324 {
325 Sub_ToFE=this;
326 assert(this->dim_which_sub_fem[this->N-1]>=0 && this->dim_which_sub_fem[this->N-1]< this->nb_sub_fem);
327 }
328
329
330 virtual R operator()(const GFElement<Mesh> & K,const RdHat & PHat,const KN_<R> & u,int componante,int op) const ;
331 virtual void FB(const What_d whatd,const Mesh & Th,const Element & K,const RdHat &PHat, KNMK_<R> & val) const =0;
set(const Mesh & Th,const Element & K,InterpolationMatrix<RdHat> & M,int ocoef,int odf,int * nump) const332 virtual void set(const Mesh & Th,const Element & K,InterpolationMatrix<RdHat> & M,int ocoef,int odf,int *nump ) const {}; // no change by deflaut
333 // ocoef is the offset of the coef in M
334 // odf is the offset in the df
335 // nump if exist give the numbering of p . 0=> no change
336
tefN(const KN<GTypeOfFE const * > & tef)337 static KN<int> tefN(const KN<GTypeOfFE const *> & tef) {
338 int n=tef.N();
339 KN<int> ntef(n);
340 for(int i=0;i<n;++i) ntef[i]=tef[i]->N;
341 return ntef;}
342
ncoeftef(const KN<GTypeOfFE const * > & tef)343 static int ncoeftef(const KN<GTypeOfFE const *> & tef) {
344 int k=0,n=tef.N();
345 for(int i=0;i<n;++i)
346 k+= tef[i]->NbcoefforInterpolation;
347 return k;}
348
349
350 private:
Count(const int * data,int n,int which)351 static int Count(const int *data,int n,int which)
352 {
353 int kk=0;
354 for (int i=0;i<n;++i)
355 if (which == data[i]) ++kk;
356 return kk;
357 }
358
LastNode(const int * data,int n)359 static int LastNode(const int *data,int n)
360 {
361 int kk=0,i0=2*n;
362 for(int i=0;i<n;i++)
363 kk=Max(kk,data[i+i0]);
364 return kk;}
365
NbNodebyWhat(const int * data,int n,int on)366 static int NbNodebyWhat(const int *data,int n,int on)
367 {
368 const int nmax=100;
369 int t[nmax];
370 for (int i=0;i<nmax;i++)// correct FH 2Oct 2015 nmax
371 t[i]=0;
372 int k=0,i0=2*n;
373 for(int i=0;i<n;i++)
374 if ( data[i]==on)
375 { int node= data[i+i0];
376 // cout << " node " << node << endl;
377 ffassert(node < nmax);
378 if( ! t[node] )
379 t[node] = 1,++k;
380 }
381
382 // cout << " on " << on << " k = " << k << endl;
383 return k;
384 }
385
386 } ;
387
388 template<class mesh>
389 struct DataFE
390 {
391 static GTypeOfFE<mesh> & P0;
392 static GTypeOfFE<mesh> & P1;
393 static GTypeOfFE<mesh> & P2;
394 static GTypeOfFE<mesh> & RT0;
395 };
396
397
398
399 template<class MMesh>
400 class GbaseFElement
401 {
402 public:
403 typedef MMesh Mesh;
404 typedef GFESpace<Mesh> FESpace;
405 typedef typename Mesh::Element Element;
406 typedef Element E;
407 typedef typename E::Rd Rd;
408 typedef typename E::RdHat RdHat;
409 typedef Fem2D::GQuadratureFormular<RdHat> QF;
410
411 const GFESpace<Mesh> &Vh;
412 const Element &T;
413 const GTypeOfFE<Mesh> * tfe;
414 const int N;
415 const int number;
416 GbaseFElement(const GFESpace<Mesh> &aVh, int k) ;
417 GbaseFElement(const GbaseFElement & K, const GTypeOfFE<Mesh> & atfe) ;
EdgeOrientation(int i) const418 R EdgeOrientation(int i) const {return T.EdgeOrientation(i);}
419
420 };
421
422 template<class Mesh>
423 class GFElement : public GbaseFElement<Mesh>
424 {
425 public:
426
427 typedef typename Mesh::Element Element;
428 typedef Element E;
429 typedef typename E::Rd Rd;
430 typedef typename E::RdHat RdHat;
431 typedef Fem2D::GQuadratureFormular<RdHat> QF;
432
433
434 friend class GFESpace<Mesh>;
435 const int *p;
436 const int nb;
437
438 GFElement(const GFESpace<Mesh> * VVh,int k) ;
439
NbOfNodes() const440 int NbOfNodes()const {return nb;}
441 int operator[](int i) const ;//{ return p ? p+i : ((&T[i])-Vh.Th.vertices);} N\u00c9 du noeud
442 int NbDoF(int i) const ; // number of DF
443 int operator()(int i,int df) const ;// { N\u00c9 du DoF du noeud i de df local df
operator ()(int df) const444 int operator()(int df) const { return operator()(NodeOfDF(df),DFOfNode(df));}
445 // void Draw(const KN_<R> & U, const KN_<R> & VIso,int j=0) const ;
446 //void Drawfill(const KN_<R> & U, const KN_<R> & VIso,int j=0,double rapz=1) const ;
447 //void Draw(const RN_& U,const RN_& V, const KN_<R> & Viso,R coef,int i0,int i1) const;
448 //Rd MinMax(const RN_& U,const RN_& V,int i0,int i1) const ;
449 //Rd MinMax(const RN_& U,int i0) const ;
450 void BF(const RdHat & PHat,RNMK_ & val) const;// { tfe->FB(Vh.Th,T,P,val);}
451 void BF(const What_d whatd, const RdHat & PHat,RNMK_ & val) const;// { tfe->FB(Vh.Th,T,P,val);}
set(InterpolationMatrix<RdHat> & M) const452 void set(InterpolationMatrix<RdHat> &M) const {this->tfe->set(this->Vh.Th,this->T,M,0,0,0);}
453 // add april 08 begin end number for df of the componante ic
dfcbegin(int ic) const454 int dfcbegin(int ic) const { return this->tfe->begin_dfcomp[ic];}
dfcend(int ic) const455 int dfcend(int ic) const { return this->tfe->end_dfcomp[ic];}
456 // the fist and last composant of a sub finite element
457 // int firstcomp(int isfe) const {return this->tfe->first_comp[isfe];}
458 // int lastcomp(int isfe) const {return this->tfe->last_comp[isfe];}
subFE(int df) const459 int subFE(int df) const {return this->tfe->fromASubFE[df] ;}
460
461 template<class RR>
Pi_h(KNM_<RR> vpt,KN_<RR> & vdf,InterpolationMatrix<RdHat> & M) const462 KN_<RR> & Pi_h(KNM_<RR> vpt,KN_<RR> & vdf,InterpolationMatrix<RdHat> &M) const
463 {
464 // compute the interpolation
465 // in : vpt value of componant j at point p : vpt(p,j)
466 // out: vdf value du the degre of freedom
467 vdf=RR();
468
469 if(M.k != this->number)
470 M.set( (const GFElement &) *this);
471
472 for (int k=0;k<M.ncoef;++k)
473 vdf[M.dofe[k]] += M.coef[k]*vpt(M.p[k],M.comp[k]);
474
475 return vdf;
476 }
477
478
NbDoF() const479 int NbDoF() const {return this->tfe->NbDoF;}
DFOnWhat(int i) const480 int DFOnWhat(int i) const { return this->tfe->DFOnWhat[i];}
FromDF(int i) const481 int FromDF(int i) const { return this->tfe->fromDF[i];}
FromFE(int i) const482 int FromFE(int i) const { return this->tfe->fromFE[i];}
483
484 // df is the df in element
NodeOfDF(int df) const485 int NodeOfDF(int df) const { return this->tfe->NodeOfDF[df];} // a node
FromASubFE(int i) const486 int FromASubFE(int i) const { return this->tfe->fromASubFE[i];}
FromASubDF(int i) const487 int FromASubDF(int i) const { return this->tfe->fromASubDF[i];}
DFOfNode(int df) const488 int DFOfNode(int df) const { return this->tfe->DFOfNode[df];} // the df number on the node
489
490 R operator()(const RdHat & PHat,const KN_<R> & u,int i,int op) const ;
491 complex<R> operator()(const RdHat & PHat,const KN_<complex<R> > & u,int i,int op) const ;
492
493 // GFElementGlobalToLocal operator()(const KN_<R> & u ) const { return GFElementGlobalToLocal(*this,u);}
494 private:
nbsubdivision() const495 int nbsubdivision() const { return this->tfe->nbsubdivision;} // for draw
496 };
497
498
499 template<class MMesh>
500 class BuildTFE { protected:
501 GTypeOfFE<MMesh> const * const tfe;
502 };
503
504 template<class MMesh>
505 struct GFESpacePtrTFE {
506 GTypeOfFE<MMesh> const * const ptrTFE;
GFESpacePtrTFEFem2D::GFESpacePtrTFE507 GFESpacePtrTFE(GTypeOfFE<MMesh> const * const pptrTFE=0) : ptrTFE(pptrTFE) {}
~GFESpacePtrTFEFem2D::GFESpacePtrTFE508 virtual ~GFESpacePtrTFE() { if(ptrTFE) delete ptrTFE;}
509 };
510
511 template<class MMesh>
512 class GFESpace : public RefCounter,protected GFESpacePtrTFE<MMesh>, public DataFENodeDF, public UniqueffId {
513 public:
514 typedef MMesh Mesh;
515 typedef GFElement<Mesh> FElement;
516 typedef typename Mesh::Element Element;
517 typedef typename Mesh::BorderElement BorderElement;
518 typedef typename Mesh::Rd Rd;
519 typedef GTypeOfFE<Mesh> TypeOfFE;
520 typedef GQuadratureFormular<typename Element::RdHat> QFElement;
521 typedef GQuadratureFormular<typename BorderElement::RdHat> QFBorderElement;
522
523 const Mesh &Th;
524
525 KN<const GTypeOfFE<Mesh> *> TFE;
526 private:
527
528 public:
529 CountPointer<const Mesh> cmesh;
530 const int N; // dim espace d'arrive
531 const int Nproduit; // 1 if non constant Max number df par node. else Max number df par node..
532 const int nb_sub_fem; // nb de sous elements finis tensorise (independe au niveau des composantes)
533 int const* const dim_which_sub_fem;// donne les dependant des composantes liee a un meme sous element fini
534 const int maxNbPtforInterpolation;
535 const int maxNbcoefforInterpolation;
536
537 // exemple si N=5,
538 // dim_which_sub_fem[0]=0;
539 // dim_which_sub_fem[1] =1;
540 // dim_which_sub_fem[2]= 2
541 // dim_which_sub_fem[3] =2
542 // dim_which_sub_fem[4] =3
543 // =>
544 // le sous elements fini 0 est lie a la composante 0
545 // le sous elements fini 1 est lie a la composante 1
546 // le sous elements fini 2 est lie aux composantes 2,3
547 // le sous elements fini 3 est lie a la composante 4
548 // donc pour les CL. les composante 2 et 3 sont lie car elle sont utiliser pour definir un
549 // meme degre de libert\u00e9.
550
551 // par defaut P1
552
GFESpace(const Mesh & TTh,const GTypeOfFE<Mesh> & tfe=DataFE<Mesh>::P1,int nbequibe=0,int * equibe=0)553 GFESpace(const Mesh & TTh,const GTypeOfFE<Mesh> & tfe=DataFE<Mesh>::P1,int nbequibe=0,int *equibe=0)
554 :
555 GFESpacePtrTFE<MMesh>(0),
556 DataFENodeDF(TTh.BuildDFNumbering(tfe.ndfonVertex,tfe.ndfonEdge,tfe.ndfonFace,tfe.ndfonVolume,nbequibe,equibe)),
557 Th(TTh),
558 TFE(1,0,&tfe),
559 cmesh(TTh),
560 N(TFE[0]->N),
561 Nproduit(FirstDfOfNodeData ? 1 :MaxNbDFPerNode ),
562 nb_sub_fem(TFE[0]->nb_sub_fem),
563 dim_which_sub_fem(TFE[0]->dim_which_sub_fem),
564 maxNbPtforInterpolation(TFE[0]->NbPtforInterpolation),
565 maxNbcoefforInterpolation(TFE[0]->NbcoefforInterpolation)
566
567 {
568 if(!lockOrientation) {
569 cerr << " Error, lockOrientation must be true to build fespace ; must check orientation element for mesh" ;
570 assert(lockOrientation);
571 }
572 if(verbosity) cout << " -- FESpace: Nb of Nodes " << NbOfNodes
573 << " Nb of DoF " << NbOfDF <<endl;
574 }
575
576 GFESpace(const GFESpace & Vh,int kk,int nbequibe=0,int *equibe=0);
577 GFESpace(const GFESpace ** Vh,int kk,int nbequibe=0,int *equibe=0);
578
FirstDFOfNode(int i) const579 int FirstDFOfNode(int i) const {return FirstDfOfNodeData ? FirstDfOfNodeData[i] : i*Nproduit;}
LastDFOfNode(int i) const580 int LastDFOfNode(int i) const {return FirstDfOfNodeData ? FirstDfOfNodeData[i+1] : (i+1)*Nproduit;}
NbDFOfNode(int i) const581 int NbDFOfNode(int i) const {return FirstDfOfNodeData ? FirstDfOfNodeData[i+1]-FirstDfOfNodeData[i] : Nproduit;}
MaximalNbOfNodes() const582 int MaximalNbOfNodes() const {return MaxNbNodePerElement;};
MaximalNbOfDF() const583 int MaximalNbOfDF() const {return MaxNbDFPerElement;};
PtrFirstNodeOfElement(int k) const584 const int * PtrFirstNodeOfElement(int k) const {
585 return NodesOfElement
586 ? NodesOfElement + (FirstNodeOfElement ? FirstNodeOfElement[k] : k*MaxNbNodePerElement)
587 : 0;}
588
SizeToStoreAllNodeofElement() const589 int SizeToStoreAllNodeofElement() const {
590 return FirstNodeOfElement
591 ? FirstNodeOfElement[NbOfElements]
592 : MaxNbNodePerElement*NbOfElements;}
593
NbOfNodesInElement(int k) const594 int NbOfNodesInElement(int k) const {
595 return FirstNodeOfElement
596 ? FirstNodeOfElement[k+1] - FirstNodeOfElement[k]
597 : MaxNbNodePerElement ;}
esize() const598 int esize() const { return MaxNbDFPerElement*N*last_operatortype;} // size to store all value of B. function
599
600
~GFESpace()601 ~GFESpace()
602 {
603 }
604 // GFESpace(Mesh & TTh,int NbDfOnSommet,int NbDfOnEdge,int NbDfOnElement,int NN=1);
605 int renum();
606
operator [](int k) const607 FElement operator[](int k) const { return FElement(this,k);}
operator [](const Element & K) const608 FElement operator[](const Element & K) const { return FElement(this,Th(K));}
operator ()(int k) const609 int operator()(int k)const {return NbOfNodesInElement(k);}
operator ()(int k,int i) const610 int operator()(int k,int i) const { // the node i of element k
611 return NodesOfElement ? *(PtrFirstNodeOfElement(k) + i) : Th(k,i) ;}
612
613 /*
614 void Draw(const KN_<R>& U,const KN_<R>& Viso,int j=0,float *colors=0,int nbcolors=0,bool hsv=true,bool drawborder=true) const ; // Draw iso line
615 void Drawfill(const KN_<R>& U,const KN_<R>& Viso,int j=0,double rapz=1,float *colors=0,int nbcolors=0,bool hsv=true,bool drawborder=true) const ; // Draw iso line
616
617 void Draw(const KN_<R>& U,const KN_<R> & Viso, R coef,int j0=0,int j1=1,float *colors=0,int nbcolors=0,bool hsv=true,bool drawborder=true) const ; // Arrow
618 void Draw(const KN_<R>& U,const KN_<R>& V,const KN_<R> & Viso, R coef,int iu=0,int iv=0,float *colors=0,int nbcolors=0,bool hsv=true,bool drawborder=true) const ; // Arrow
619 Rd MinMax(const KN_<R>& U,const KN_<R>& V,int j0,int j1,bool bb=true) const ;
620 Rd MinMax(const KN_<R>& U,int j0, bool bb=true) const ;
621 // void destroy() {RefCounter::destroy();}
622 */
isFEMesh() const623 bool isFEMesh() const { return ! NodesOfElement && ( N==1) ;} // to make optim
624 template <class R>
625 KN<R> newSaveDraw(const KN_<R> & U,int composante,int & lg,KN<typename Mesh::RdHat> &Psub,KN<int> &Ksub,int op_U=0) const ;
626
627
628
629 private: // for gibbs
630 int gibbsv (long* ptvoi,long* vois,long* lvois,long* w,long* v);
631 };
632
633 template<class Mesh>
GbaseFElement(const GFESpace<Mesh> & aVh,int k)634 inline GbaseFElement<Mesh>::GbaseFElement( const GFESpace<Mesh> &aVh, int k)
635 : Vh(aVh),T(Vh.Th[k]),tfe(aVh.TFE[k]),N(aVh.N),number(k){}
636
637 template<class Mesh>
GbaseFElement(const GbaseFElement & K,const GTypeOfFE<Mesh> & atfe)638 inline GbaseFElement<Mesh>::GbaseFElement(const GbaseFElement & K, const GTypeOfFE<Mesh> & atfe)
639 : Vh(K.Vh),T(K.T),tfe(&atfe),N(Vh.N),number(K.number){}
640
641 template<class Mesh>
GFElement(const GFESpace<Mesh> * VVh,int k)642 GFElement<Mesh>::GFElement(const GFESpace<Mesh> * VVh,int k)
643 : GbaseFElement<Mesh>(*VVh,k) ,
644 p(this->Vh.PtrFirstNodeOfElement(k)),
645 nb(this->Vh.NbOfNodesInElement(k))
646 {}
647
648 template<class Mesh>
operator [](int i) const649 inline int GFElement<Mesh>::operator[](int i) const {
650 return p ? p[i] : ((&this->T[i])-this->Vh.Th.vertices);}
651
652 template<class Mesh>
operator ()(int i,int df) const653 inline int GFElement<Mesh>::operator()(int i,int df) const {
654 return this->Vh.FirstDFOfNode(p ? p[i] : ((&this->T[i])-this->Vh.Th.vertices)) + df;}
655
656 template<class Mesh>
NbDoF(int i) const657 inline int GFElement<Mesh>::NbDoF(int i) const {
658 int node =p ? p[i] : ((&this->T[i])-this->Vh.Th.vertices);
659 return this->Vh.LastDFOfNode(node)-this->Vh.FirstDFOfNode(node);}
660
661 template<class Mesh>
BF(const RdHat & PHat,RNMK_ & val) const662 inline void GFElement<Mesh>::BF(const RdHat & PHat,RNMK_ & val) const {
663 this->tfe->FB(Fop_D0|Fop_D1,this->Vh.Th,this->T,PHat,val);}
664
665
666 template<class Mesh>
BF(const What_d whatd,const RdHat & PHat,RNMK_ & val) const667 inline void GFElement<Mesh>::BF(const What_d whatd,const RdHat & PHat,RNMK_ & val) const { this->tfe->FB(whatd,this->Vh.Th,this->T,PHat,val);}
668
669 //template<class Mesh>
670 // inline void GFElement<Mesh>::set(InterpolationMatrix<typename Mesh::Element::RdHat> &M) const { this->tfe->set(this->Vh.Th,this->T,&M);}
671
672 template<class Mesh>
operator ()(const RdHat & PHat,const KN_<R> & u,int i,int op) const673 inline R GFElement<Mesh>::operator()(const RdHat & PHat,
674 const KN_<R> & u,int i,int op) const
675 {
676 return (*this->tfe)(*this,PHat,u,i,op);
677 }
678
679 template<class Mesh>
operator ()(const RdHat & PHat,const KN_<complex<R>> & u,int i,int op) const680 inline complex<R> GFElement<Mesh>::operator()(const RdHat & PHat,const KN_<complex<R> > & u,int i,int op) const
681 {
682 complex<double> * pu=u; // pointeur du tableau
683 double *pr = static_cast<double*>(static_cast<void*>(pu));
684
685 const KN_<R> ur(pr,u.n,u.step*2);
686 const KN_<R> ui(pr+1,u.n,u.step*2);
687
688 return complex<R>((*this->tfe)(*this,PHat,ur,i,op),(*this->tfe)(*this,PHat,ui,i,op));
689 }
690
691
692
693 template<class Mesh>
operator ()(const GFElement<Mesh> & K,const RdHat & PHat,const KN_<R> & u,int componante,int op) const694 R GTypeOfFE<Mesh>::operator()(const GFElement<Mesh> & K,const RdHat & PHat,const KN_<R> & u,int componante,int op) const
695 {
696 KNMK<R> fb(NbDoF,N,last_operatortype); // the value for basic fonction
697 KN<R> fk(NbDoF);
698 for (int i=0;i<NbDoF;i++) // get the local value
699 fk[i] = u[K(i)];
700 // get value of basic function
701 FB( 1 << op ,K.Vh.Th,K.T,PHat,fb);
702 R r = (fb('.',componante,op),fk);
703 return r;
704 }
705
706 int nbdf_d(const int ndfitem[4],const int nd[4]);
707 int nbnode_d(const int ndfitem[4],const int nd[4]);
708 int *builddata_d(const int ndfitem[4],const int nd[4]);
709
710
711
712
713 template<class Mesh>
714 class GTypeOfFESum: public GTypeOfFE<Mesh> {
715
716 public:
717 typedef GFElement<Mesh> FElement;
718 typedef typename Mesh::Element Element;
719 typedef typename Element::RdHat RdHat;
720 typedef typename Mesh::Rd Rd;
721 const int k;
722 KN<const GTypeOfFE<Mesh> *> teb;
723 KN<int> NN,DF,comp,numPtInterpolation;
724
725 GTypeOfFESum(const KN< GTypeOfFE<Mesh> const *> & t);
726 GTypeOfFESum(const GFESpace<Mesh> **,int kk);
727 GTypeOfFESum(const GFESpace<Mesh> &,int kk);
728
729 void Build(); // the true constructor
730
731 void init(InterpolationMatrix<RdHat> & M,FElement * pK=0,int odf=0,int ocomp=0,int *pp=0) const;
732 void set(const Mesh & Th,const Element & K,InterpolationMatrix<RdHat> & M,int ocoef,int odf,int *nump ) const; // no change by deflaut
733 void FB(const What_d whatd,const Mesh & Th,const Element & K,const RdHat &PHat, KNMK_<R> & val) const ;
~GTypeOfFESum()734 ~GTypeOfFESum(){}
735 } ;
736
737
738 template<class Mesh>
FB(const What_d whatd,const Mesh & Th,const Element & K,const RdHat & PHat,KNMK_<R> & val) const739 void GTypeOfFESum<Mesh>::FB(const What_d whatd,const Mesh & Th,const Element & K,const RdHat &PHat, KNMK_<R> & val) const
740 {
741 val=0.0;
742 SubArray t(val.K());
743 for (int i=0;i<k;i++)
744 {
745 int j=comp[i];
746 int ni=NN[i];
747 int di=DF[i];
748 int i1=i+1;
749 int nii=NN[i1];
750 int dii=DF[i1];
751 assert(ni<nii && di < dii);
752 RNMK_ v(val(SubArray(dii-di,di),SubArray(nii-ni,ni),t));
753 if (j<=i)
754 teb[i]->FB(whatd,Th,K,PHat,v);
755 else
756 v=val(SubArray(DF[j+1]-DF[j],DF[j]),SubArray(NN[j+1]-NN[j],NN[j]),t);
757 }
758 }
759
760
761
762 template<class RdHat>
763 template<class Mesh>
InterpolationMatrix(const GFESpace<Mesh> & Vh)764 InterpolationMatrix<RdHat>::InterpolationMatrix(const GFESpace<Mesh> &Vh)
765 :
766 N(Vh.N),np(Vh.maxNbPtforInterpolation),ncoef(Vh.maxNbcoefforInterpolation),
767 invariant(Vh.TFE.constant() ? Vh.TFE[0]->invariantinterpolationMatrix: false),
768 k(-1),
769 P(np),
770 comp(ncoef),
771 p(ncoef),
772 dofe(ncoef)
773 {
774 Vh.TFE[0]->GTypeOfFE<Mesh>::init(*this,0,0,0,0);
775 }
776
777 template<class RdHat>
778 template<class Mesh>
InterpolationMatrix(const GTypeOfFE<Mesh> & tef)779 InterpolationMatrix<RdHat>::InterpolationMatrix(const GTypeOfFE<Mesh> & tef)
780 :
781 N(tef.N),np(tef.NbPtforInterpolation),ncoef(tef.NbcoefforInterpolation),
782 invariant(tef.invariantinterpolationMatrix),
783 k(-1),
784 P(np),
785 comp(ncoef),
786 p(ncoef),
787 dofe(ncoef)
788 {
789 // virtual void init(InterpolationMatrix<RdHat> & M,FElement * pK=0,int odf=0,int ocomp=0,int *pp=0) const
790 tef.GTypeOfFE<Mesh>::init(*this,0,0,0,0);
791 }
792
793 template<class RdHat> template<class Mesh>
set(const GFElement<Mesh> & FK)794 void InterpolationMatrix<RdHat>::set(const GFElement<Mesh> & FK)
795 {
796 if(k==FK.number) return;
797 k=FK.number;
798 if(invariant) return;
799 FK.set(*this);
800 //assert(invariant);
801 // a faire ...
802 }
803
804 typedef GTypeOfFE<Mesh3> TypeOfFE3;
805 typedef GTypeOfFE<MeshS> TypeOfFES;
806 typedef GTypeOfFE<MeshL> TypeOfFEL;
807 typedef GFESpace<Mesh3> FESpace3;
808 typedef GFESpace<MeshS> FESpaceS;
809 typedef GFESpace<MeshL> FESpaceL;
810 typedef GFESpace<Mesh2> FESpace2;
811 typedef GFElement<MeshL> FElementL;
812 typedef GFElement<Mesh3> FElement3;
813 typedef GFElement<MeshS> FElementS;
814 typedef GFElement<Mesh2> FElement2;
815 typedef GbaseFElement<Mesh2> baseFElement2;
816 typedef GbaseFElement<Mesh3> baseFElement3;
817 typedef GbaseFElement<MeshS> baseFElementS;
818 typedef GbaseFElement<MeshL> baseFElementL;
819
820 }
821
822
823 #endif
824