1 /****************************************************************************/
2 /* This file is part of FreeFEM. */
3 /* */
4 /* FreeFEM is free software: you can redistribute it and/or modify */
5 /* it under the terms of the GNU Lesser General Public License as */
6 /* published by the Free Software Foundation, either version 3 of */
7 /* the License, or (at your option) any later version. */
8 /* */
9 /* FreeFEM is distributed in the hope that it will be useful, */
10 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
11 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
12 /* GNU Lesser General Public License for more details. */
13 /* */
14 /* You should have received a copy of the GNU Lesser General Public License */
15 /* along with FreeFEM. If not, see <http://www.gnu.org/licenses/>. */
16 /****************************************************************************/
17 // SUMMARY : ...
18 // LICENSE : LGPLv3
19 // ORG : LJLL Universite Pierre et Marie Curie, Paris, FRANCE
20 // AUTHORS : ...
21 // E-MAIL : ...
22
23 /* clang-format off */
24 //ff-c++-LIBRARY-dep:
25 //ff-c++-cpp-dep:
26 /* clang-format on */
27
28 #include "ff++.hpp"
29 #include "AddNewFE.h"
30
31 // Attention probleme de numerotation des inconnues
32 // -------------------------------------------------
33 // dans freefem, il y a un noeud par objets sommet, arete, element.
34 // et donc la numerotation des dl dans l'element depend
35 // de l'orientation des aretes
36 //
37 /// ---------------------------------------------------------------
38 namespace Fem2D {
39 // -------------------
40 // ttdc1_ finite element fully discontinue.
41 // -------------------
42 class TypeOfFE_P1ttdc1_ : public TypeOfFE {
43 public:
44 static int Data[];
45 static double Pi_h_coef[];
46 static const R2 G;
47 static const R cshrink;
48 static const R cshrink1;
49 // (1 -1/3)*
Shrink(const R2 & P)50 static R2 Shrink(const R2 &P) { return (P - G) * cshrink + G; }
51
Shrink1(const R2 & P)52 static R2 Shrink1(const R2 &P) { return (P - G) * cshrink1 + G; }
53
TypeOfFE_P1ttdc1_()54 TypeOfFE_P1ttdc1_( ) : TypeOfFE(0, 0, 3, 1, Data, 1, 1, 3, 3, Pi_h_coef) {
55 const R2 Pt[] = {Shrink(R2(0, 0)), Shrink(R2(1, 0)), Shrink(R2(0, 1))};
56
57 for (int i = 0; i < NbDoF; i++) {
58 pij_alpha[i] = IPJ(i, i, 0);
59 P_Pi_h[i] = Pt[i];
60 }
61 }
62
63 void FB(const bool *whatd, const Mesh &Th, const Triangle &K, const RdHat &PHat,
64 RNMK_ &val) const;
65 virtual R operator( )(const FElement &K, const R2 &PHat, const KN_< R > &u, int componante,
66 int op) const;
67 };
68
69 const R2 TypeOfFE_P1ttdc1_::G(1. / 3., 1. / 3.);
70 const R TypeOfFE_P1ttdc1_::cshrink = 1;
71 const R TypeOfFE_P1ttdc1_::cshrink1 = 1. / TypeOfFE_P1ttdc1_::cshrink;
72
73 class TypeOfFE_P2ttdc1_ : public TypeOfFE {
74 public:
75 static int Data[];
76 static double Pi_h_coef[];
77 static const R2 G;
78 static const R cshrink;
79 static const R cshrink1;
Shrink(const R2 & P)80 static R2 Shrink(const R2 &P) { return (P - G) * cshrink + G; }
81
Shrink1(const R2 & P)82 static R2 Shrink1(const R2 &P) { return (P - G) * cshrink1 + G; }
83
TypeOfFE_P2ttdc1_()84 TypeOfFE_P2ttdc1_( ) : TypeOfFE(0, 0, 6, 1, Data, 3, 1, 6, 6, Pi_h_coef) {
85 const R2 Pt[] = {Shrink(R2(0, 0)), Shrink(R2(1, 0)), Shrink(R2(0, 1)),
86 Shrink(R2(0.5, 0.5)), Shrink(R2(0, 0.5)), Shrink(R2(0.5, 0))};
87
88 for (int i = 0; i < NbDoF; i++) {
89 pij_alpha[i] = IPJ(i, i, 0);
90 P_Pi_h[i] = Pt[i];
91 }
92 }
93
94 void FB(const bool *whatd, const Mesh &Th, const Triangle &K, const RdHat &PHat,
95 RNMK_ &val) const;
96 };
97
98 // on what nu df on node node of df
99 int TypeOfFE_P1ttdc1_::Data[] = {6, 6, 6, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 1, 2, 0, 0, 3};
100 int TypeOfFE_P2ttdc1_::Data[] = {6, 6, 6, 6, 6, 6, 0, 1, 2, 3, 4, 5, 0, 0, 0, 0, 0,
101 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 4, 5, 0, 0, 6};
102 double TypeOfFE_P1ttdc1_::Pi_h_coef[] = {1., 1., 1.};
103 double TypeOfFE_P2ttdc1_::Pi_h_coef[] = {1., 1., 1., 1., 1., 1.};
104 const R2 TypeOfFE_P2ttdc1_::G(1. / 3., 1. / 3.);
105 const R TypeOfFE_P2ttdc1_::cshrink = 1;
106 const R TypeOfFE_P2ttdc1_::cshrink1 = 1. / TypeOfFE_P2ttdc1_::cshrink;
operator ( )(const FElement & K,const R2 & P1Hat,const KN_<R> & u,int componante,int op) const107 R TypeOfFE_P1ttdc1_::operator( )(const FElement &K, const R2 &P1Hat, const KN_< R > &u,
108 int componante, int op) const {
109 R2 PHat = Shrink1(P1Hat);
110 R u0(u(K(0))), u1(u(K(1))), u2(u(K(2)));
111 R r = 0;
112
113 if (op == 0) {
114 R l0 = 1 - PHat.x - PHat.y, l1 = PHat.x, l2 = PHat.y;
115 r = u0 * l0 + u1 * l1 + l2 * u2;
116 } else {
117 const Triangle &T = K.T;
118 R2 D0 = T.H(0) * cshrink1, D1 = T.H(1) * cshrink1, D2 = T.H(2) * cshrink1;
119 if (op == 1) {
120 r = D0.x * u0 + D1.x * u1 + D2.x * u2;
121 } else {
122 r = D0.y * u0 + D1.y * u1 + D2.y * u2;
123 }
124 }
125
126 return r;
127 }
128
FB(const bool * whatd,const Mesh &,const Triangle & K,const RdHat & PHat,RNMK_ & val) const129 void TypeOfFE_P1ttdc1_::FB(const bool *whatd, const Mesh &, const Triangle &K, const RdHat &PHat,
130 RNMK_ &val) const {
131 R2 P = Shrink1(PHat);
132
133 // const Triangle & K(FE.T);
134 R2 A(K[0]), B(K[1]), C(K[2]);
135 R l0 = 1 - P.x - P.y, l1 = P.x, l2 = P.y;
136
137 if (val.N( ) < 3) {
138 throwassert(val.N( ) >= 3);
139 }
140
141 throwassert(val.M( ) == 1);
142
143 val = 0;
144 RN_ f0(val('.', 0, op_id));
145
146 if (whatd[op_id]) {
147 f0[0] = l0;
148 f0[1] = l1;
149 f0[2] = l2;
150 }
151
152 if (whatd[op_dx] || whatd[op_dy]) {
153 R2 Dl0(K.H(0) * cshrink1), Dl1(K.H(1) * cshrink1), Dl2(K.H(2) * cshrink1);
154
155 if (whatd[op_dx]) {
156 RN_ f0x(val('.', 0, op_dx));
157 f0x[0] = Dl0.x;
158 f0x[1] = Dl1.x;
159 f0x[2] = Dl2.x;
160 }
161
162 if (whatd[op_dy]) {
163 RN_ f0y(val('.', 0, op_dy));
164 f0y[0] = Dl0.y;
165 f0y[1] = Dl1.y;
166 f0y[2] = Dl2.y;
167 }
168 }
169 }
170
FB(const bool * whatd,const Mesh &,const Triangle & K,const RdHat & PHat,RNMK_ & val) const171 void TypeOfFE_P2ttdc1_::FB(const bool *whatd, const Mesh &, const Triangle &K, const RdHat &PHat,
172 RNMK_ &val) const {
173 R2 P = Shrink1(PHat);
174
175 // const Triangle & K(FE.T);
176 R2 A(K[0]), B(K[1]), C(K[2]);
177 R l0 = 1 - P.x - P.y, l1 = P.x, l2 = P.y;
178 R l4_0 = (4 * l0 - 1), l4_1 = (4 * l1 - 1), l4_2 = (4 * l2 - 1);
179
180 throwassert(val.N( ) >= 6);
181 throwassert(val.M( ) == 1);
182
183 val = 0;
184 if (whatd[op_id]) {
185 RN_ f0(val('.', 0, op_id));
186 f0[0] = l0 * (2 * l0 - 1);
187 f0[1] = l1 * (2 * l1 - 1);
188 f0[2] = l2 * (2 * l2 - 1);
189 f0[3] = 4 * l1 * l2; // oppose au sommet 0
190 f0[4] = 4 * l0 * l2; // oppose au sommet 1
191 f0[5] = 4 * l1 * l0; // oppose au sommet 3
192 }
193
194 if (whatd[op_dx] || whatd[op_dy] || whatd[op_dxx] || whatd[op_dyy] || whatd[op_dxy]) {
195 R2 Dl0(K.H(0) * cshrink1), Dl1(K.H(1) * cshrink1), Dl2(K.H(2) * cshrink1);
196 if (whatd[op_dx]) {
197 RN_ f0x(val('.', 0, op_dx));
198 f0x[0] = Dl0.x * l4_0;
199 f0x[1] = Dl1.x * l4_1;
200 f0x[2] = Dl2.x * l4_2;
201 f0x[3] = 4 * (Dl1.x * l2 + Dl2.x * l1);
202 f0x[4] = 4 * (Dl2.x * l0 + Dl0.x * l2);
203 f0x[5] = 4 * (Dl0.x * l1 + Dl1.x * l0);
204 }
205
206 if (whatd[op_dy]) {
207 RN_ f0y(val('.', 0, op_dy));
208 f0y[0] = Dl0.y * l4_0;
209 f0y[1] = Dl1.y * l4_1;
210 f0y[2] = Dl2.y * l4_2;
211 f0y[3] = 4 * (Dl1.y * l2 + Dl2.y * l1);
212 f0y[4] = 4 * (Dl2.y * l0 + Dl0.y * l2);
213 f0y[5] = 4 * (Dl0.y * l1 + Dl1.y * l0);
214 }
215
216 if (whatd[op_dxx]) {
217 RN_ fxx(val('.', 0, op_dxx));
218
219 fxx[0] = 4 * Dl0.x * Dl0.x;
220 fxx[1] = 4 * Dl1.x * Dl1.x;
221 fxx[2] = 4 * Dl2.x * Dl2.x;
222 fxx[3] = 8 * Dl1.x * Dl2.x;
223 fxx[4] = 8 * Dl0.x * Dl2.x;
224 fxx[5] = 8 * Dl0.x * Dl1.x;
225 }
226
227 if (whatd[op_dyy]) {
228 RN_ fyy(val('.', 0, op_dyy));
229 fyy[0] = 4 * Dl0.y * Dl0.y;
230 fyy[1] = 4 * Dl1.y * Dl1.y;
231 fyy[2] = 4 * Dl2.y * Dl2.y;
232 fyy[3] = 8 * Dl1.y * Dl2.y;
233 fyy[4] = 8 * Dl0.y * Dl2.y;
234 fyy[5] = 8 * Dl0.y * Dl1.y;
235 }
236
237 if (whatd[op_dxy]) {
238 assert(val.K( ) > op_dxy);
239 RN_ fxy(val('.', 0, op_dxy));
240
241 fxy[0] = 4 * Dl0.x * Dl0.y;
242 fxy[1] = 4 * Dl1.x * Dl1.y;
243 fxy[2] = 4 * Dl2.x * Dl2.y;
244 fxy[3] = 4 * (Dl1.x * Dl2.y + Dl1.y * Dl2.x);
245 fxy[4] = 4 * (Dl0.x * Dl2.y + Dl0.y * Dl2.x);
246 fxy[5] = 4 * (Dl0.x * Dl1.y + Dl0.y * Dl1.x);
247 }
248 }
249 }
250
251 //
252 // end ttdc1_
253
254
255
256 // 3d
257 template<class MMesh>
258 void SetPtPkDC(typename MMesh::Element::RdHat *Pt, int kk, int nn, R cc = 1) ;
259 //volume
260 template<>
SetPtPkDC(R3 * Pt,int kk,int nn,R cc)261 void SetPtPkDC<Mesh3>(R3 *Pt, int kk, int nn, R cc) { // P0 P1 et P2 , P1b
262 const int d = 3;
263 int n = 0;
264 double dK = kk;
265 double cc1 = 1 - cc; //
266 const R3 G = R3::diag(1. / (d + 1)); // barycenter
267
268 for (int i = 0; i <= kk; ++i) {
269 for (int j = 0; j <= kk - i; ++j) {
270 for (int k = 0; k <= kk - i - j; ++k) {
271 int l = kk - i - j - k;
272 ffassert(l >= 0 && l <= kk);
273 Pt[n++] = R3(k / dK, j / dK, i / dK) * cc + G * cc1;
274 }
275 }
276 }
277
278 ffassert(n == nn);
279 if (verbosity > 9) {
280 cout << " Pkdc = " << KN_< R3 >(Pt, nn) << "\n";
281 }
282 }
283
284 //surface
285 template<>
SetPtPkDC(R2 * Pt,int kk,int nn,R cc)286 void SetPtPkDC<MeshS>(R2 *Pt, int kk, int nn, R cc) {
287 const int dHat = 2;
288 int n = 0;
289 double dK = kk;
290 double cc1 = 1 - cc;
291 const R2 G = R2::diag(1. / (dHat + 1)); // barycenter
292
293 for (int i = 0; i <= kk; ++i)
294 for (int j = 0; j <= kk - i; ++j)
295 Pt[n++] = R2(j / dK, i / dK) * cc + G * cc1;
296 ffassert(n == nn);
297 if (verbosity > 9) {
298 cout << " Pkdc = " << KN_< R2 >(Pt, nn) << "\n";
299 }
300 }
301
302 //curve
303 template<>
SetPtPkDC(R1 * Pt,int kk,int nn,R cc)304 void SetPtPkDC<MeshL>(R1 *Pt, int kk, int nn, R cc) {
305 const int dHat = 2;
306 int n = 0;
307 double dK = kk;
308 double cc1 = 1 - cc;
309 const R G = R(1./2.); // barycenter
310
311 for (int i = 0; i <= kk; ++i)
312 Pt[n++] = R1(i / dK) * cc + G * cc1;
313 ffassert(n == nn);
314 if (verbosity > 9) {
315 cout << " Pkdc = " << KN_< R1 >(Pt, nn) << "\n";
316 }
317 }
318
319
320
321 // 3d
322 template<class MMesh>
323 class TypeOfFE_LagrangeDC3d : public GTypeOfFE< MMesh > {
324 // typedef typename MMesh Mesh;
325
326 public:
327 typedef GFElement<MMesh> FElement;
328 typedef typename MMesh::Element Element;
329 typedef typename Element::Rd Rd;
330 typedef typename Element::RdHat RdHat;
331 static const int d = Rd::d;
332 static const int dHat = RdHat::d;
333 const R cshrink;
334 const R cshrink1;
335 static const RdHat G;
336
Shrink(const RdHat & P) const337 RdHat Shrink(const RdHat &P) const { return (P - G) * cshrink + G; }
338
Shrink1(const RdHat & P) const339 RdHat Shrink1(const RdHat &P) const { return (P - G) * cshrink1 + G; }
340
341 const int k;
342 struct A4 {
343 int dfon[4];
344
A4Fem2D::TypeOfFE_LagrangeDC3d::A4345 A4(int k) {
346 int ndf = (dHat == 3) ? ((k + 3) * (k + 2) * (k + 1) / 6)
347 : ((dHat == 2) ? ((k + 2) * (k + 1) / 2) : k + 1);
348
349 dfon[0] = dfon[1] = dfon[2] = dfon[3] = 0;
350 dfon[dHat] = ndf;
351
352 if (verbosity > 9) {
353 cout << "A4 " << k << " " << dfon[0] << dfon[1] << dfon[2] << dfon[3] << endl;
354 }
355 }
356
operator const int*Fem2D::TypeOfFE_LagrangeDC3d::A4357 operator const int *( ) const { return dfon; }
358 };
359 RdHat *Pt;
TypeOfFE_LagrangeDC3d(int kk,R cc)360 TypeOfFE_LagrangeDC3d(int kk, R cc)
361 : // dfon ,N,nsub(graphique) , const mat interpolation , discontinuous
362 GTypeOfFE< MMesh >(A4(kk), 1, Max(kk, 1), true, true), cshrink(cc), cshrink1(1. / cc),
363 k(kk) {
364 int n = this->NbDoF;
365
366 if (verbosity > 9) {
367 cout << "\n +++ Pdc" << k << " : ndof : " << n << endl;
368 }
369
370 SetPtPkDC<MMesh>(this->PtInterpolation, k, this->NbDoF, cc);
371 if (verbosity > 9)
372 cout << " ppppp " << this->PtInterpolation << endl;
373
374
375 {
376 for (int i = 0; i < n; i++) {
377 this->pInterpolation[i] = i;
378 this->cInterpolation[i] = 0;
379 this->dofInterpolation[i] = i;
380 this->coefInterpolation[i] = 1.;
381 }
382 }
383 }
384
~TypeOfFE_LagrangeDC3d()385 ~TypeOfFE_LagrangeDC3d( ) {}
386
387
388 void FB(const What_d whatd, const MMesh &Th, const Element &K, const RdHat &PHat,RNMK_ &val) const;
389 R operator( )(const FElement &K, const RdHat &PHat, const KN_< R > &u, int componante,int op) const;
390
391 private:
392 TypeOfFE_LagrangeDC3d(const TypeOfFE_LagrangeDC3d &);
393 void operator=(const TypeOfFE_LagrangeDC3d &);
394 };
395
396
397
398 template<>
FB(const What_d whatd,const Mesh3 & Th,const Mesh3::Element & K,const RdHat & PHat,RNMK_ & val) const399 void TypeOfFE_LagrangeDC3d<Mesh3>::FB(const What_d whatd, const Mesh3 &Th, const Mesh3::Element &K,const RdHat &PHat, RNMK_ &val) const {
400
401 R3 P = this->Shrink1(PHat);
402 R l[] = {1. - P.sum( ), P.x, P.y, P.z};
403
404 assert(val.N( ) >= Element::nv);
405 assert(val.M( ) == 1);
406
407 val = 0;
408 RN_ f0(val('.', 0, op_id));
409
410 if (whatd & Fop_D0) {
411 f0[0] = l[0];
412 f0[1] = l[1];
413 f0[2] = l[2];
414 f0[3] = l[3];
415 }
416
417 if (whatd & Fop_D1) {
418 R3 Dl[4];
419 K.Gradlambda(Dl);
420
421 for (int i = 0; i < 4; ++i) {
422 Dl[i] *= cshrink1;
423 }
424
425 if (whatd & Fop_dx) {
426 RN_ f0x(val('.', 0, op_dx));
427 f0x[0] = Dl[0].x;
428 f0x[1] = Dl[1].x;
429 f0x[2] = Dl[2].x;
430 f0x[3] = Dl[3].x;
431 }
432
433 if (whatd & Fop_dy) {
434 RN_ f0y(val('.', 0, op_dy));
435 f0y[0] = Dl[0].y;
436 f0y[1] = Dl[1].y;
437 f0y[2] = Dl[2].y;
438 f0y[3] = Dl[3].y;
439 }
440
441 if (whatd & Fop_dz) {
442 RN_ f0z(val('.', 0, op_dz));
443 f0z[0] = Dl[0].z;
444 f0z[1] = Dl[1].z;
445 f0z[2] = Dl[2].z;
446 f0z[3] = Dl[3].z;
447 }
448 }
449 }
450
451
452
453 template<>
operator ( )(const FElement & K,const R3 & PHat1,const KN_<R> & u,int componante,int op) const454 R TypeOfFE_LagrangeDC3d<Mesh3>::operator( )(const FElement &K, const R3 &PHat1, const KN_< R > &u,int componante, int op) const {
455 R3 PHat = Shrink1(PHat1);
456 R r = 0;
457
458 if (k == 1) {
459 R u0(u(K(0))), u1(u(K(1))), u2(u(K(2))), u3(u(K(3)));
460 if (op == 0) {
461 R l[4];
462 PHat.toBary(l);
463 r = u0 * l[0] + u1 * l[1] + l[2] * u2 + l[3] * u3;
464 } else if (op == op_dx || op == op_dy || op == op_dz) {
465 const Element &T = K.T;
466 R3 D[4];
467 T.Gradlambda(D);
468
469 for (int i = 0; i < 4; ++i) {
470 D[i] *= cshrink1;
471 }
472
473 if (op == op_dx) {
474 r = D[0].x * u0 + D[1].x * u1 + D[2].x * u2 + D[3].x * u3;
475 } else if (op == op_dy) {
476 r = D[0].y * u0 + D[1].y * u1 + D[2].y * u2 + D[3].y * u3;
477 } else {
478 r = D[0].z * u0 + D[1].z * u1 + D[2].z * u2 + D[3].z * u3;
479 }
480 }
481 } else {
482 ffassert(0); // to do ..
483 }
484
485 return r;
486 }
487
488 template<>
489 const R3 TypeOfFE_LagrangeDC3d<Mesh3>::G(1. / 4., 1. / 4., 1. / 4.);
490
491
492
493 template<>
FB(const What_d whatd,const MeshS & Th,const MeshS::Element & K,const RdHat & PHat,RNMK_ & val) const494 void TypeOfFE_LagrangeDC3d<MeshS>::FB(const What_d whatd, const MeshS &Th, const MeshS::Element &K,const RdHat &PHat, RNMK_ &val) const {
495 R2 P = this->Shrink1(PHat);
496 R l[] = {1. - P.sum( ), P.x, P.y};
497
498 assert(val.N( ) >= Element::nv);
499 assert(val.M( ) == 1);
500
501 val = 0;
502 RN_ f0(val('.', 0, op_id));
503
504 if (whatd & Fop_D0) {
505 f0[0] = l[0];
506 f0[1] = l[1];
507 f0[2] = l[2];
508 }
509
510 if (whatd & Fop_D1) {
511 R3 Dl[3];
512 K.Gradlambda(Dl);
513
514 for (int i = 0; i < 3; ++i) {
515 Dl[i] *= cshrink1;
516 }
517
518 if (whatd & Fop_dx) {
519 RN_ f0x(val('.', 0, op_dx));
520 f0x[0] = Dl[0].x;
521 f0x[1] = Dl[1].x;
522 f0x[2] = Dl[2].x;
523 }
524
525 if (whatd & Fop_dy) {
526 RN_ f0y(val('.', 0, op_dy));
527 f0y[0] = Dl[0].y;
528 f0y[1] = Dl[1].y;
529 f0y[2] = Dl[2].y;
530 }
531
532 if (whatd & Fop_dz) {
533 RN_ f0z(val('.', 0, op_dz));
534 f0z[0] = Dl[0].z;
535 f0z[1] = Dl[1].z;
536 f0z[2] = Dl[2].z;
537 }
538 }
539 }
540
541 template<>
operator ( )(const FElement & K,const R2 & PHat1,const KN_<R> & u,int componante,int op) const542 R TypeOfFE_LagrangeDC3d<MeshS>::operator( )(const FElement &K, const R2 &PHat1, const KN_< R > &u,int componante, int op) const {
543 R2 PHat = Shrink1(PHat1);
544 R r = 0;
545
546 if (k == 1) {
547 R u0(u(K(0))), u1(u(K(1))), u2(u(K(2)));
548
549 if (op == 0) {
550 R l[3];
551 PHat.toBary(l);
552 r = u0 * l[0] + u1 * l[1] + l[2] * u2 ;
553 } else if (op == op_dx || op == op_dy || op == op_dz) {
554 const Element &T = K.T;
555 R3 D[3];
556 T.Gradlambda(D);
557
558 for (int i = 0; i < 3; ++i) {
559 D[i] *= cshrink1;
560 }
561
562 if (op == op_dx) {
563 r = D[0].x * u0 + D[1].x * u1 + D[2].x * u2 ;
564 } else if (op == op_dy) {
565 r = D[0].y * u0 + D[1].y * u1 + D[2].y * u2 ;
566 } else {
567 r = D[0].z * u0 + D[1].z * u1 + D[2].z * u2 ;
568 }
569 }
570 } else {
571 ffassert(0); // to do ..
572 }
573
574 return r;
575 }
576 template<>
577 const R2 TypeOfFE_LagrangeDC3d<MeshS>::G(1. / 3., 1. / 3.);
578
579
580 template<>
FB(const What_d whatd,const MeshL & Th,const MeshL::Element & K,const RdHat & PHat,RNMK_ & val) const581 void TypeOfFE_LagrangeDC3d<MeshL>::FB(const What_d whatd, const MeshL &Th, const MeshL::Element &K,const RdHat &PHat, RNMK_ &val) const {
582 R1 P = this->Shrink1(PHat);
583 R l[] = {1. - P.x, P.x};
584
585 assert(val.N( ) >= Element::nv);
586 assert(val.M( ) == 1);
587
588 val = 0;
589 RN_ f0(val('.', 0, op_id));
590
591 if (whatd & Fop_D0) {
592 f0[0] = l[0];
593 f0[1] = l[1];
594 }
595
596 if (whatd & Fop_D1) {
597 R3 Dl[2];
598 K.Gradlambda(Dl);
599
600 for (int i = 0; i < 2; ++i) {
601 Dl[i] *= cshrink1;
602 }
603
604 if (whatd & Fop_dx) {
605 RN_ f0x(val('.', 0, op_dx));
606 f0x[0] = Dl[0].x;
607 f0x[1] = Dl[1].x;
608 }
609
610 if (whatd & Fop_dy) {
611 RN_ f0y(val('.', 0, op_dy));
612 f0y[0] = Dl[0].y;
613 f0y[1] = Dl[1].y;
614 }
615
616 if (whatd & Fop_dz) {
617 RN_ f0z(val('.', 0, op_dz));
618 f0z[0] = Dl[0].z;
619 f0z[1] = Dl[1].z;
620 }
621 }
622 }
623
624 template<>
operator ( )(const FElement & K,const R1 & PHat1,const KN_<R> & u,int componante,int op) const625 R TypeOfFE_LagrangeDC3d<MeshL>::operator( )(const FElement &K, const R1 &PHat1, const KN_< R > &u,int componante, int op) const {
626 R1 PHat = Shrink1(PHat1);
627 R r = 0;
628
629 if (k == 1) {
630 R u0(u(K(0))), u1(u(K(1)));
631
632 if (op == 0) {
633 R l[2];
634 PHat.toBary(l);
635 r = u0 * l[0] + u1 * l[1] ;
636 } else if (op == op_dx || op == op_dy || op == op_dz) {
637 const Element &T = K.T;
638 R3 D[2];
639 T.Gradlambda(D);
640
641 for (int i = 0; i < 2; ++i) {
642 D[i] *= cshrink1;
643 }
644
645 if (op == op_dx) {
646 r = D[0].x * u0 + D[1].x * u1 ;
647 } else if (op == op_dy) {
648 r = D[0].y * u0 + D[1].y * u1 ;
649 } else {
650 r = D[0].z * u0 + D[1].z * u1 ;
651 }
652 }
653 } else {
654 ffassert(0); // to do ..
655 }
656
657 return r;
658 }
659
660 template<>
661 const R1 TypeOfFE_LagrangeDC3d<MeshL>::G(1. / 2.);
662
663
664 } // namespace Fem2D
665
finit()666 static void finit( ) {
667 // link with FreeFem++
668 static TypeOfFE_P1ttdc1_ P1dc1LagrangeP1dc1;
669 static TypeOfFE_P2ttdc1_ P2dc1LagrangeP2dc1;
670
671 static TypeOfFE_LagrangeDC3d<Mesh3> TypeOfFE_LagrangeDC3dtt(1, 0.999);
672 static TypeOfFE_LagrangeDC3d<Mesh3> TypeOfFE_LagrangeDC3dtt1(1, 1.);
673
674 static TypeOfFE_LagrangeDC3d<MeshS> TypeOfFE_LagrangeDC3dStt(1, 0.999);
675 static TypeOfFE_LagrangeDC3d<MeshS> TypeOfFE_LagrangeDC3dStt1(1, 1.);
676
677 static TypeOfFE_LagrangeDC3d<MeshL> TypeOfFE_LagrangeDC3dLtt(1, 0.999);
678 static TypeOfFE_LagrangeDC3d<MeshL> TypeOfFE_LagrangeDC3dLtt1(1, 1.);
679
680 // a static variable to add the finite element to freefem++
681 static AddNewFE P1dcLagrange("P1dc1", &P1dc1LagrangeP1dc1);
682 static AddNewFE P2dcLagrange("P2dc1", &P2dc1LagrangeP2dc1);
683
684 static AddNewFE3 P1dttLagrange3d("P1dc3d", &TypeOfFE_LagrangeDC3dtt, "P1dc");
685 static AddNewFE3 P1dttLagrange3d1("P1dc3d1", &TypeOfFE_LagrangeDC3dtt1);
686
687 static AddNewFES P1dttLagrange3dS("P1dc3dS", &TypeOfFE_LagrangeDC3dStt, "P1dc");
688 static AddNewFES P1dttLagrange3dS1("P1dc3dS1", &TypeOfFE_LagrangeDC3dStt1);
689
690 static AddNewFEL P1dttLagrange3dL("P1dc3dL", &TypeOfFE_LagrangeDC3dLtt, "P1dc");
691 static AddNewFEL P1dttLagrange3dL1("P1dc3dL1", &TypeOfFE_LagrangeDC3dLtt1);
692
693 }
694
695 LOADFUNC(finit)
696
697 // --- fin --
698