1 // -*- Mode : c++ -*-
2 //
3 // SUMMARY :
4 // USAGE :
5 // ORG :
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 #include <cmath>
29 #include <cstdlib>
30 #include "error.hpp"
31 #include <iostream>
32 #include <fstream>
33 #include <map>
34 #include "rgraph.hpp"
35 using namespace std;
36
37 #include "RNM.hpp"
38 #include "fem.hpp"
39 #include "FESpace.hpp"
40
41 namespace Fem2D {
42
43 class TypeOfFE_RT : public TypeOfFE { public:
44 static int Data[];
TypeOfFE_RT()45 TypeOfFE_RT(): TypeOfFE(0,1,0,2,Data,1,1,6,3)
46 {const R2 Pt[] = { R2(0.5,0.5), R2(0.0,0.5), R2(0.5,0.0) };
47 for (int p=0,kk=0;p<3;p++)
48 { P_Pi_h[p]=Pt[p];
49 for (int j=0;j<2;j++)
50 pij_alpha[kk++]= IPJ(p,p,j);
51 }}
52 // void FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
53 void FB(const bool * watdd, const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
54 // void D2_FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
55 // void Pi_h(const baseFElement & K,RN_ & val, InterpolFunction f, R* v,int, void *) const;
56 void Pi_h_alpha(const baseFElement & K,KN_<double> & v) const ;
57 } ;
58 // on what nu df on node node of df
59 int TypeOfFE_RT::Data[]={3,4,5, 0,0,0, 0,1,2, 0,0,0, 0,1,2, 0,0, 0,0,3,3};
60
61 /* void TypeOfFE_RT::D2_FB(const Mesh & ,const Triangle & ,const R2 & ,RNMK_ & val) const
62 { //
63 val=0;
64 }*/
65 /*
66 void TypeOfFE_RT::FB(const Mesh & Th,const Triangle & K,const R2 & PHat,RNMK_ & val) const
67 { //
68 // const Triangle & K(FE.T);
69 R2 P(K(PHat));
70 R2 A(K[0]), B(K[1]),C(K[2]);
71 R l0=1-P.x-P.y,l1=P.x,l2=P.y;
72 // R2 Dl0(K.H(0)), Dl1(K.H(1)), Dl2(K.H(2));
73 if (val.N() <3)
74 throwassert(val.N() >=3);
75 throwassert(val.M()==2 );
76 throwassert(val.K()==3 );
77 RN_ f0(val('.',0,0));
78 RN_ f1(val('.',1,0));
79 val=0;
80 // RN_ df0(val(0,'.',0));
81 // RN_ fy(val('.','.',2));
82 // a_i ([x,y]-c_i) , ou c_i = A,B , C si i= 0,1,2
83 // int_T a_i div([x,y]-c_i) = 1
84 // div div([x,y]-c_i) = 2
85 // donc a_i = 1/(2 area T)
86
87 R a=1./(2*K.area);
88 R a0= K.EdgeOrientation(0) * a ;
89 R a1= K.EdgeOrientation(1) * a ;
90 R a2= K.EdgeOrientation(2) * a ;
91 // if (Th(K)< 2) cout << Th(K) << " " << A << " " << B << " " << C << "; " << a0 << " " << a1 << " "<< a2 << endl;;
92
93 // ------------
94 f0[0] = (P.x-A.x)*a0;
95 f1[0] = (P.y-A.y)*a0;
96
97 f0[1] = (P.x-B.x)*a1;
98 f1[1] = (P.y-B.y)*a1;
99
100 f0[2] = (P.x-C.x)*a2;
101 f1[2] = (P.y-C.y)*a2;
102 // ----------------
103 // ----------------
104 // BUG dans RT correct FH le 17 sept 2002
105 // dx [x,y] = [1,0] et non [1,1]
106 // dy [x,y] = [0,1] et non [1,1]
107 // -------------------------------------
108
109 val(0,0,1) = a0;
110 val(1,0,1) = a1;
111 val(2,0,1) = a2;
112 val(0,1,2) = a0;
113 val(1,1,2) = a1;
114 val(2,1,2) = a2;
115
116 }
117 */
FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 & PHat,RNMK_ & val) const118 void TypeOfFE_RT::FB(const bool *whatd,const Mesh & Th,const Triangle & K,const R2 & PHat,RNMK_ & val) const
119 { //
120 // const Triangle & K(FE.T);
121 R2 P(K(PHat));
122 R2 A(K[0]), B(K[1]),C(K[2]);
123 // R l0=1-P.x-P.y,l1=P.x,l2=P.y;
124 // R2 Dl0(K.H(0)), Dl1(K.H(1)), Dl2(K.H(2));
125 if (val.N() <3)
126 throwassert(val.N() >=3);
127 throwassert(val.M()==2 );
128 // throwassert(val.K()==3 );
129 val=0;
130 R a=1./(2*K.area);
131 R a0= K.EdgeOrientation(0) * a ;
132 R a1= K.EdgeOrientation(1) * a ;
133 R a2= K.EdgeOrientation(2) * a ;
134 // if (Th(K)< 2) cout << Th(K) << " " << A << " " << B << " " << C << "; " << a0 << " " << a1 << " "<< a2 << endl;;
135
136 // ------------
137 if (whatd[op_id])
138 {
139 assert(val.K()>op_id);
140 RN_ f0(val('.',0,0));
141 RN_ f1(val('.',1,0));
142 f0[0] = (P.x-A.x)*a0;
143 f1[0] = (P.y-A.y)*a0;
144
145 f0[1] = (P.x-B.x)*a1;
146 f1[1] = (P.y-B.y)*a1;
147
148 f0[2] = (P.x-C.x)*a2;
149 f1[2] = (P.y-C.y)*a2;
150 }
151 // ----------------
152 // BUG dans RT correct FH le 17 sept 2002
153 // dx [x,y] = [1,0] et non [1,1]
154 // dy [x,y] = [0,1] et non [1,1]
155 // -------------------------------------
156 if (whatd[op_dx])
157 {
158 assert(val.K()>op_dx);
159 val(0,0,op_dx) = a0;
160 val(1,0,op_dx) = a1;
161 val(2,0,op_dx) = a2;
162 }
163 if (whatd[op_dy])
164 {
165 assert(val.K()>op_dy);
166 val(0,1,op_dy) = a0;
167 val(1,1,op_dy) = a1;
168 val(2,1,op_dy) = a2;
169 }
170 }
171 /*
172 void TypeOfFE_RT::Pi_h(const baseFElement & K,RN_ & val, InterpolFunction f, R* v,int j, void * arg) const
173 {
174 const R2 Pt[] = { R2(0.5,0.5), R2(0.0,0.5), R2(0.5,0.0) };
175 const Triangle & T(K.T);
176 // if (K.number<2) cout << K.number << ": " ;
177 for (int i=0;i<3;i++)
178 {
179 f(v,T(Pt[i]),K,i,Pt[i],arg);
180 R2 E(T.Edge(i));
181 R signe = T.EdgeOrientation(i) ;
182 val[i]= signe*(v[j]*E.y-v[j+1]*E.x);
183 // if (K.number<2) cout << val[i] << " " ;
184 }
185 // if (K.number<2) cout << endl;
186 }
187 */
Pi_h_alpha(const baseFElement & K,KN_<double> & v) const188 void TypeOfFE_RT::Pi_h_alpha(const baseFElement & K,KN_<double> & v) const
189 {
190 const Triangle & T(K.T);
191
192 for (int i=0,k=0;i<3;i++)
193 {
194 R2 E(T.Edge(i));
195 R signe = T.EdgeOrientation(i) ;
196 v[k++]= signe*E.y;
197 v[k++]=-signe*E.x;
198 }
199 }
200 // -------------------
201
202
203
204 class TypeOfFE_RTmodif : public TypeOfFE { public:
205 static int Data[];
TypeOfFE_RTmodif()206 TypeOfFE_RTmodif(): TypeOfFE(0,1,0,2,Data,1,1,6,3)
207 {const R2 Pt[] = { R2(0.5,0.5), R2(0.0,0.5), R2(0.5,0.0) };
208 for (int p=0,kk=0;p<3;p++)
209 { P_Pi_h[p]=Pt[p];
210 for (int j=0;j<2;j++)
211 pij_alpha[kk++]= IPJ(p,p,j);
212 }}
213 // void FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
214 void FB(const bool * whatd, const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
215 // void D2_FB(const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
216 // void Pi_h(const baseFElement & K,RN_ & val, InterpolFunction f, R* v,int, void *) const;
217 void Pi_h_alpha(const baseFElement & K,KN_<double> & v) const ;
218 } ;
219 // on what nu df on node node of df
220 int TypeOfFE_RTmodif::Data[]={3,4,5, 0,0,0, 0,1,2, 0,0,0, 0,1,2, 0,0, 0,0, 3,3};
221
222
FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 & PHat,RNMK_ & val) const223 void TypeOfFE_RTmodif::FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 & PHat,RNMK_ & val) const
224 { //
225 // const Triangle & K(FE.T);
226 R2 P(K(PHat));
227 R2 A(K[0]), B(K[1]),C(K[2]);
228 R la=1-PHat.x-PHat.y,lb=PHat.x,lc=PHat.y;
229 R2 Dla(K.H(0)), Dlb(K.H(1)), Dlc(K.H(2));
230 if (val.N() <3)
231 throwassert(val.N() >=3);
232 throwassert(val.M()==2 );
233
234 R2 AB(A,B),AC(A,C),BA(B,A),BC(B,C),CA(C,A),CB(C,B);
235
236 R aa0= 1./(((AB,Dlb) + (AC,Dlc))*K.area);
237 R aa1= 1./(((BA,Dla) + (BC,Dlc))*K.area);
238 R aa2= 1./(((CA,Dla) + (CB,Dlb))*K.area);
239 int i=0;
240 R a0= &K[ (i+1)%3] < &K[ (i+2)%3] ? aa0 : -aa0 ;
241 i=1;
242 R a1= &K[ (i+1)%3] < &K[ (i+2)%3] ? aa1 : -aa1 ;
243 i=2;
244 R a2= &K[ (i+1)%3] < &K[ (i+2)%3] ? aa2 : -aa2 ;
245 // if (Th(K)< 2) cout << Th(K) << " " << A << " " << B << " " << C << "; " << a0 << " " << a1 << " "<< a2 << endl;;
246
247 R2 Va= AB*(lb*a0) + AC*(lc*a0);
248 R2 Vb= BA*(la*a1) + BC*(lc*a1);
249 R2 Vc= CA*(la*a2) + CB*(lb*a2);
250 R2 Va_x= AB*(Dlb.x*a0) + AC*(Dlc.x*a0);
251 R2 Vb_x= BA*(Dla.x*a1) + BC*(Dlc.x*a1);
252 R2 Vc_x= CA*(Dla.x*a2) + CB*(Dlb.x*a2);
253 R2 Va_y= AB*(Dlb.y*a0) + AC*(Dlc.y*a0);
254 R2 Vb_y= BA*(Dla.y*a1) + BC*(Dlc.y*a1);
255 R2 Vc_y= CA*(Dla.y*a2) + CB*(Dlb.y*a2);
256
257 if( whatd[op_id])
258 {
259 RN_ f0(val('.',0,0));
260 RN_ f1(val('.',1,0));
261
262 f0[0] = Va.x;
263 f1[0] = Va.y;
264
265 f0[1] = Vb.x;
266 f1[1] = Vb.y;
267
268 f0[2] = Vc.x;
269 f1[2] = Vc.y;
270 }
271 // ----------------
272 if( whatd[op_dx])
273 {
274 val(0,0,1) = Va_x.x;
275 val(0,1,1) = Va_x.y;
276
277 val(1,0,1) = Vb_x.x;
278 val(1,1,1) = Vb_x.y;
279
280 val(2,0,1) = Vc_x.x;
281 val(2,1,1) = Vc_x.y;
282 }
283
284 if( whatd[op_dy])
285 {
286 val(0,0,2) = Va_y.x;
287 val(0,1,2) = Va_y.y;
288
289 val(1,0,2) = Vb_y.x;
290 val(1,1,2) = Vb_y.y;
291
292 val(2,0,2) = Vc_y.x;
293 val(2,1,2) = Vc_y.y;
294 }
295
296 }
297
Pi_h_alpha(const baseFElement & K,KN_<double> & v) const298 void TypeOfFE_RTmodif::Pi_h_alpha(const baseFElement & K,KN_<double> & v) const
299 {
300 const Triangle & T(K.T);
301
302 for (int i=0,k=0;i<3;i++)
303 {
304 R2 E(T.Edge(i));
305 R signe = &T[ (i+1)%3] < &T[ (i+2)%3] ? 1.0 : -1.0 ;
306 v[k++]= signe*E.y;
307 v[k++]=-signe*E.x;
308 }
309 }
310
311
312 // ---------------------
313 class TypeOfFE_P0 : public TypeOfFE { public:
314 static int Data[];
315 static double Pi_h_coef[];
316
TypeOfFE_P0()317 TypeOfFE_P0(): TypeOfFE(0,0,1,1,Data,1,1,1,1,Pi_h_coef){
318 pij_alpha[0]=IPJ(0,0,0);
319 P_Pi_h[0]=R2(1./3.,1./3.);
320 }
321 void FB(const bool * watdd, const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
322
323 };
324
325 // on what nu df on node node of df
326 int TypeOfFE_P0::Data[]={6, 0, 0, 0 , 0 ,0, 0, 1 };
327 double TypeOfFE_P0::Pi_h_coef[]={1.0};
328
329
FB(const bool * whatd,const Mesh &,const Triangle & K,const R2 & PHat,RNMK_ & val) const330 void TypeOfFE_P0::FB(const bool* whatd,const Mesh & ,const Triangle & K,const R2 & PHat,RNMK_ & val) const
331 { //
332 // const Triangle & K(FE.T);
333 R2 P(K(PHat));
334 R2 A(K[0]), B(K[1]),C(K[2]);
335 //R l0=1-P.x-P.y,l1=P.x,l2=P.y;
336 R2 Dl0(K.H(0)), Dl1(K.H(1)), Dl2(K.H(2));
337 throwassert(val.N() >=1);
338 throwassert(val.M()==1 );
339 // throwassert(val.K()==3 );
340 val=0;
341 if ( whatd[op_id])
342 val(0,0,0) =1;
343 }
344
345 // ------ P1 non conforme --------
346 class TypeOfFE_P1ncLagrange : public TypeOfFE { public:
347 static int Data[];
348 static double Pi_h_coef[];
TypeOfFE_P1ncLagrange()349 TypeOfFE_P1ncLagrange(): TypeOfFE(0,1,0,1,Data,1,1,3,3,Pi_h_coef)
350 {
351 const R2 Pt[] = { R2(0.5,0.5), R2(0.0,0.5), R2(0.5,0.0) };
352 for (int i=0;i<NbDoF;i++) {
353 pij_alpha[i]= IPJ(i,i,0);
354 P_Pi_h[i]=Pt[i]; }
355 }
356
357 void FB(const bool * whatd, const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
358
359 } ;
360
361 // on what nu df on node node of df
362 int TypeOfFE_P1ncLagrange::Data[]={3,4,5, 0,0,0, 0,1,2, 0,0,0, 0,1,2, 0, 0,3};
363 double TypeOfFE_P1ncLagrange::Pi_h_coef[]={1.,1.,1.};
364
365
366 class TypeOfFE_ConsEdge : public TypeOfFE { public:
367 static int Data[];
368 static double Pi_h_coef[];
369
TypeOfFE_ConsEdge()370 TypeOfFE_ConsEdge(): TypeOfFE(0,1,0,1,Data,3,1,3,3,Pi_h_coef)
371 { const R2 Pt[] = { R2(0.5,0.5), R2(0.0,0.5), R2(0.5,0.0) };
372 for (int i=0;i<NbDoF;i++) {
373 pij_alpha[i]= IPJ(i,i,0);
374 P_Pi_h[i]=Pt[i]; }
375 }
376
377 void FB(const bool * whatd, const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
378
379 } ;
380 // on what nu df on node node of df
381 int TypeOfFE_ConsEdge::Data[]={3,4,5, 0,0,0, 0,1,2, 0,0,0, 0,1,2, 0, 0,3};
382 double TypeOfFE_ConsEdge::Pi_h_coef[]={1.,1.,1.};
FB(const bool * whatd,const Mesh &,const Triangle & K,const R2 & P,RNMK_ & val) const383 void TypeOfFE_ConsEdge::FB(const bool * whatd,const Mesh & ,const Triangle & K,const R2 & P,RNMK_ & val) const
384 {
385 // const Triangle & K(FE.T);
386 R2 A(K[0]), B(K[1]),C(K[2]);
387 R l0=1-P.x-P.y,l1=P.x,l2=P.y;
388
389 if (val.N() <3)
390 throwassert(val.N() >=3);
391 throwassert(val.M()==1 );
392
393 val=0;
394 if (whatd[op_id])
395 {
396
397 RN_ f0(val('.',0,0));
398 //
399 f0[0] = double(l0 <= min(l1,l2) ); // arete
400 f0[1] = double(l1 <= min(l0,l2) );
401 f0[2] = double(l2 <= min(l0,l1) );
402 }
403
404 }
405 /*
406 class TypeOfFE_P1Edge : public TypeOfFE { public:
407 static int Data[];
408 static double Pi_h_coef[];
409
410 TypeOfFE_P1Edge(): TypeOfFE(0,2,0,1,Data,3,1,12,6,Pi_h_coef)
411 { R2 Pt[6] ;
412
413 int kk=0;
414 for(int i=0;i<3;++i)
415 for(int j=0;i<QF_GaussLegendre2.n;++j)
416 { R2 A(TriangleHat[VerticesOfTriangularEdge[i][0]]);
417 R2 B(TriangleHat[VerticesOfTriangularEdge[i][1]]);
418 Pt[k++]=A*(QF_GaussLegendre2[j].x)+ B*(1-QF_GaussLegendre2[j].x)
419 }
420
421 int other[6]= { 1,0, 3,2,5,4 };
422 k=0;
423 for (int i=0;i<NbDoF;i++) {
424 pij_alpha[i]= IPJ(i,i,0);
425 if(other[i]>=0)
426 pij_alpha[kk++]= IPJ(i,other[i],0);
427 P_Pi_h[i]=Pt[i]; }
428 }
429
430 void FB(const bool * whatd, const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
431
432 } ;
433 // on what nu df on node node of df
434 int TypeOfFE_ConsEdge::Data[]={3,4,5, 0,0,0, 0,1,2, 0,0,0, 0,1,2, 0, 0,3};
435 double TypeOfFE_ConsEdge::Pi_h_coef[]={1.,1.,1.};
436 void TypeOfFE_ConsEdge::FB(const bool * whatd,const Mesh & ,const Triangle & K,const R2 & P,RNMK_ & val) const
437 {
438 // const Triangle & K(FE.T);
439 R2 A(K[0]), B(K[1]),C(K[2]);
440 R l0=1-P.x-P.y,l1=P.x,l2=P.y;
441
442 if (val.N() <3)
443 throwassert(val.N() >=3);
444 throwassert(val.M()==1 );
445
446 val=0;
447 if (whatd[op_id])
448 {
449
450 RN_ f0(val('.',0,0));
451 //
452 f0[0] = double(l0 <= min(l1,l2) ); // arete
453 f0[1] = double(l1 <= min(l0,l2) );
454 f0[2] = double(l2 <= min(l0,l1) );
455 }
456
457 }
458
459 */
460
FB(const bool * whatd,const Mesh &,const Triangle & K,const R2 & P,RNMK_ & val) const461 void TypeOfFE_P1ncLagrange::FB(const bool * whatd,const Mesh & ,const Triangle & K,const R2 & P,RNMK_ & val) const
462 {
463 // const Triangle & K(FE.T);
464 R2 A(K[0]), B(K[1]),C(K[2]);
465 // l1( cshrink1*(cshrink*((1,0)-G)+G)-G)+G = 1
466 R l0=1-P.x-P.y,l1=P.x,l2=P.y;
467
468 if (val.N() <3)
469 throwassert(val.N() >=3);
470 throwassert(val.M()==1 );
471 // throwassert(val.K()==3 );
472
473 val=0;
474 if (whatd[op_id])
475 {
476 RN_ f0(val('.',0,0));
477 f0[0] = 1-l0*2;
478 f0[1] = 1-l1*2;
479 f0[2] = 1-l2*2;
480 }
481 if (whatd[op_dx] || whatd[op_dy] )
482 {
483 R2 Dl0(K.H(0)), Dl1(K.H(1)), Dl2(K.H(2));
484 if (whatd[op_dx])
485 {
486 RN_ f0x(val('.',0,op_dx));
487 f0x[0] = -Dl0.x*2;
488 f0x[1] = -Dl1.x*2;
489 f0x[2] = -Dl2.x*2;
490 }
491 if (whatd[op_dy])
492 {
493 RN_ f0y(val('.',0,op_dy));
494 f0y[0] = -Dl0.y*2;
495 f0y[1] = -Dl1.y*2;
496 f0y[2] = -Dl2.y*2;
497 }
498 }
499 }
500
501 // The RT orthogonal FE
502
503 class TypeOfFE_RTortho : public TypeOfFE { public:
504 static int Data[];
TypeOfFE_RTortho()505 TypeOfFE_RTortho(): TypeOfFE(0,1,0,2,Data,1,1,6,3)
506 {const R2 Pt[] = { R2(0.5,0.5), R2(0.0,0.5), R2(0.5,0.0) };
507 for (int p=0,kk=0;p<3;p++)
508 { P_Pi_h[p]=Pt[p];
509 for (int j=0;j<2;j++)
510 pij_alpha[kk++]= IPJ(p,p,j);
511 }}
512
513 void FB(const bool * watdd, const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
514 void Pi_h_alpha(const baseFElement & K,KN_<double> & v) const ;
515 } ;
516 // on what nu df on node node of df
517 int TypeOfFE_RTortho::Data[]={3,4,5, 0,0,0, 0,1,2, 0,0,0, 0,1,2, 0,0, 0,0, 3,3};
518
519
520
FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 & PHat,RNMK_ & val) const521 void TypeOfFE_RTortho::FB(const bool *whatd,const Mesh & Th,const Triangle & K,const R2 & PHat,RNMK_ & val) const
522 { //
523 // const Triangle & K(FE.T);
524 R2 P(K(PHat));
525 R2 A(K[0]), B(K[1]),C(K[2]);
526 //R l0=1-P.x-P.y,l1=P.x,l2=P.y;
527 // R2 Dl0(K.H(0)), Dl1(K.H(1)), Dl2(K.H(2));
528 if (val.N() <3)
529 throwassert(val.N() >=3);
530 throwassert(val.M()==2 );
531 // throwassert(val.K()==3 );
532 val=0;
533 R a=1./(2*K.area);
534 R a0= K.EdgeOrientation(0) * a ;
535 R a1= K.EdgeOrientation(1) * a ;
536 R a2= K.EdgeOrientation(2) * a ;
537 // if (Th(K)< 2) cout << Th(K) << " " << A << " " << B << " " << C << "; " << a0 << " " << a1 << " "<< a2 << endl;;
538
539 // ------------
540 if (whatd[op_id])
541 {
542 assert(val.K()>op_id);
543 RN_ f0(val('.',0,0));
544 RN_ f1(val('.',1,0));
545 f1[0] = (P.x-A.x)*a0;
546 f0[0] = -(P.y-A.y)*a0;
547
548 f1[1] = (P.x-B.x)*a1;
549 f0[1] = -(P.y-B.y)*a1;
550
551 f1[2] = (P.x-C.x)*a2;
552 f0[2] = -(P.y-C.y)*a2;
553 }
554 // ----------------
555 if (whatd[op_dx])
556 {
557 assert(val.K()>op_dx);
558 val(0,1,op_dx) = a0;
559 val(1,1,op_dx) = a1;
560 val(2,1,op_dx) = a2;
561 }
562 if (whatd[op_dy])
563 {
564 assert(val.K()>op_dy);
565 val(0,0,op_dy) = -a0;
566 val(1,0,op_dy) = -a1;
567 val(2,0,op_dy) = -a2;
568 }
569 }
570
Pi_h_alpha(const baseFElement & K,KN_<double> & v) const571 void TypeOfFE_RTortho::Pi_h_alpha(const baseFElement & K,KN_<double> & v) const
572 {
573 const Triangle & T(K.T);
574
575 for (int i=0,k=0;i<3;i++)
576 {
577 R2 E(T.Edge(i));
578 R signe = &T[ (i+1)%3] < &T[ (i+2)%3] ? 1.0 : -1.0 ;
579 v[k++]= signe*E.x;
580 v[k++]= signe*E.y;
581 }
582 }
583 // -------------------
584 // ttdc finite element fully discontinue.
585 // -------------------
586 class TypeOfFE_P1ttdc : public TypeOfFE { public:
587 static int Data[];
588 static double Pi_h_coef[];
589 static const R2 G;
590 static const R cshrink;
591 static const R cshrink1;
592 // (1 -1/3)*
593
Shrink(const R2 & P)594 static R2 Shrink(const R2& P){ return (P-G)*cshrink+G;}
Shrink1(const R2 & P)595 static R2 Shrink1(const R2& P){ return (P-G)*cshrink1+G;}
596
TypeOfFE_P1ttdc()597 TypeOfFE_P1ttdc(): TypeOfFE(0,0,3,1,Data,1,1,3,3,Pi_h_coef)
598 { const R2 Pt[] = { Shrink(R2(0,0)), Shrink(R2(1,0)), Shrink(R2(0,1)) };
599 for (int i=0;i<NbDoF;i++) {
600 pij_alpha[i]= IPJ(i,i,0);
601 P_Pi_h[i]=Pt[i];
602 // cout << Pt[i] << " " ;
603 }
604 // cout <<" cshrink: " << cshrink << " cshrink1 : "<< cshrink1 <<endl;
605 }
606
607 void FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
608
609
610 virtual R operator()(const FElement & K,const R2 & PHat,const KN_<R> & u,int componante,int op) const ;
611
612 } ;
613 const R2 TypeOfFE_P1ttdc::G(1./3.,1./3.);
614 const R TypeOfFE_P1ttdc::cshrink=1-1e-2;
615 const R TypeOfFE_P1ttdc::cshrink1=1./TypeOfFE_P1ttdc::cshrink;
616
617 class TypeOfFE_P2ttdc : public TypeOfFE { public:
618 static int Data[];
619 static double Pi_h_coef[];
620 static const R2 G;
621 static const R cshrink;
622 static const R cshrink1;
623
Shrink(const R2 & P)624 static R2 Shrink(const R2& P){ return (P-G)*cshrink+G;}
Shrink1(const R2 & P)625 static R2 Shrink1(const R2& P){ return (P-G)*cshrink1+G;}
626
TypeOfFE_P2ttdc()627 TypeOfFE_P2ttdc(): TypeOfFE(0,0,6,1,Data,3,1,6,6,Pi_h_coef)
628 { const R2 Pt[] = { Shrink(R2(0,0)), Shrink(R2(1,0)), Shrink(R2(0,1)),
629 Shrink(R2(0.5,0.5)),Shrink(R2(0,0.5)),Shrink(R2(0.5,0)) };
630 for (int i=0;i<NbDoF;i++) {
631 pij_alpha[i]= IPJ(i,i,0);
632 P_Pi_h[i]=Pt[i]; }
633 }
634
635
636 void FB(const bool * whatd,const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
637
638
639 } ;
640 // on what nu df on node node of df
641 int TypeOfFE_P1ttdc::Data[]={6,6,6, 0,1,2, 0,0,0, 0,0,0, 0,1,2, 0, 0,3};
642 int TypeOfFE_P2ttdc::Data[]={6,6,6,6,6,6, 0,1,2,3,4,5, 0,0,0,0,0,0, 0,0,0,0,0,0, 0,1,2,3,4,5, 0, 0,6};
643 double TypeOfFE_P1ttdc::Pi_h_coef[]={1.,1.,1.};
644 double TypeOfFE_P2ttdc::Pi_h_coef[]={1.,1.,1.,1.,1.,1.};
645
646 const R2 TypeOfFE_P2ttdc::G(1./3.,1./3.);
647 const R TypeOfFE_P2ttdc::cshrink=1-1e-2;
648 const R TypeOfFE_P2ttdc::cshrink1=1./TypeOfFE_P2ttdc::cshrink;
649
650
operator ()(const FElement & K,const R2 & P1Hat,const KN_<R> & u,int componante,int op) const651 R TypeOfFE_P1ttdc::operator()(const FElement & K,const R2 & P1Hat,const KN_<R> & u,int componante,int op) const
652 {
653
654 R2 PHat=Shrink1(P1Hat);
655 R u0(u(K(0))), u1(u(K(1))), u2(u(K(2)));
656 R r=0;
657 if (op==0)
658 {
659 R l0=1-PHat.x-PHat.y,l1=PHat.x,l2=PHat.y;
660 r = u0*l0+u1*l1+l2*u2;
661 }
662 else
663 {
664 const Triangle & T=K.T;
665 R2 D0 = T.H(0)*cshrink1 , D1 = T.H(1)*cshrink1 , D2 = T.H(2)*cshrink1 ;
666 if (op==1)
667 r = D0.x*u0 + D1.x*u1 + D2.x*u2 ;
668 else
669 r = D0.y*u0 + D1.y*u1 + D2.y*u2 ;
670 }
671 // cout << r << "\t";
672 return r;
673 }
674
FB(const bool * whatd,const Mesh &,const Triangle & K,const R2 & P1,RNMK_ & val) const675 void TypeOfFE_P1ttdc::FB(const bool *whatd,const Mesh & ,const Triangle & K,const R2 & P1,RNMK_ & val) const
676 {
677 R2 P=Shrink1(P1);
678
679 // const Triangle & K(FE.T);
680 R2 A(K[0]), B(K[1]),C(K[2]);
681 R l0=1-P.x-P.y,l1=P.x,l2=P.y;
682
683 if (val.N() <3)
684 throwassert(val.N() >=3);
685 throwassert(val.M()==1 );
686 // throwassert(val.K()==3 );
687
688 val=0;
689 RN_ f0(val('.',0,op_id));
690
691 if (whatd[op_id])
692 {
693 f0[0] = l0;
694 f0[1] = l1;
695 f0[2] = l2;}
696 if (whatd[op_dx] || whatd[op_dy])
697 {
698 R2 Dl0(K.H(0)*cshrink1), Dl1(K.H(1)*cshrink1), Dl2(K.H(2)*cshrink1);
699
700 if (whatd[op_dx])
701 {
702 RN_ f0x(val('.',0,op_dx));
703 f0x[0] = Dl0.x;
704 f0x[1] = Dl1.x;
705 f0x[2] = Dl2.x;
706 }
707
708 if (whatd[op_dy]) {
709 RN_ f0y(val('.',0,op_dy));
710 f0y[0] = Dl0.y;
711 f0y[1] = Dl1.y;
712 f0y[2] = Dl2.y;
713 }
714 }
715 }
716
717
718
FB(const bool * whatd,const Mesh &,const Triangle & K,const R2 & P1,RNMK_ & val) const719 void TypeOfFE_P2ttdc::FB(const bool *whatd,const Mesh & ,const Triangle & K,const R2 & P1,RNMK_ & val) const
720 {
721 R2 P=Shrink1(P1);
722
723 // const Triangle & K(FE.T);
724 R2 A(K[0]), B(K[1]),C(K[2]);
725 R l0=1-P.x-P.y,l1=P.x,l2=P.y;
726 R l4_0=(4*l0-1),l4_1=(4*l1-1),l4_2=(4*l2-1);
727
728 // throwassert(FE.N == 1);
729 throwassert( val.N()>=6);
730 throwassert(val.M()==1);
731 // throwassert(val.K()==3 );
732
733 val=0;
734 // --
735 if (whatd[op_id])
736 {
737 RN_ f0(val('.',0,op_id));
738 f0[0] = l0*(2*l0-1);
739 f0[1] = l1*(2*l1-1);
740 f0[2] = l2*(2*l2-1);
741 f0[3] = 4*l1*l2; // oppose au sommet 0
742 f0[4] = 4*l0*l2; // oppose au sommet 1
743 f0[5] = 4*l1*l0; // oppose au sommet 3
744 }
745 if( whatd[op_dx] || whatd[op_dy] || whatd[op_dxx] || whatd[op_dyy] || whatd[op_dxy])
746 {
747 R2 Dl0(K.H(0)*cshrink1), Dl1(K.H(1)*cshrink1), Dl2(K.H(2)*cshrink1);
748 if (whatd[op_dx])
749 {
750 RN_ f0x(val('.',0,op_dx));
751 f0x[0] = Dl0.x*l4_0;
752 f0x[1] = Dl1.x*l4_1;
753 f0x[2] = Dl2.x*l4_2;
754 f0x[3] = 4*(Dl1.x*l2 + Dl2.x*l1) ;
755 f0x[4] = 4*(Dl2.x*l0 + Dl0.x*l2) ;
756 f0x[5] = 4*(Dl0.x*l1 + Dl1.x*l0) ;
757 }
758
759 if (whatd[op_dy])
760 {
761 RN_ f0y(val('.',0,op_dy));
762 f0y[0] = Dl0.y*l4_0;
763 f0y[1] = Dl1.y*l4_1;
764 f0y[2] = Dl2.y*l4_2;
765 f0y[3] = 4*(Dl1.y*l2 + Dl2.y*l1) ;
766 f0y[4] = 4*(Dl2.y*l0 + Dl0.y*l2) ;
767 f0y[5] = 4*(Dl0.y*l1 + Dl1.y*l0) ;
768 }
769
770 if (whatd[op_dxx])
771 {
772 RN_ fxx(val('.',0,op_dxx));
773
774 fxx[0] = 4*Dl0.x*Dl0.x;
775 fxx[1] = 4*Dl1.x*Dl1.x;
776 fxx[2] = 4*Dl2.x*Dl2.x;
777 fxx[3] = 8*Dl1.x*Dl2.x;
778 fxx[4] = 8*Dl0.x*Dl2.x;
779 fxx[5] = 8*Dl0.x*Dl1.x;
780 }
781
782 if (whatd[op_dyy])
783 {
784 RN_ fyy(val('.',0,op_dyy));
785 fyy[0] = 4*Dl0.y*Dl0.y;
786 fyy[1] = 4*Dl1.y*Dl1.y;
787 fyy[2] = 4*Dl2.y*Dl2.y;
788 fyy[3] = 8*Dl1.y*Dl2.y;
789 fyy[4] = 8*Dl0.y*Dl2.y;
790 fyy[5] = 8*Dl0.y*Dl1.y;
791 }
792 if (whatd[op_dxy])
793 {
794 assert(val.K()>op_dxy);
795 RN_ fxy(val('.',0,op_dxy));
796
797 fxy[0] = 4*Dl0.x*Dl0.y;
798 fxy[1] = 4*Dl1.x*Dl1.y;
799 fxy[2] = 4*Dl2.x*Dl2.y;
800 fxy[3] = 4*(Dl1.x*Dl2.y + Dl1.y*Dl2.x);
801 fxy[4] = 4*(Dl0.x*Dl2.y + Dl0.y*Dl2.x);
802 fxy[5] = 4*(Dl0.x*Dl1.y + Dl0.y*Dl1.x);
803 }
804
805 }
806
807 }
808
809
810 //
811 // end ttdc
812 // ------------------
813
814 static TypeOfFE_RTortho The_TypeOfFE_RTortho;
815 static TypeOfFE_RT The_TypeOfFE_RT;
816 static TypeOfFE_P0 The_TypeOfFE_P0;
817 static TypeOfFE_P1ttdc The_TypeOfFE_P1ttdc;
818 static TypeOfFE_P2ttdc The_TypeOfFE_P2ttdc;
819 static TypeOfFE_RTmodif The_TypeOfFE_RTmodif;
820 static TypeOfFE_P1ncLagrange The_TypeOfFE_P1nc;
821 static TypeOfFE_ConsEdge The_TypeOfFE_ConsEdge; // add FH
822 TypeOfFE & RTLagrangeOrtho(The_TypeOfFE_RTortho);
823 TypeOfFE & RTLagrange(The_TypeOfFE_RT);
824 TypeOfFE & RTmodifLagrange(The_TypeOfFE_RTmodif);
825 TypeOfFE & P0Lagrange(The_TypeOfFE_P0);
826 TypeOfFE & P1ncLagrange(The_TypeOfFE_P1nc);
827 TypeOfFE & P1ttdc(The_TypeOfFE_P1ttdc);
828 TypeOfFE & P2ttdc(The_TypeOfFE_P2ttdc);
829 TypeOfFE & P0edge(The_TypeOfFE_ConsEdge);
830
831 }
832