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 "error.hpp"
29 #include "rgraph.hpp"
30 using namespace std;
31 #include "RNM.hpp"
32 #include "fem.hpp"
33 #include "FESpace.hpp"
34
35 namespace Fem2D {
36
37 // ------ P2h Hierarchical (just remove P1 node of the P2 finite element) --------
38 class TypeOfFE_P2hLagrange : public TypeOfFE { public:
39 static int Data[];
40 static double Pi_h_coef[];
41
TypeOfFE_P2hLagrange()42 TypeOfFE_P2hLagrange(): TypeOfFE(0,1,0,1,Data,4,1,3,3,Pi_h_coef)
43 {
44
45 const R2 Pt[] = { R2(0.5,0.5), R2(0.0,0.5), R2(0.5,0.0) };
46 for (int i=0;i<NbDoF;i++) {
47 pij_alpha[i]= IPJ(i,i,0);
48 P_Pi_h[i]=Pt[i]; }
49 }
50 void FB(const bool * whatd, const Mesh & Th,const Triangle & K,const R2 &P, RNMK_ & val) const;
51 } ;
52 // on what nu df on node node of df
53 int TypeOfFE_P2hLagrange::Data[]={3,4,5, 0,0,0, 0,1,2, 0,0,0, 0,1,2, 0, 0,3 };
54 double TypeOfFE_P2hLagrange::Pi_h_coef[]={1.,1.,1.};
55
FB(const bool * whatd,const Mesh &,const Triangle & K,const R2 & P,RNMK_ & val) const56 void TypeOfFE_P2hLagrange::FB(const bool * whatd,const Mesh & ,const Triangle & K,const R2 & P,RNMK_ & val) const
57 {
58 R2 A(K[0]), B(K[1]),C(K[2]);
59 R l0=1-P.x-P.y,l1=P.x,l2=P.y;
60 // R l4_0=(4*l0-1),l4_1=(4*l1-1),l4_2=(4*l2-1);
61
62 throwassert( val.N()>=3);
63 throwassert(val.M()==1);
64
65 val=0;
66 // --
67 if (whatd[op_id])
68 {
69 RN_ f0(val('.',0,op_id));
70 f0[0] = 4*l1*l2; // oppose au sommet 0
71 f0[1] = 4*l0*l2; // oppose au sommet 1
72 f0[2] = 4*l1*l0; // oppose au sommet 3
73 }
74 if( whatd[op_dx] || whatd[op_dy] || whatd[op_dxx] || whatd[op_dyy] || whatd[op_dxy])
75 {
76 R2 Dl0(K.H(0)), Dl1(K.H(1)), Dl2(K.H(2));
77 if (whatd[op_dx])
78 {
79 RN_ f0x(val('.',0,op_dx));
80 f0x[0] = 4*(Dl1.x*l2 + Dl2.x*l1) ;
81 f0x[1] = 4*(Dl2.x*l0 + Dl0.x*l2) ;
82 f0x[2] = 4*(Dl0.x*l1 + Dl1.x*l0) ;
83 }
84
85 if (whatd[op_dy])
86 {
87 RN_ f0y(val('.',0,op_dy));
88 f0y[0] = 4*(Dl1.y*l2 + Dl2.y*l1) ;
89 f0y[1] = 4*(Dl2.y*l0 + Dl0.y*l2) ;
90 f0y[2] = 4*(Dl0.y*l1 + Dl1.y*l0) ;
91 }
92
93 if (whatd[op_dxx])
94 {
95 RN_ fxx(val('.',0,op_dxx));
96
97 fxx[0] = 8*Dl1.x*Dl2.x;
98 fxx[1] = 8*Dl0.x*Dl2.x;
99 fxx[2] = 8*Dl0.x*Dl1.x;
100 }
101
102 if (whatd[op_dyy])
103 {
104 RN_ fyy(val('.',0,op_dyy));
105 fyy[0] = 8*Dl1.y*Dl2.y;
106 fyy[1] = 8*Dl0.y*Dl2.y;
107 fyy[2] = 8*Dl0.y*Dl1.y;
108 }
109 if (whatd[op_dxy])
110 {
111 assert(val.K()>op_dxy);
112 RN_ fxy(val('.',0,op_dxy));
113
114 fxy[0] = 4*(Dl1.x*Dl2.y + Dl1.y*Dl2.x);
115 fxy[1] = 4*(Dl0.x*Dl2.y + Dl0.y*Dl2.x);
116 fxy[2] = 4*(Dl0.x*Dl1.y + Dl0.y*Dl1.x);
117 }
118
119 }
120
121 }
122 // link with FreeFem++ do not work with static library .a
123 // FH so add a extern name to call in init_static_FE (see end of FESpace.cpp)
init_FE_P2h()124 void init_FE_P2h() { };
125
126 extern ListOfTFE typefem_P2h;
127
128 static TypeOfFE_P2hLagrange P2LagrangeP2h;
129 // given the name of the finite element in FreeFem++
130 ListOfTFE typefem_P2h("P2h", &P2LagrangeP2h);
131
132
133 // --- fin --
134 } // FEM2d namespace
135