1 // -*- Mode : c++ -*-
2 //
3 // SUMMARY :
4 // USAGE :
5 // ORG : pfes3_tefk
6 // AUTHOR : Frederic Hecht
7 // E-MAIL : hecht@ann.jussieu.fr
8 //
9
10 /*
11
12 This file is part of Freefem++
13
14 Freefem++ is free software; you can redistribute it and/or modify
15 it under the terms of the GNU Lesser General Public License as published by
16 the Free Software Foundation; either version 2.1 of the License, or
17 (at your option) any later version.
18
19 Freefem++ is distributed in the hope that it will be useful,
20 but WITHOUT ANY WARRANTY; without even the implied warranty of
21 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 GNU Lesser General Public License for more details.
23
24 You should have received a copy of the GNU Lesser General Public License
25 along with Freefem++; if not, write to the Free Software
26 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
27 */
28 // ATTENTION pfes est la classe qui modelise un pointeur
29 // sur un FESpace (donc un espace d'element fini
30 // mais si la maillage change alors
31 // l'espace change
32 // cf la fonction set qui reconstruit FESpace
33
34 #ifndef lgfem_hpp_
35 #define lgfem_hpp_
36 #include <array_resize.hpp>
37
38 extern Block *currentblock;
39
40 void init_lgmat( ); // initialisation for sparse mat functionnallity
41
42 class v_fes;
43 class v_fes3;
44 class v_fesS;
45 class v_fesL;
46 typedef v_fes *pfes;
47 typedef v_fes3 *pfes3;
48 typedef v_fesS *pfesS;
49 typedef v_fesL *pfesL;
50
51 namespace Fem2D {
52 class Mesh3;
53
54 }
55 using Fem2D::Mesh;
56 using Fem2D::Mesh3;
57 using Fem2D::MeshL;
58 using Fem2D::MeshS;
59
60 typedef const Mesh *pmesh;
61 typedef const Mesh3 *pmesh3;
62 typedef const MeshS *pmeshS;
63 typedef const MeshL *pmeshL;
64
65 using Fem2D::FESpace;
66 using Fem2D::R;
67 using Fem2D::TypeOfFE;
68 // typedef double R;
69 namespace {
70 using namespace Fem2D;
71 using Fem2D::Vertex;
72
73 class lgVertex {
74 public:
75 typedef double R;
76 CountPointer< const Mesh > pTh;
77 Vertex *v;
Check() const78 void Check( ) const {
79 if (!v || !pTh) {
80 ExecError("Too bad! Unset Vertex!");
81 }
82 }
init()83 void init( ) {
84 v = 0;
85 pTh.init( );
86 }
lgVertex(const Mesh * Th,long kk)87 lgVertex(const Mesh *Th, long kk) : pTh(Th), v(&(*pTh)(kk)) {}
lgVertex(const Mesh * Th,Vertex * kk)88 lgVertex(const Mesh *Th, Vertex *kk) : pTh(Th), v(kk) {}
operator int() const89 operator int( ) const {
90 Check( );
91 return (*pTh)(v);
92 }
operator R2*()93 operator R2 *( ) {
94 Check( );
95 return v;
96 }
x() const97 R x( ) const {
98 Check( );
99 return v->x;
100 }
y() const101 R y( ) const {
102 Check( );
103 return v->y;
104 }
105 // R z() const {Check() ; return v->z;}
lab() const106 long lab( ) const {
107 Check( );
108 return v->lab;
109 }
destroy()110 void destroy( ) { pTh.destroy( ); }
111 };
112
113 class lgElement {
114 public:
115 struct Adj { //
116 const Mesh *pTh;
117 Triangle *k;
Adj__anon30b793eb0111::lgElement::Adj118 Adj(const lgElement &pp) : pTh(pp.pTh), k(pp.k) {}
adj__anon30b793eb0111::lgElement::Adj119 lgElement adj(long &e) const {
120 int ee;
121 ffassert(pTh && k && e >= 0 && e < 3);
122 long kk = pTh->ElementAdj((*pTh)(k), ee = e);
123 e = ee;
124 return lgElement(pTh, kk);
125 }
126 };
127 CountPointer< const Mesh > pTh;
128 Triangle *k;
129
lgElement()130 lgElement( ) : k(0) {}
Check() const131 void Check( ) const {
132 if (!k || !pTh) {
133 ExecError("Unset Triangle,Sorry!");
134 }
135 }
init()136 void init( ) {
137 k = 0;
138 pTh.init( );
139 }
destroy()140 void destroy( ) { pTh.destroy( ); }
lgElement(const Mesh * Th,long kk)141 lgElement(const Mesh *Th, long kk) : pTh(Th), k(&(*Th)[kk]) {}
lgElement(const Mesh * Th,Triangle * kk)142 lgElement(const Mesh *Th, Triangle *kk) : pTh(Th), k(kk) {}
operator int() const143 operator int( ) const {
144 Check( );
145 return (*pTh)(k);
146 }
operator [](const long & i) const147 lgVertex operator[](const long &i) const {
148 Check( );
149 return lgVertex(pTh, &(*k)[i]);
150 }
lab() const151 long lab( ) const {
152 Check( );
153 return k ? k->lab : LONG_MIN;
154 }
area() const155 double area( ) const {
156 Check( );
157 return k->area;
158 }
n() const159 long n( ) const { return k ? 3 : 0; }
operator ==(const lgElement & l) const160 bool operator==(const lgElement &l) const { return pTh == l.pTh && k == l.k; }
operator !=(const lgElement & l) const161 bool operator!=(const lgElement &l) const { return pTh != l.pTh || k != l.k; }
operator <(const lgElement & l) const162 bool operator<(const lgElement &l) const { return pTh == l.pTh && k < l.k; }
operator <=(const lgElement & l) const163 bool operator<=(const lgElement &l) const { return pTh == l.pTh && k <= l.k; }
164 };
165 // add FH August 2009 ...
166
167 class lgBoundaryEdge {
168 public:
169 struct BE {
170 const Mesh *p;
BE__anon30b793eb0111::lgBoundaryEdge::BE171 BE(const Mesh *pp) : p(pp) {}
BE__anon30b793eb0111::lgBoundaryEdge::BE172 BE(const Mesh **pp) : p(*pp) {}
operator const Mesh*__anon30b793eb0111::lgBoundaryEdge::BE173 operator const Mesh *( ) const { return p; }
174 };
175
176 CountPointer< const Mesh > pTh;
177 BoundaryEdge *k;
178
lgBoundaryEdge()179 lgBoundaryEdge( ) : k(0) {}
Check() const180 void Check( ) const {
181 if (!k || !pTh) {
182 ExecError("Unset BoundaryEdge,Sorry!");
183 }
184 }
init()185 void init( ) {
186 k = 0;
187 pTh.init( );
188 }
destroy()189 void destroy( ) { pTh.destroy( ); }
lgBoundaryEdge(const Mesh * Th,long kk)190 lgBoundaryEdge(const Mesh *Th, long kk) : pTh(Th), k(&(*pTh).be(kk)) {}
lgBoundaryEdge(const Mesh * Th,BoundaryEdge * kk)191 lgBoundaryEdge(const Mesh *Th, BoundaryEdge *kk) : pTh(Th), k(kk) {}
lgBoundaryEdge(const BE & be,long kk)192 lgBoundaryEdge(const BE &be, long kk) : pTh(be.p), k(&(*pTh).be(kk)) {}
lgBoundaryEdge(const BE & be,BoundaryEdge * kk)193 lgBoundaryEdge(const BE &be, BoundaryEdge *kk) : pTh(be.p), k(kk) {}
operator int() const194 operator int( ) const {
195 Check( );
196 return (*pTh)(k);
197 }
operator [](const long & i) const198 lgVertex operator[](const long &i) const {
199 Check( );
200 return lgVertex(pTh, &(*k)[i]);
201 }
lab() const202 long lab( ) const {
203 Check( );
204 return k ? k->lab : 0;
205 }
length() const206 double length( ) const {
207 Check( );
208 return k->length( );
209 }
n() const210 long n( ) const { return k ? 2 : 0; }
Element() const211 lgElement Element( ) const {
212 Check( );
213 int ee;
214 return lgElement(pTh, (*pTh).BoundaryElement((*pTh)(k), ee));
215 }
EdgeElement() const216 long EdgeElement( ) const {
217 Check( );
218 int ee;
219 (*pTh).BoundaryElement((*pTh)(k), ee);
220 return ee;
221 }
222 };
223
224 } // namespace
225
226 void GetPeriodic(const int d, Expression perio, int &nbcperiodic, Expression *&periodic);
227
228 bool BuildPeriodic(int nbcperiodic, Expression *periodic, const Mesh &Th, Stack stack, int &nbdfv,
229 KN< int > &ndfv, int &nbdfe, KN< int > &ndfe);
230
231 // <<v_fes>> uses [[file:~/ff/src/femlib/RefCounter.hpp::RefCounter]]
232
233 // 2d
234 class v_fes : public RefCounter {
235 public:
236 typedef ::pfes pfes;
237 typedef ::FESpace FESpace;
238
239 static const int dHat = 2, d = 2;
240 const int N;
241 const pmesh *ppTh; // adr du maillage
242 CountPointer< FESpace > pVh;
243 Stack stack; // the stack is use whith periodique expression
244
245 int nbcperiodic;
246 Expression *periodic;
247
operator FESpace*()248 operator FESpace *( ) {
249 throwassert(dHat == 2);
250 if (!pVh || *ppTh != &pVh->Th) pVh = CountPointer< FESpace >(update( ), true);
251 return pVh;
252 }
253
254 FESpace *update( );
255
v_fes(int NN,const pmesh * t,Stack s,int n,Expression * p)256 v_fes(int NN, const pmesh *t, Stack s, int n, Expression *p)
257 : N(NN), ppTh(t), pVh(0), stack(s), nbcperiodic(n), periodic(p) {}
v_fes(int NN,const v_fes * f,Stack s,int n,Expression * p)258 v_fes(int NN, const v_fes *f, Stack s, int n, Expression *p)
259 : N(NN), ppTh(f->ppTh), pVh(0), stack(s), nbcperiodic(n), periodic(p) {}
260
destroy()261 void destroy( ) {
262 ppTh = 0;
263 pVh = 0;
264 delete this;
265 }
~v_fes()266 virtual ~v_fes( ) {}
267 bool buildperiodic(Stack stack, int &nbdfv, KN< int > &ndfv, int &nbdfe, KN< int > &ndfe);
268 virtual FESpace *buildupdate(int &nbdfv, KN< int > &ndfv, int &nbdfe, KN< int > &ndfe) = 0;
269 virtual FESpace *buildupdate( ) = 0;
270 };
271
272 // 3D volume
273 class v_fes3 : public RefCounter {
274 public:
275 typedef pfes3 pfes;
276 typedef FESpace3 FESpace;
277
278 static const int dHat = 3, d = 3;
279 const int N;
280 const pmesh3 *ppTh; // adr du maillage
281 CountPointer< FESpace3 > pVh;
282
283 Stack stack; // the stack is use whith periodique expression
284
285 int nbcperiodic;
286 Expression *periodic;
287
operator FESpace3*()288 operator FESpace3 *( ) {
289 throwassert(dHat == 3);
290 if (!pVh || *ppTh != &pVh->Th) pVh = CountPointer< FESpace3 >(update( ), true);
291 return pVh;
292 }
293 FESpace3 *update( );
294
v_fes3(int NN,const pmesh3 * t,Stack s,int n,Expression * p)295 v_fes3(int NN, const pmesh3 *t, Stack s, int n, Expression *p)
296 : N(NN), ppTh(t), pVh(0), stack(s), nbcperiodic(n), periodic(p) {}
v_fes3(int NN,const v_fes3 * f,Stack s,int n,Expression * p)297 v_fes3(int NN, const v_fes3 *f, Stack s, int n, Expression *p)
298 : N(NN), ppTh(f->ppTh), pVh(0), stack(s), nbcperiodic(n), periodic(p) {}
299
300 // void destroy(){ ppTh=0;pVh=0; delete this;}
~v_fes3()301 virtual ~v_fes3( ) {}
302 bool buildperiodic(Stack stack, KN< int > &ndfe);
buildupdate(KN<int> & ndfe)303 virtual FESpace3 *buildupdate(KN< int > &ndfe) { return 0; }
buildupdate()304 virtual FESpace3 *buildupdate( ) { return 0; };
305 };
306
307 // 3D surface
308 class v_fesS : public RefCounter {
309 public:
310 typedef pfesS pfes;
311 typedef FESpaceS FESpace;
312
313 static const int dHat = 2, d = 3;
314 const int N;
315 const pmeshS *ppTh; // adr du maillage
316
317 CountPointer< FESpaceS > pVh;
318
319 Stack stack; // the stack is use whith periodique expression
320
321 int nbcperiodic;
322 Expression *periodic;
323
operator FESpaceS*()324 operator FESpaceS *( ) {
325 throwassert(dHat == 2);
326 if (!pVh || *ppTh != &pVh->Th) pVh = CountPointer< FESpaceS >(update( ), true);
327 return pVh;
328 }
329 FESpaceS *update( );
330
v_fesS(int NN,const pmeshS * t,Stack s,int n,Expression * p)331 v_fesS(int NN, const pmeshS *t, Stack s, int n, Expression *p) /// TODO
332 : N(NN), ppTh(t), pVh(0), stack(s), nbcperiodic(n), periodic(p) {
333 } // take a pmesh3 and use the pmeshS
v_fesS(int NN,const v_fesS * f,Stack s,int n,Expression * p)334 v_fesS(int NN, const v_fesS *f, Stack s, int n, Expression *p)
335 : N(NN), ppTh(f->ppTh), pVh(0), stack(s), nbcperiodic(n), periodic(p) {}
336
337 // void destroy(){ ppTh=0;pVh=0; delete this;}
~v_fesS()338 virtual ~v_fesS( ) {}
339 bool buildperiodic(Stack stack, KN< int > &ndfe);
buildupdate(KN<int> & ndfe)340 virtual FESpaceS *buildupdate(KN< int > &ndfe) { return 0; }
buildupdate()341 virtual FESpaceS *buildupdate( ) { return 0; };
342 };
343
344 // 3D curve
345 class v_fesL : public RefCounter {
346 public:
347 typedef pfesL pfes;
348 typedef FESpaceL FESpace;
349
350 static const int dHat = 1, d = 3;
351 const int N;
352 const pmeshL *ppTh; // adr du maillage
353
354 CountPointer< FESpaceL > pVh;
355
356 Stack stack; // the stack is use whith periodique expression
357
358 int nbcperiodic;
359 Expression *periodic;
360
operator FESpaceL*()361 operator FESpaceL *( ) {
362 throwassert(dHat == 1);
363 if (!pVh || *ppTh != &pVh->Th) pVh = CountPointer< FESpaceL >(update( ), true);
364 return pVh;
365 }
366 FESpaceL *update( );
367
v_fesL(int NN,const pmeshL * t,Stack s,int n,Expression * p)368 v_fesL(int NN, const pmeshL *t, Stack s, int n, Expression *p) /// TODO
369 : N(NN), ppTh(t), pVh(0), stack(s), nbcperiodic(n), periodic(p) {
370 } // take a pmesh3 and use the pmeshS
v_fesL(int NN,const v_fesL * f,Stack s,int n,Expression * p)371 v_fesL(int NN, const v_fesL *f, Stack s, int n, Expression *p)
372 : N(NN), ppTh(f->ppTh), pVh(0), stack(s), nbcperiodic(n), periodic(p) {}
373
374 // void destroy(){ ppTh=0;pVh=0; delete this;}
~v_fesL()375 virtual ~v_fesL( ) {}
376 bool buildperiodic(Stack stack, KN< int > &ndfe);
buildupdate(KN<int> & ndfe)377 virtual FESpaceL *buildupdate(KN< int > &ndfe) { return 0; }
buildupdate()378 virtual FESpaceL *buildupdate( ) { return 0; };
379 };
380
381 // 2d
382 class pfes_tef : public v_fes {
383 public:
384 const TypeOfFE *tef;
pfes_tef(const pmesh * t,const TypeOfFE * tt,Stack s=NullStack,int n=0,Expression * p=0)385 pfes_tef(const pmesh *t, const TypeOfFE *tt, Stack s = NullStack, int n = 0, Expression *p = 0)
386 : v_fes(tt->N, t, s, n, p), tef(tt) {
387 operator FESpace *( );
388 }
buildupdate(int & nbdfv,KN<int> & ndfv,int & nbdfe,KN<int> & ndfe)389 FESpace *buildupdate(int &nbdfv, KN< int > &ndfv, int &nbdfe, KN< int > &ndfe) {
390 return *ppTh ? new FESpace(**ppTh, *tef, nbdfv, (int *)ndfv, nbdfe, (int *)ndfe) : 0;
391 }
buildupdate()392 FESpace *buildupdate( ) { return *ppTh ? new FESpace(**ppTh, *tef) : 0; }
393 };
394
395 class pfes_tefk : public v_fes {
396 public:
397 const TypeOfFE **tef;
398 const int k;
pfes_tefk(const pmesh * t,const TypeOfFE ** tt,int kk,Stack s=NullStack,int n=0,Expression * p=0)399 pfes_tefk(const pmesh *t, const TypeOfFE **tt, int kk, Stack s = NullStack, int n = 0,
400 Expression *p = 0)
401 : v_fes(sum(tt, &Fem2D::TypeOfFE::N, kk), t, s, n, p), tef(tt), k(kk) {
402 // cout << "pfes_tefk const" << tef << " " << this << endl;
403 operator FESpace *( );
404 }
buildupdate()405 FESpace *buildupdate( ) {
406 // cout << "pfes_tefk upd:" << tef << " " << this << endl;
407 assert(tef);
408 return *ppTh ? new FESpace(**ppTh, tef, k) : 0;
409 }
~pfes_tefk()410 virtual ~pfes_tefk( ) { delete[] tef; }
buildupdate(int & nbdfv,KN<int> & ndfv,int & nbdfe,KN<int> & ndfe)411 FESpace *buildupdate(int &nbdfv, KN< int > &ndfv, int &nbdfe, KN< int > &ndfe) {
412 assert(tef);
413 return *ppTh ? new FESpace(**ppTh, tef, k, nbdfv, ndfv, nbdfe, ndfe) : 0;
414 }
415 };
416
417 // 3D volume
418 class pfes3_tef : public v_fes3 {
419 public:
420 const TypeOfFE3 *tef;
pfes3_tef(const pmesh3 * t,const TypeOfFE3 * tt,Stack s=NullStack,int n=0,Expression * p=0)421 pfes3_tef(const pmesh3 *t, const TypeOfFE3 *tt, Stack s = NullStack, int n = 0, Expression *p = 0)
422 : v_fes3(tt->N, t, s, n, p), tef(tt) {
423 operator FESpace3 *( );
424 }
buildupdate(KN<int> & ndfe)425 FESpace3 *buildupdate(KN< int > &ndfe) {
426 return *ppTh ? new FESpace3(**ppTh, *tef, ndfe.size( ) / 2, ndfe) : 0;
427 }
buildupdate()428 FESpace3 *buildupdate( ) { return *ppTh ? new FESpace3(**ppTh, *tef) : 0; }
429 };
430
431 class pfes3_tefk : public v_fes3 {
432 public:
433 const TypeOfFE3 **tef;
434 const int k;
435 KN< GTypeOfFE< Mesh3 > const * > atef;
436 GTypeOfFESum< Mesh3 > tefs;
437
sum(const Fem2D::TypeOfFE3 ** l,int const Fem2D::TypeOfFE3::* p,int n)438 static int sum(const Fem2D::TypeOfFE3 **l, int const Fem2D::TypeOfFE3::*p, int n) {
439 int r = 0;
440 for (int i = 0; i < n; i++) r += l[i]->*p;
441 return r;
442 }
443
pfes3_tefk(const pmesh3 * t,const Fem2D::TypeOfFE3 ** tt,int kk,Stack s=NullStack,int n=0,Expression * p=0)444 pfes3_tefk(const pmesh3 *t, const Fem2D::TypeOfFE3 **tt, int kk, Stack s = NullStack, int n = 0,
445 Expression *p = 0)
446 : v_fes3(sum((const Fem2D::TypeOfFE3 **)tt, &Fem2D::TypeOfFE3::N, kk), t, s, n, p), tef(tt),
447 k(kk), atef(kk, tt), tefs(atef)
448
449 {
450 // cout << "pfes_tefk const" << tef << " " << this << endl;
451 operator FESpace3 *( );
452 }
buildupdate()453 FESpace3 *buildupdate( ) {
454 // cout << "pfes_tefk upd:" << tef << " " << this << endl;
455 // assert(tef);
456 return *ppTh ? new FESpace3(**ppTh, tefs) : 0;
457 }
~pfes3_tefk()458 virtual ~pfes3_tefk( ) { delete[] tef; }
buildupdate(KN<int> & ndfe)459 FESpace3 *buildupdate(KN< int > &ndfe) {
460 return *ppTh ? new FESpace3(**ppTh, tefs, ndfe.size( ) / 2, ndfe) : 0;
461 }
462 };
463
464 // 3D surface
465 class pfesS_tef : public v_fesS {
466 public:
467 const TypeOfFES *tef;
pfesS_tef(const pmeshS * t,const TypeOfFES * tt,Stack s=NullStack,int n=0,Expression * p=0)468 pfesS_tef(const pmeshS *t, const TypeOfFES *tt, Stack s = NullStack, int n = 0, Expression *p = 0)
469 : v_fesS(tt->N, t, s, n, p), tef(tt) {
470 operator FESpaceS *( );
471 }
buildupdate(KN<int> & ndfe)472 FESpaceS *buildupdate(KN< int > &ndfe) {
473 return *ppTh ? new FESpaceS(**ppTh, *tef, ndfe.size( ) / 2, ndfe) : 0;
474 }
buildupdate()475 FESpaceS *buildupdate( ) { return *ppTh ? new FESpaceS(**ppTh, *tef) : 0; }
476 };
477
478 class pfesS_tefk : public v_fesS {
479 public:
480 const TypeOfFES **tef;
481 const int k;
482 KN< GTypeOfFE< MeshS > const * > atef;
483 GTypeOfFESum< MeshS > tefs;
484
sum(const Fem2D::TypeOfFES ** l,int const Fem2D::TypeOfFES::* p,int n)485 static int sum(const Fem2D::TypeOfFES **l, int const Fem2D::TypeOfFES::*p, int n) {
486 int r = 0;
487 for (int i = 0; i < n; i++) r += l[i]->*p;
488 return r;
489 }
490
pfesS_tefk(const pmeshS * t,const Fem2D::TypeOfFES ** tt,int kk,Stack s=NullStack,int n=0,Expression * p=0)491 pfesS_tefk(const pmeshS *t, const Fem2D::TypeOfFES **tt, int kk, Stack s = NullStack, int n = 0,
492 Expression *p = 0)
493 : v_fesS(sum((const Fem2D::TypeOfFES **)tt, &Fem2D::TypeOfFES::N, kk), t, s, n, p), tef(tt),
494 k(kk), atef(kk, tt), tefs(atef)
495
496 {
497 // cout << "pfes_tefk const" << tef << " " << this << endl;
498 operator FESpaceS *( );
499 }
buildupdate()500 FESpaceS *buildupdate( ) {
501 // cout << "pfes_tefk upd:" << tef << " " << this << endl;
502 // assert(tef);
503 return *ppTh ? new FESpaceS(**ppTh, tefs) : 0;
504 }
~pfesS_tefk()505 virtual ~pfesS_tefk( ) { delete[] tef; }
buildupdate(KN<int> & ndfe)506 FESpaceS *buildupdate(KN< int > &ndfe) {
507 return *ppTh ? new FESpaceS(**ppTh, tefs, ndfe.size( ) / 2, ndfe) : 0;
508 }
509 };
510
511 // 3D curve
512 class pfesL_tef : public v_fesL {
513 public:
514 const TypeOfFEL *tef;
pfesL_tef(const pmeshL * t,const TypeOfFEL * tt,Stack s=NullStack,int n=0,Expression * p=0)515 pfesL_tef(const pmeshL *t, const TypeOfFEL *tt, Stack s = NullStack, int n = 0, Expression *p = 0)
516 : v_fesL(tt->N, t, s, n, p), tef(tt) {
517 operator FESpaceL *( );
518 }
buildupdate(KN<int> & ndfe)519 FESpaceL *buildupdate(KN< int > &ndfe) {
520 return *ppTh ? new FESpaceL(**ppTh, *tef, ndfe.size( ) / 2, ndfe) : 0;
521 }
buildupdate()522 FESpaceL *buildupdate( ) { return *ppTh ? new FESpaceL(**ppTh, *tef) : 0; }
523 };
524
525 class pfesL_tefk : public v_fesL {
526 public:
527 const TypeOfFEL **tef;
528 const int k;
529 KN< GTypeOfFE< MeshL > const * > atef;
530 GTypeOfFESum< MeshL > tefs;
531
sum(const Fem2D::TypeOfFEL ** l,int const Fem2D::TypeOfFEL::* p,int n)532 static int sum(const Fem2D::TypeOfFEL **l, int const Fem2D::TypeOfFEL::*p, int n) {
533 int r = 0;
534 for (int i = 0; i < n; i++) r += l[i]->*p;
535 return r;
536 }
537
pfesL_tefk(const pmeshL * t,const Fem2D::TypeOfFEL ** tt,int kk,Stack s=NullStack,int n=0,Expression * p=0)538 pfesL_tefk(const pmeshL *t, const Fem2D::TypeOfFEL **tt, int kk, Stack s = NullStack, int n = 0,
539 Expression *p = 0)
540 : v_fesL(sum((const Fem2D::TypeOfFEL **)tt, &Fem2D::TypeOfFEL::N, kk), t, s, n, p), tef(tt),
541 k(kk), atef(kk, tt), tefs(atef)
542
543 {
544 // cout << "pfes_tefk const" << tef << " " << this << endl;
545 operator FESpaceL *( );
546 }
buildupdate()547 FESpaceL *buildupdate( ) {
548 // cout << "pfes_tefk upd:" << tef << " " << this << endl;
549 // assert(tef);
550 return *ppTh ? new FESpaceL(**ppTh, tefs) : 0;
551 }
~pfesL_tefk()552 virtual ~pfesL_tefk( ) { delete[] tef; }
buildupdate(KN<int> & ndfe)553 FESpaceL *buildupdate(KN< int > &ndfe) {
554 return *ppTh ? new FESpaceL(**ppTh, tefs, ndfe.size( ) / 2, ndfe) : 0;
555 }
556 };
557
558 class pfes_fes : public v_fes {
559 public:
560 pfes *Vh;
561 int n;
pfes_fes(pfes * Vhh,int nn,Stack s=NullStack,int n=0,Expression * p=0)562 pfes_fes(pfes *Vhh, int nn, Stack s = NullStack, int n = 0, Expression *p = 0)
563 : v_fes((**Vhh).N * nn, static_cast< const v_fes * >(*Vhh), s, n, p), Vh(Vhh), n(nn) {
564 operator FESpace *( );
565 };
buildupdate()566 FESpace *buildupdate( ) { return new FESpace(*(FESpace *)**Vh, n); }
buildupdate(int & nbdfv,KN<int> & ndfv,int & nbdfe,KN<int> & ndfe)567 FESpace *buildupdate(int &nbdfv, KN< int > &ndfv, int &nbdfe, KN< int > &ndfe) {
568 InternalError(" No way to define a periodic BC in this case: tensorisation of FEspace ");
569 // return new FESpace(***Vh,n,nbdfv,ndfv,nbdfe,ndfe);
570 }
571 };
572
573 template< class K, class v_fes >
574 class FEbase; // <<FEbase>>
575 template< class K, class v_fes >
576 class FEcomp {
577 public:
578 typedef typename v_fes::pfes pfes;
579 typedef typename v_fes::FESpace FESpace;
580
581 friend class FEbase< K, v_fes >;
582 FEbase< K, v_fes > *base;
583 int comp;
FEcomp(FEbase<K,v_fes> * b,int c)584 FEcomp(FEbase< K, v_fes > *b, int c) : base(b), comp(c){};
585
586 private: // rule of programming
587 FEcomp(const FEcomp &);
588 void operator=(const FEcomp &);
589 };
590
591 template< class K, class v_fes >
592 class FEbase {
593 public:
594 typedef typename v_fes::pfes pfes;
595 typedef typename v_fes::FESpace FESpace;
596
597 v_fes *const *pVh; // pointeur sur la variable stockant FESpace;
598 KN< K > *xx; // value
599 CountPointer< FESpace > Vh; // espace courant
600
x()601 KN< K > *x( ) { return xx; }
602
FEbase(const pfes * ppVh)603 FEbase(const pfes *ppVh) : pVh(ppVh), xx(0), Vh(0) {}
604
~FEbase()605 ~FEbase( ) { delete xx; }
destroy()606 void destroy( ) { // cout << "~FEbase destroy " << this << endl;
607 delete this;
608 }
609
operator =(KN<K> * y)610 void operator=(KN< K > *y) {
611 Vh = **pVh;
612 throwassert((bool)Vh);
613 if (xx) delete xx;
614 xx = y;
615 ffassert(y->N( ) == Vh->NbOfDF);
616 }
617
operator =(KN_<K> & y)618 void operator=(KN_< K > &y) {
619 Vh = **pVh;
620 throwassert((bool)Vh);
621 if (xx) { // resize if need
622 if (xx->N( ) != Vh->NbOfDF) delete xx;
623 xx = 0;
624 }
625 if (!xx) xx = new KN< K >(Vh->NbOfDF);
626 ffassert(SameShape(y, *xx));
627 *xx = y;
628 }
629
newVh()630 FESpace *newVh( ) {
631 throwassert(pVh);
632 const pfes pp = *pVh;
633 // cout << pVh << " " << *pVh << endl;
634 return *pp;
635 }
636
operator FESpace&()637 operator FESpace &( ) {
638 throwassert(Vh);
639 return *Vh;
640 }
641
642 private: // rule of programming
643 FEbase(const FEbase &);
644 void operator=(const FEbase &);
645 };
646
647 template< class K >
648 class FEbaseArrayKn {
649 public: // for eigen value
650 int N;
FEbaseArrayKn(int NN)651 FEbaseArrayKn(int NN) : N(NN) {}
652 virtual void set(int i, KN_< K >) = 0;
653 virtual KN< K > *get(int i) const = 0; // for P. Jolivet
654 virtual void resize(int i) = 0; // for P. Jolivet
655 };
656
657 template< class K, class v_fes >
658 class FEbaseArray : public FEbaseArrayKn< K > {
659 public:
660 typedef typename v_fes::pfes pfes;
661 typedef typename v_fes::FESpace FESpace;
662
663 // int N;
664 FEbase< K, v_fes > **xx;
FEbaseArray(const pfes * ppVh,int NN)665 FEbaseArray(const pfes *ppVh, int NN)
666 : FEbaseArrayKn< K >(NN), xx(new FEbase< K, v_fes > *[std::max(NN, 1)]) {
667 for (int i = 0; i < std::max(this->N, 1); i++) xx[i] = new FEbase< K, v_fes >(ppVh);
668 }
~FEbaseArray()669 ~FEbaseArray( ) {
670 // cout << " ~FEbaseArray " << endl;
671 for (int i = 0; i < std::max(this->N, 1); i++) xx[i]->destroy( );
672 delete[] xx;
673 }
destroy()674 void destroy( ) { // cout << " destroy ~FEbaseArray " << endl;
675 delete this;
676 }
operator [](int i) const677 FEbase< K, v_fes > **operator[](int i) const {
678 if (xx == 0 || i < 0 || i >= this->N) ExecError("Out of bound in FEbaseArray");
679 return xx + i;
680 }
681
resize(int i)682 void resize(int i) {
683 if (xx != 0 && i > 0 && i != this->N) {
684 FEbase< K, v_fes > **yy = new FEbase< K, v_fes > *[i];
685 if (i > this->N) {
686 for (unsigned int j = 0; j < std::max(this->N, 1); ++j) yy[j] = xx[j];
687 for (unsigned int j = std::max(this->N, 1); j < i; ++j)
688 yy[j] = new FEbase< K, v_fes >(xx[0]->pVh);
689 } else {
690 for (unsigned int j = 0; j < i; ++j) yy[j] = xx[j];
691 for (unsigned int j = i; j < this->N; ++j) xx[j]->destroy( );
692 }
693 FEbase< K, v_fes > **oldXx = this->xx;
694 this->xx = yy;
695 delete[] oldXx;
696 this->N = i;
697 }
698 }
699
set(int i,KN_<K> v)700 void set(int i, KN_< K > v) { **(operator[](i)) = v; }
701
get(int i) const702 KN< K > *get(int i) const { return (**(operator[](i))).xx; }
703
704 private: // rule of programming
705 FEbaseArray(const FEbaseArray &);
706 void operator=(const FEbaseArray &);
707 };
708
709 void GetPeriodic(const int d, Expression perio, int &nbcperiodic, Expression *&periodic);
710 int GetPeriodic(Expression bb, Expression &b, Expression &f);
711 int GetPeriodic(Expression bb, Expression &b, Expression &f1, Expression &f2);
712
713 C_F0 NewFEarray(const char *id, Block *currentblock, C_F0 &fespacetype, CC_F0 init, bool cplx,
714 int dim);
715 C_F0 NewFEarray(ListOfId *ids, Block *currentblock, C_F0 &fespacetype, CC_F0 init, bool cplx,
716 int dim);
717 C_F0 NewFEvariable(const char *id, Block *currentblock, C_F0 &fespacetype, CC_F0 init, bool cplx,
718 int dim);
719 C_F0 NewFEvariable(ListOfId *ids, Block *currentblock, C_F0 &fespacetype, CC_F0 init, bool cplx,
720 int dim);
NewFEvariable(const char * id,Block * currentblock,C_F0 & fespacetype,bool cplx,int dim)721 inline C_F0 NewFEvariable(const char *id, Block *currentblock, C_F0 &fespacetype, bool cplx,
722 int dim) {
723 CC_F0 init;
724 init = 0;
725 return NewFEvariable(id, currentblock, fespacetype, init, cplx, dim);
726 }
727
NewFEvariable(ListOfId * ids,Block * currentblock,C_F0 & fespacetype,bool cplx,int dim)728 inline C_F0 NewFEvariable(ListOfId *ids, Block *currentblock, C_F0 &fespacetype, bool cplx,
729 int dim) {
730 CC_F0 init;
731 init = 0;
732 return NewFEvariable(ids, currentblock, fespacetype, init, cplx, dim);
733 }
734
735 size_t dimFESpaceImage(const basicAC_F0 &args);
736 aType typeFESpace(const basicAC_F0 &args);
737
738 template< class K, class vv_fes, class FE = FEbase< K, vv_fes > >
739 class E_FEcomp : public E_F0mps {
740 public:
741 // typedef FEbase<K,vv_fes> FE;
742 typedef vv_fes v_fes;
743 typedef pair< FE *, int > Result;
744 Expression a0;
745 const int comp, N;
746
operator ( )(Stack s) const747 AnyType operator( )(Stack s) const {
748 return SetAny< Result >(Result(*GetAny< FE ** >((*a0)(s)), comp));
749 }
E_FEcomp(const C_F0 & x,const int cc,int NN)750 E_FEcomp(const C_F0 &x, const int cc, int NN) : a0(x.LeftValue( )), comp(cc), N(NN) {
751 if (x.left( ) != atype< FE ** >( ))
752 cout << "E_FEcomp: Bug " << *x.left( ) << " != " << *atype< FE ** >( ) << " case "
753 << typeid(K).name( ) << endl;
754 // CompileError("E_FEcomp: Bug ?");
755 throwassert(x.left( ) == atype< FE ** >( ) && a0);
756 }
operator aType() const757 operator aType( ) const { return atype< Result >( ); }
758 };
759
760 // typedef double R;
761 typedef pair< FEbase< double, v_fes > *, int > aFEvarR;
762 typedef pair< FEbaseArray< double, v_fes > *, int > aFEArrayR;
763 typedef pair< FEbase< Complex, v_fes > *, int > aFEvarC;
764
765 template< class K >
766 class interpolate_f_X_1 : public OneOperator {
767 public:
768 // to interpolate f o X^{-1}
769 typedef FEbase< K, v_fes > *pferbase;
770 typedef FEbaseArray< K, v_fes > *pferbasearray;
771 typedef pair< pferbase, int > pfer;
772 typedef pair< pferbasearray, int > pferarray;
773 class type {};
774 class CODE : public E_F0mps {
775 public:
776 Expression f;
777 Expression x, y;
778
CODE(const basicAC_F0 & args)779 CODE(const basicAC_F0 &args) : f(args[0].LeftValue( )) {
780 const E_Array *X(dynamic_cast< const E_Array * >(args[1].LeftValue( )));
781 if (!X || X->size( ) != 2) {
782 CompileError("array of 2 double x,y f([xx,yy]) = ... ");
783 }
784 x = to< double >((*X)[0]);
785 y = to< double >((*X)[1]);
786 assert(f && x && y);
787 }
operator ( )(Stack) const788 AnyType operator( )(Stack) const {
789 ExecError("No evaluation");
790 return Nothing;
791 }
operator aType() const792 operator aType( ) const { return atype< void >( ); }
793
794 }; // end class CODE
795
code(const basicAC_F0 & args) const796 E_F0 *code(const basicAC_F0 &args) const { return new CODE(args); }
797
interpolate_f_X_1()798 interpolate_f_X_1( )
799 : OneOperator(map_type[typeid(type).name( )], map_type[typeid(pfer).name( )],
800 map_type[typeid(E_Array).name( )]) {}
801 };
802
update()803 inline FESpace *v_fes::update( ) {
804 assert(dHat == 2);
805 if (!*ppTh) return 0;
806 if (nbcperiodic) {
807 assert(periodic);
808 const Mesh &Th(**ppTh);
809 KN< int > ndfv(Th.nv);
810 KN< int > ndfe(Th.neb);
811 int nbdfv, nbdfe;
812 bool ret = buildperiodic(stack, nbdfv, ndfv, nbdfe, ndfe);
813 return ret ? buildupdate(nbdfv, ndfv, nbdfe, ndfe) : buildupdate( );
814 } else
815 return buildupdate( );
816 }
817
update()818 inline FESpace3 *v_fes3::update( ) {
819 assert(dHat == 3);
820 if (!*ppTh) return 0;
821 if (nbcperiodic) {
822 assert(periodic);
823 KN< int > ndfe;
824 bool ret = buildperiodic(stack, ndfe);
825 return ret ? buildupdate(ndfe) : buildupdate( );
826 } else
827 return buildupdate( );
828 }
829
update()830 inline FESpaceS *v_fesS::update( ) {
831 assert(dHat == 2);
832 if (!*ppTh) return 0;
833 if (nbcperiodic) {
834 assert(periodic);
835 KN< int > ndfe;
836 bool ret = buildperiodic(stack, ndfe);
837 return ret ? buildupdate(ndfe) : buildupdate( );
838 } else
839 return buildupdate( );
840 }
841
update()842 inline FESpaceL *v_fesL::update( ) {
843 assert(dHat == 1);
844 if (!*ppTh) return 0;
845 if (nbcperiodic) {
846 assert(periodic);
847 KN< int > ndfe;
848 bool ret = buildperiodic(stack, ndfe);
849 return ret ? buildupdate(ndfe) : buildupdate( );
850 } else
851 return buildupdate( );
852 }
853
854 class TabFuncArg {
855 public:
856 typedef double R;
857 typedef Fem2D::R2 R2;
858 Stack s;
859 int nb;
860 Expression *e;
eval(R * v) const861 void eval(R *v) const {
862 for (int i = 0; i < nb; i++)
863 if (e[i]) {
864 v[i] = GetAny< R >((*e[i])(s));
865 // cout <<" (" << i << " " << *v << ") ";
866 } else
867 v[i] = 0;
868 }
eval_X() const869 R2 eval_X( ) const { return R2(GetAny< R >((*e[nb - 2])(s)), GetAny< R >((*e[nb - 1])(s))); }
eval_2(R * v) const870 void eval_2(R *v) const {
871 for (int i = 0; i < nb - 2; i++)
872 if (e[i])
873 v[i] = GetAny< R >((*e[i])(s));
874 else
875 v[i] = 0;
876 }
TabFuncArg(Stack ss,int n)877 TabFuncArg(Stack ss, int n) : s(ss), nb(n), e(new Expression[n]) {}
operator =(int j)878 void operator=(int j) {
879 ffassert(j == 0);
880 for (int i = 0; i < nb; i++) e[i] = 0;
881 } // resert
operator [](int i)882 Expression &operator[](int i) { return e[i]; }
~TabFuncArg()883 ~TabFuncArg( ) { delete[] e; }
884
885 private: // no copy
886 TabFuncArg(const TabFuncArg &);
887 void operator=(const TabFuncArg &);
888 };
889
890 template< class K >
891 class E_F_StackF0F0opt2 : public E_F0mps {
892 public:
893 typedef AnyType (*func)(Stack, Expression, Expression);
894 func f;
895 Expression a0, a1;
E_F_StackF0F0opt2(func ff,Expression aa0,Expression aa1)896 E_F_StackF0F0opt2(func ff, Expression aa0, Expression aa1) : f(ff), a0(aa0), a1(aa1) {
897
898 const E_FEcomp< K, v_fes > *e = dynamic_cast< const E_FEcomp< K, v_fes > * >(a0);
899 if (e && e->N != 1) {
900 cerr << " err interpolation of no scalar FE function componant (" << e->comp << " in 0.."
901 << e->N << ") << with scalar function \n";
902 CompileError("interpolation of no scalar FE function componant with scalar function ");
903 }
904
905 /* // inutil car a0 est un composante d'un vecteur ????
906 // pour l'instant on a juste une erreur a l'execution
907 // et non a la compilation.
908
909 if (a0->nbitem() != a1->nbitem() )
910 { // bofbof
911 cerr << " err interpolation of no scalar FE function ("<< a0->nbitem() << " != " <<
912 a1->nbitem() << ") << with scalar function \n" ; CompileError("interpolation of no scalar
913 FE function with scalar function ");
914 } */
915 deque< pair< Expression, int > > ll;
916 MapOfE_F0 m;
917 size_t top = currentblock->OffSet(0), topbb = top; // FH. bofbof ???
918 int ret = aa1->Optimize(ll, m, top);
919 a1 = new E_F0_Optimize(ll, m, ret);
920 currentblock->OffSet(top - topbb);
921 }
operator ( )(Stack s) const922 AnyType operator( )(Stack s) const { return (*f)(s, a0, a1); }
923 };
924
925 template< class R >
ShowBound(const KN<R> & y,ostream & f)926 inline ostream &ShowBound(const KN< R > &y, ostream &f) {
927 f << " -- vector function's bound " << y.min( ) << " " << y.max( );
928 return f;
929 }
930 template<>
ShowBound(const KN<Complex> & y,ostream & f)931 inline ostream &ShowBound< Complex >(const KN< Complex > &y, ostream &f) {
932 f << " -- vector function's bound : (no complex Value) ";
933 return f;
934 }
935
936 template< class K >
937 class Op3_K2R : public ternary_function< K, R, R, K > {
938 public:
939 class Op : public E_F0mps {
940 public:
941 Expression a, b, c;
Op(Expression aa,Expression bb,Expression cc)942 Op(Expression aa, Expression bb, Expression cc) : a(aa), b(bb), c(cc) {}
operator ( )(Stack s) const943 AnyType operator( )(Stack s) const {
944 R xx(GetAny< R >((*b)(s)));
945 R yy(GetAny< R >((*c)(s)));
946 MeshPoint &mp = *MeshPointStack(s), mps = mp;
947 mp.set(xx, yy, 0.0);
948 AnyType ret = (*a)(s);
949 mp = mps;
950 return ret;
951 }
952 };
953 };
954 template< class K >
955 class Op4_K2R : public quad_function< K, R, R, R, K > {
956 public:
957 class Op : public E_F0mps {
958 public:
959 Expression a, b, c, d;
Op(Expression aa,Expression bb,Expression cc,Expression dd)960 Op(Expression aa, Expression bb, Expression cc, Expression dd) : a(aa), b(bb), c(cc), d(dd) {}
operator ( )(Stack s) const961 AnyType operator( )(Stack s) const {
962 // cout <<"class Op4_K2R : public quad_function<K,R,R,R,K>" << endl;
963 R xx(GetAny< R >((*b)(s)));
964 R yy(GetAny< R >((*c)(s)));
965 R zz(GetAny< R >((*d)(s)));
966 MeshPoint &mp = *MeshPointStack(s), mps = mp;
967 mp.set(xx, yy, zz);
968 AnyType ret = (*a)(s);
969 mp = mps;
970 return ret;
971 }
972 };
973 };
974
975 // 3D
976 template< class K, class v_fes >
977 class E_set_fev3 : public E_F0mps {
978 public:
979 typedef typename v_fes::pfes pfes;
980 typedef typename v_fes::FESpace FESpace;
981 typedef typename FESpace::Mesh Mesh;
982 typedef typename FESpace::FElement FElement;
983 typedef typename Mesh::Element Element;
984 typedef typename Mesh::Vertex Vertex;
985 typedef typename Mesh::RdHat RdHat;
986 typedef typename Mesh::Rd Rd;
987
988 E_Array aa;
989 Expression ppfe;
990 bool optimize, optimizecheck;
991
992 vector< size_t > where_in_stack_opt;
993 Expression optiexp0, optiexpK;
994
995 E_set_fev3(const E_Array *a, Expression pp);
996
997 AnyType operator( )(Stack) const;
operator aType() const998 operator aType( ) const { return atype< void >( ); }
999 };
1000
1001 // 2d
1002 template< class K >
1003 class E_set_fev : public E_F0mps {
1004 public:
1005 const int dim;
1006 E_Array aa;
1007 Expression ppfe;
1008 bool optimize, optimizecheck;
1009 vector< size_t > where_in_stack_opt;
1010 Expression optiexp0, optiexpK;
1011
1012 E_set_fev(const E_Array *a, Expression pp, int ddim = 2);
1013
1014 AnyType operator( )(Stack) const;
1015 AnyType Op2d(Stack) const;
1016 AnyType Op3d(Stack) const;
1017 AnyType OpS(Stack) const;
1018
operator aType() const1019 operator aType( ) const { return atype< void >( ); }
1020 };
1021
1022 bool InCircularList(const int *p, int i, int k);
1023 template< class T >
1024 int numeroteclink(KN_< T > &ndfv);
1025
1026 // for ... vectorial FE array ..
1027
1028 template< class K, class v_fes >
1029 class E_FEcomp_get_elmnt_array : public E_F0 {
1030 public:
1031 typedef FEbaseArray< K, v_fes > *pfekbasearray;
1032 typedef FEbase< K, v_fes > *pfekbase;
1033 typedef pair< pfekbase, int > pfek;
1034 typedef pair< pfekbasearray, int > pfekarray;
1035 typedef pfek R;
1036 typedef pfekarray A;
1037 typedef long B;
1038 typedef FEbaseArray< K, v_fes > CFE;
1039 typedef E_FEcomp< K, v_fes, CFE > E_KFEArray;
1040
1041 Expression a0, a1;
1042 const E_KFEArray *a00;
1043 const int comp, N;
1044
E_FEcomp_get_elmnt_array(Expression aa0,Expression aa1,int compp,int NN,const E_KFEArray * aa00)1045 E_FEcomp_get_elmnt_array(Expression aa0, Expression aa1, int compp, int NN,
1046 const E_KFEArray *aa00)
1047 : a0(aa0), a1(aa1), a00(aa00), comp(compp), N(NN) {}
operator ( )(Stack s) const1048 AnyType operator( )(Stack s) const {
1049 return SetAny< R >(get_element(GetAny< A >((*a0)(s)), GetAny< B >((*a1)(s))));
1050 }
MeshIndependent() const1051 bool MeshIndependent( ) const { return a0->MeshIndependent( ) && a1->MeshIndependent( ); } //
1052 };
1053
1054 template< class K, class v_fes >
1055 class OneOperator2_FE_get_elmnt : public OneOperator {
1056 // /// Add<pferbasearray*>("[","",new
1057 // OneOperator2_FEcomp<pferbase*,pferbasearray*,long>(get_element)); // not used ...
1058 // Add<pferarray>("[","",new OneOperator2_<pfer,pferarray,long>(get_element));
1059 typedef FEbase< K, v_fes > *pfekbase;
1060 typedef FEbaseArray< K, v_fes > *pfekbasearray;
1061 typedef pair< pfekbase, int > pfek;
1062 typedef pair< pfekbasearray, int > pfekarray;
1063
1064 typedef pfek R;
1065 typedef pfekarray A;
1066 typedef long B;
1067 typedef E_FEcomp_get_elmnt_array< K, v_fes > CODE;
1068 typedef FEbaseArray< K, v_fes > CFE;
1069 typedef E_FEcomp< K, v_fes, CFE > E_KFEArray;
1070 typedef E_FEcomp< K, v_fes > E_KFEi;
1071
1072 aType t0, t1; // return type type de f, f(t1, t2)
1073 public:
code(const basicAC_F0 & args) const1074 E_F0 *code(const basicAC_F0 &args) const {
1075 const E_KFEArray *afe = dynamic_cast< const E_KFEArray * >(args[0].LeftValue( ));
1076 // cout << " build xxx " << afe << " " << *args[0].left() << endl;
1077 // afe=0;// E_KFEi n'est pas le bon type on mame le vieux code ?????
1078 ffassert(afe);
1079 return new CODE(t0->CastTo(args[0]), t1->CastTo(args[1]), afe->comp, afe->N, afe);
1080 }
OneOperator2_FE_get_elmnt()1081 OneOperator2_FE_get_elmnt( )
1082 : OneOperator(map_type[typeid(R).name( )], map_type[typeid(A).name( )],
1083 map_type[typeid(B).name( )]),
1084 t0(map_type[typeid(A).name( )]), t1(map_type[typeid(B).name( )]) {}
1085 };
1086
1087 template< typename KK, typename vv_fes, typename CC >
1088 struct FFset3 {
1089 typedef KK K;
1090 typedef vv_fes v_fes;
1091 typedef vv_fes *pfes;
1092 typedef FEbase< K, vv_fes > **R;
1093 typedef R A;
1094 typedef pfes *B;
1095 typedef CC C;
1096 typedef Expression MC;
CloneFFset31097 static Expression Clone(Expression e) { return e; }
CheckFFset31098 static bool Check(MC l) { return true; }
fFFset31099 static void f(KN< K > *x, MC a, A pp, Stack stack) {
1100 CC at = GetAny< C >((*a)(stack));
1101 *x = at;
1102 }
1103 };
1104 template< typename KK, typename vv_fes, typename CC >
1105 struct FFset3call {
1106 typedef KK K;
1107 typedef vv_fes v_fes;
1108 typedef vv_fes *pfes;
1109 typedef FEbase< K, vv_fes > **R;
1110 typedef R A;
1111 typedef pfes *B;
1112 typedef CC C;
1113 typedef Expression MC;
CloneFFset3call1114 static Expression Clone(Expression e) { return e; }
CheckFFset3call1115 static bool Check(MC l) { return true; }
fFFset3call1116 static void f(KN< K > *x, MC a, A pp, Stack stack) {
1117 CC at = GetAny< C >((*a)(stack));
1118 at.call(*x);
1119 }
1120 };
1121
1122 template< typename FF >
1123 class init_FE_eqarray : public OneOperator {
1124 public:
1125 // <pferbase*,pferbase*,pfes*
1126 typedef typename FF::K K;
1127 typedef typename FF::v_fes v_fes;
1128 typedef v_fes *pfes;
1129 typedef typename v_fes::FESpace FESpace;
1130 typedef FEbase< K, v_fes > **R;
1131 typedef R A;
1132 typedef pfes *B;
1133 typedef typename FF::C C;
1134 typedef typename FF::MC MC;
1135 class CODE : public E_F0 {
1136 public:
1137 Expression a0, a1;
1138 MC a2;
CODE(Expression aa0,Expression aa1,Expression aa2)1139 CODE(Expression aa0, Expression aa1, Expression aa2) : a0(aa0), a1(aa1), a2(FF::Clone(aa2)) {
1140 ffassert(FF::Check(a2));
1141 }
operator ( )(Stack s) const1142 AnyType operator( )(Stack s) const {
1143
1144 AnyType aa0 = (*a0)(s);
1145 AnyType aa1 = (*a1)(s);
1146
1147 A pp = GetAny< A >(aa0);
1148 B a = GetAny< B >(aa1);
1149 *pp = new FEbase< K, v_fes >(a);
1150 KN< K > *x = (**pp).x( );
1151 if (!x) { // defined
1152 FESpace *Vh = (*pp)->newVh( );
1153 throwassert(Vh);
1154 **pp = x = new KN< K >(Vh->NbOfDF);
1155 *x = K( );
1156 }
1157 int n = x->N( );
1158 FF::f(x, a2, pp, s);
1159 int nn = x->N( );
1160 ffassert(n == nn);
1161 return aa0;
1162 }
nbitem() const1163 virtual size_t nbitem( ) const { return a2->nbitem( ); }
MeshIndependent() const1164 bool MeshIndependent( ) const {
1165 return a0->MeshIndependent( ) && a1->MeshIndependent( ) && a2->MeshIndependent( );
1166 } //
1167 };
init_FE_eqarray(int ppref)1168 init_FE_eqarray(int ppref)
1169 : OneOperator(map_type[typeid(R).name( )], map_type[typeid(A).name( )],
1170 map_type[typeid(B).name( )], map_type[typeid(C).name( )]) {
1171 pref = ppref;
1172 }
code(const basicAC_F0 & args) const1173 E_F0 *code(const basicAC_F0 &args) const {
1174
1175 return new CODE(args[0], args[1], map_type[typeid(C).name( )]->CastTo(args[2]));
1176 }
1177 };
1178
1179 #endif
1180