v0.9.0
fem_tools.h
Go to the documentation of this file.
1 /** \file fem_tools.h
2  * \brief Loose implementation of some useful functions
3  *
4  * FIXME: Implementation here is very unstructured, need cleaning and pruning
5  *
6  */
7 
8 /*
9  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
12  * License for more details.
13  *
14  * You should have received a copy of the GNU Lesser General Public
15  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>
16 */
17 
18 #ifndef __FEM_H__
19 #define __FEM_H__
20 
21 #include<petscsys.h>
22 
23 #ifdef __cplusplus
24 extern "C" {
25 #endif
26  #include <cblas.h>
27  #include <lapack_wrap.h>
28 #ifdef __cplusplus
29 }
30 #endif
31 
32 #define LAMBDA(E,NU) (E*NU/((1.+NU)*(1.-2.*NU)))
33 #define MU(E,NU) (0.5*E/(1.+NU))
34 #define DELTA(NU_P,NU_PZ,E_P,E_Z) (((1+NU_P)*(1-NU_P-2*NU_PZ*(NU_PZ*E_Z/E_P)))/(E_P*E_P*E_Z))
35 
36 #define N_MBTET0(x, y, z) ( 1.-x-y-z ) ///< tetrahedral shape function
37 #define N_MBTET1(x, y, z) ( x ) ///< tetrahedral shape function
38 #define N_MBTET2(x, y, z) ( y ) ///< tetrahedral shape function
39 #define N_MBTET3(x, y, z) ( z ) ///< tetrahedral shape function
40 #define diffN_MBTET0x ( -1. ) ///< derivative of tetrahedral shape function
41 #define diffN_MBTET0y ( -1. ) ///< derivative of tetrahedral shape function
42 #define diffN_MBTET0z ( -1. ) ///< derivative of tetrahedral shape function
43 #define diffN_MBTET1x ( 1. ) ///< derivative of tetrahedral shape function
44 #define diffN_MBTET1y ( 0. ) ///< derivative of tetrahedral shape function
45 #define diffN_MBTET1z ( 0. ) ///< derivative of tetrahedral shape function
46 #define diffN_MBTET2x ( 0. ) ///< derivative of tetrahedral shape function
47 #define diffN_MBTET2y ( 1. ) ///< derivative of tetrahedral shape function
48 #define diffN_MBTET2z ( 0. ) ///< derivative of tetrahedral shape function
49 #define diffN_MBTET3x ( 0. ) ///< derivative of tetrahedral shape function
50 #define diffN_MBTET3y ( 0. ) ///< derivative of tetrahedral shape function
51 #define diffN_MBTET3z ( 1. ) ///< derivative of tetrahedral shape function
52 
53 //MBTRI
54 #define N_MBTRI0(x, y) ( 1.-x-y ) ///< triangle shape function
55 #define N_MBTRI1(x, y) ( x ) ///< triangle shape function
56 #define N_MBTRI2(x, y) ( y ) ///< triangle shape function
57 #define diffN_MBTRI0x ( -1. ) ///< derivative of triangle shape function
58 #define diffN_MBTRI0y ( -1. ) ///< derivative of triangle shape function
59 #define diffN_MBTRI1x ( 1. ) ///< derivative of triangle shape function
60 #define diffN_MBTRI1y ( 0. ) ///< derivative of triangle shape function
61 #define diffN_MBTRI2x ( 0. ) ///< derivative of triangle shape function
62 #define diffN_MBTRI2y ( 1. ) ///< derivative of triangle shape function
63 
64 //MBQUAD
65 #define N_MBQUAD0(x, y) ( (1.-x)*(1.-y) ) ///< quad shape function
66 #define N_MBQUAD1(x, y) ( x*(1-y) ) ///< quad shape function
67 #define N_MBQUAD2(x, y) ( x*y ) ///< quad shape function
68 #define N_MBQUAD3(x, y) ( (1.-x)*(y) ) ///< quad shape function
69 #define diffN_MBQUAD0x(y) ( -(1.-y) )
70 #define diffN_MBQUAD0y(x) ( -(1.-x) )
71 #define diffN_MBQUAD1x(y) ( (1.-y) )
72 #define diffN_MBQUAD1y(x) ( -x )
73 #define diffN_MBQUAD2x(y) ( y )
74 #define diffN_MBQUAD2y(x) ( x )
75 #define diffN_MBQUAD3x(y) ( -y )
76 #define diffN_MBQUAD3y(x) ( (1.-x) )
77 
78 //MBEDGE
79 #define N_MBEDGE0(x) ( 1.-(x) ) ///< edge shape function
80 #define N_MBEDGE1(x) ( x ) ///< edge shape function
81 #define diffN_MBEDGE0 ( -1. ) ///< derivative of edge shape function
82 #define diffN_MBEDGE1 ( 1. ) ///< derivative of edge shape function
83 
84 //MBTRIQ
85 #define N_MBTRIQ0(x, y) ( (1.-x-y)*(2*(1.-x-y)-1.) )
86 #define N_MBTRIQ1(x, y) ( x*(2.*x-1.) )
87 #define N_MBTRIQ2(x, y) ( y*(2.*y-1.) )
88 #define N_MBTRIQ3(x, y) ( 4.*(1.-x-y)*x )
89 #define N_MBTRIQ4(x, y) ( 4.*x*y )
90 #define N_MBTRIQ5(x, y) ( 4.*(1.-x-y)*y )
91 #define diffN_MBTRIQ0x(x, y) ( x+y-3.*(1.-x-y) )
92 #define diffN_MBTRIQ0y(x, y) ( x+y-3.*(1.-x-y) )
93 #define diffN_MBTRIQ1x(x, y) ( -1.+4.*x )
94 #define diffN_MBTRIQ1y(x, y) ( 0. )
95 #define diffN_MBTRIQ2x(x, y) ( 0. )
96 #define diffN_MBTRIQ2y(x, y) ( -1.+4.*y )
97 #define diffN_MBTRIQ3x(x, y) ( 4.*((1.-x-y)-x) )
98 #define diffN_MBTRIQ3y(x, y) ( -4.*x )
99 #define diffN_MBTRIQ4x(x, y) ( 4.*y )
100 #define diffN_MBTRIQ4y(x, y) ( 4.*x )
101 #define diffN_MBTRIQ5x(x, y) ( -4.*y )
102 #define diffN_MBTRIQ5y(x, y) ( 4.*((1.-x-y)-y) )
103 
104 #ifdef __cplusplus
105 extern "C" {
106 #endif
107 
108 /// print matric M
109 void print_mat(double *M,int m,int n);
110 /// print upper part of the symmetric matrix
111 void print_mat_sym_upper(double *M,int m,int n);
112 /// priint complex matrix
113 void print_mat_complex(__CLPK_doublecomplex *M,int m,int n);
114 
115 /// \brief calculate shape functions on triangle
116 /// \param N shape function array
117 /// \param X array of Gauss X coordinates
118 /// \param Y array of Gauss Y coordinates
119 /// \param G_DIM number of Gauss points
120 PetscErrorCode ShapeMBTRI(double *N,const double *X,const double *Y,const int G_DIM);
121 /// calculate derivatives of shape functions
122 PetscErrorCode ShapeDiffMBTRI(double *diffN);
123 
124 /// calculate face normal
125 /// \param diffN derivatives of shape functions
126 /// \param coords is position of the nodes
127 /// \param normal vector
128 PetscErrorCode ShapeFaceNormalMBTRI(double *diffN,const double *coords,double *normal);
129 PetscErrorCode ShapeFaceBaseMBTRI(
130  double *diffN,const double *coords,
131  double *normal,double *s1,double *s2);
132 
133 /// calculate derivative of normal in respect to nodal positions
134 PetscErrorCode ShapeFaceDiffNormalMBTRI(double *diffN,const double *coords,double *diff_normal);
135 /// calculate jacobioan
136 void ShapeJacMBTRI(double *diffN,const double *coords,double *Jac);
137 /// calculate derivatives of shape functions in space
138 void ShapeDiffMBTRIinvJ(double *diffN,double *invJac,double *diffNinvJac);
139 /// calculate shape functions
140 PetscErrorCode ShapeMBTET(double *N,const double *G_X,const double *G_Y,const double *G_Z,int DIM);
141 /// calculate derivatives of shape functions
142 PetscErrorCode ShapeDiffMBTET(double *diffN);
143 /// determined of jacobian
144 double ShapeDetJacVolume(double *jac);
145 /// calculate jacobian
146 PetscErrorCode ShapeJacMBTET(double *diffN,const double *coords,double *jac);
147 // calculate inverse of jacobian
148 PetscErrorCode ShapeInvJacVolume(double *jac);
149 /// calculate TET volume
150 double ShapeVolumeMBTET(double *diffN,const double *coords);
151 /// calculate shape functions derivatives in space
152 PetscErrorCode ShapeDiffMBTETinvJ(double *diffN,double *invJac,double *diffNinvJac);
153 
154 /// calculate spin matrix from vector
155 // \param spinOmega is a spin matrix
156 // \param vecOmega is a spin vector
157 PetscErrorCode Spin(double *spinOmega,double *vecOmega);
158 
159 /// Compose complex matrix (3x3) from two real matrices
160 PetscErrorCode make_complex_matrix(double *reA,double *imA,__CLPK_doublecomplex *xA);
161 /// Complex normal
162 PetscErrorCode Normal_hierarchical(
163  int order_approx,int *order_edge_approx,
164  int order,int *order_edge,
165  double *diffN,double *diffN_face,double *diffN_edge[],
166  double *dofs,double *dofs_edge[],double *dofs_face,
167  double *idofs,double *idofs_edge[],double *idofs_face,
168  __CLPK_doublecomplex *xnormal,
171  int gg);
172 PetscErrorCode Base_scale(
174 
175 /**
176  * \brief calculate local coordinates for given global coordinates
177  *
178  * new version for multiple points need to be implemented
179  */
180 PetscErrorCode ShapeMBTET_inverse(
181  double *N,double *diffN,const double *elem_coords,const double *glob_coords,double *loc_coords
182 );
183 
184 /**
185  * \brief calculate local coordinates of triangle element for given global coordinates in 2D (Assume e.g. z=0)
186  \f[
187  \left[\begin{array}{cc}
188 \frac{\partial N_{1}}{\partial\xi}x_{N_{1}}+\frac{\partial N_{2}}{\partial\xi}x_{N_{2}}+\frac{\partial N_{3}}{\partial\xi}x_{N_{3}} & \frac{\partial N_{1}}{\partial\eta}x_{N_{1}}+\frac{\partial N_{2}}{\partial\eta}x_{N_{2}}+\frac{\partial N_{3}}{\partial\eta}x_{N_{3}}\\
189 \frac{\partial N_{1}}{\partial\xi}y_{N_{1}}+\frac{\partial N_{2}}{\partial\xi}y_{N_{2}}+\frac{\partial N_{3}}{\partial\xi}y_{N_{3}} & \frac{\partial N_{1}}{\partial\eta}y_{N_{1}}+\frac{\partial N_{2}}{\partial\eta}y_{N_{2}}+\frac{\partial N_{3}}{\partial\eta}y_{N_{3}}
190 \end{array}\right]\left\{ \begin{array}{c}
191 \xi\\
192 \eta
193 \end{array}\right\} =\left\{ \begin{array}{c}
194 x_{gp}-\left(N_{1}x_{N_{1}}+N_{2}x_{N_{2}}+N_{3}x_{N_{3}}\right)\\
195 y_{gp}-\left(N_{1}y_{N_{1}}+N_{2}y_{N_{2}}+N_{3}y_{N_{3}}\right)
196 \end{array}\right\}
197  \f]
198 
199  /// \param N shape function array
200  /// \param diffN array of shape function derivative w.r.t to local coordinates
201  /// \param elem_coords global coordinates of element
202  /// \param glob_coords global coordinates of required point
203  /// \param loc_coords local coordinates of required point
204  */
205 PetscErrorCode ShapeMBTRI_inverse(
206  double *N,double *diffN,const double *elem_coords,const double *glob_coords,double *loc_coords
207 );
208 
209 /// calculate gradient of deformation
210 PetscErrorCode GradientOfDeformation(double *diffN,double *dofs,double *F);
211 
212 //2 Node edge
213 PetscErrorCode ShapeMBEDGE(double *N,const double *G_X,int DIM);
214 PetscErrorCode ShapeDiffMBEDGE(double *diffN);
215 
216 //10 Node Tet
217 PetscErrorCode ShapeMBTRIQ(double *N,const double *X,const double *Y,const int G_DIM);
218 PetscErrorCode ShapeDiffMBTRIQ(double *diffN,const double *X,const double *Y,const int G_DIM);
219 PetscErrorCode ShapeMBTETQ(double *N,const double x,const double y,const double z);
220 PetscErrorCode ShapeDiffMBTETQ(double *diffN,const double x,const double y,const double z);
221 PetscErrorCode ShapeMBTETQ_GAUSS(double *N,const double *X,const double *Y,const double *Z,const int G_DIM);
222 PetscErrorCode ShapeDiffMBTETQ_GAUSS(double *diffN,const double *X,const double *Y,const double *Z,const int G_DIM);
223 PetscErrorCode ShapeJacMBTETQ(const double *diffN,const double *coords,double *Jac);
224 PetscErrorCode ShapeMBTETQ_detJac_at_Gauss_Points(double *detJac_at_Gauss_Points,const double *diffN,const double *coords,int G_DIM);
225 double ShapeVolumeMBTETQ(const double *diffN,const double *coords,int G_DIM,double *G_TET_W);
226 PetscErrorCode ShapeMBTETQ_inverse(
227  double *N,double *diffN,const double *elem_coords,const double *glob_coords,double *loc_coords,const double eps
228 );
229 
230 //complex part
231 void ShapeDiffMBTETinvJ_complex(double *diffN,__CLPK_doublecomplex *invJac,__CLPK_doublecomplex *diffNinvJac,const enum CBLAS_TRANSPOSE Trans);
232 PetscErrorCode ShapeFaceNormalMBTRI_complex(double *diffN,__CLPK_doublecomplex *xcoords,__CLPK_doublecomplex *xnormal);
233 PetscErrorCode MakeComplexTensor(double *reA,double *imA,__CLPK_doublecomplex *xA);
234 PetscErrorCode InvertComplexGradient(__CLPK_doublecomplex *xF);
236 
237 //integration
238 ///Compute weights and integration points for edge using Grundmann_Moeller rule
239 PetscErrorCode Grundmann_Moeller_integration_points_1D_EDGE(int rule,double *G_TRI_X,double *G_TRI_W);
240 ///Compute weights and integration points for 2D Triangle using Grundmann_Moeller rule
241 PetscErrorCode Grundmann_Moeller_integration_points_2D_TRI(int rule,double *G_TRI_X,double *G_TRI_Y,double *G_TRI_W);
242 ///Compute weights and integration points for 3D Tet using Grundmann_Moeller rule
243 PetscErrorCode Grundmann_Moeller_integration_points_3D_TET(int rule,double *G_TET_X,double *G_TET_Y,double *G_TET_Z, double *G_TET_W);
244 
245 //http://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_tet/quadrature_rules_tet.html
246 //TRI
247 static const double G_TRI_X1[] = {
248 3.3333333333333331e-01
249 };
250 static const double G_TRI_Y1[] = {
251 3.3333333333333331e-01
252 };
253 static const double G_TRI_W1[] = {
254 1.
255 };
256 static const double G_TRI_X3[] = {
257 0.5, 0., 0.5
258 };
259 static const double G_TRI_Y3[] = {
260 0., 0.5, 0.5
261 };
262 static const double G_TRI_W3[] = {
263 3.3333333333333331e-01, 3.3333333333333331e-01, 3.3333333333333331e-01
264 };
265 static const double G_TRI_X4[] = {
266 7.503111022260811058e-02, 1.785587282636164064e-01,
267 2.800199154990741235e-01, 6.663902460147014262e-01
268 };
269 static const double G_TRI_Y4[] = {
270 2.800199154990741235e-01, 6.663902460147014262e-01,
271 7.503111022260811058e-02, 1.785587282636164064e-01
272 };
273 static const double G_TRI_W4[] = {
274 1.8195861825602258066e-01, 3.1804138174397683647e-01,
275 1.8195861825602258066e-01, 3.1804138174397683647e-01
276 };
277 static const double G_TRI_X7[] = {
278 0.333333333333333, 0.736712498968435, 0.736712498968435,
279 0.237932366472434, 0.237932366472434, 0.025355134551932,
280 0.025355134551932
281 };
282 static const double G_TRI_Y7[] = {
283 0.333333333333333, 0.237932366472434, 0.025355134551932,
284 0.736712498968435, 0.025355134551932, 0.736712498968435,
285 0.237932366472434
286 };
287 static const double G_TRI_W7[] = {
288  0.375000000000000, 0.104166666666667, 0.104166666666667, 0.104166666666667, 0.104166666666667, 0.104166666666667, 0.104166666666667
289 };
290 static const double G_TRI_X13[] = {
291 0.333333333333333, 0.479308067841923, 0.260345966079038, 0.260345966079038, 0.869739794195568,
292 0.065130102902216, 0.065130102902216, 0.638444188569809, 0.638444188569809, 0.312865496004875,
293 0.312865496004875, 0.048690315425316, 0.048690315425316
294 };
295 static const double G_TRI_Y13[] = {
296 0.333333333333333, 0.260345966079038, 0.479308067841923, 0.260345966079038, 0.065130102902216,
297 0.869739794195568, 0.065130102902216, 0.312865496004875, 0.048690315425316, 0.638444188569809,
298 0.048690315425316, 0.638444188569809, 0.312865496004875
299 };
300 static const double G_TRI_W13[] = {
301  -0.149570044467670, 0.175615257433204, 0.175615257433204, 0.175615257433204, 0.053347235608839,
302  0.053347235608839, 0.053347235608839, 0.077113760890257, 0.077113760890257, 0.077113760890257,
303  0.077113760890257, 0.077113760890257, 0.077113760890257 };
304 
305 static const double G_TRI_X19[] = {
306 0.333333333333333, 0.797426985353087, 0.101286507323456, 0.101286507323456, 0.059715871789770,
307 0.470142064105115, 0.470142064105115, 0.535795346449899, 0.232102326775050, 0.232102326775050,
308 0.941038278231121, 0.029480860884440, 0.029480860884440, 0.738416812340510, 0.738416812340510,
309 0.232102326775050, 0.232102326775050, 0.029480860884440, 0.029480860884440};
310 static const double G_TRI_Y19[] = {
311 0.333333333333333, 0.101286507323456, 0.797426985353087, 0.101286507323456, 0.470142064105115,
312 0.059715871789770, 0.470142064105115, 0.232102326775050, 0.535795346449899, 0.232102326775050,
313 0.029480860884440, 0.941038278231121, 0.029480860884440, 0.232102326775050, 0.029480860884440,
314 0.738416812340510, 0.029480860884440, 0.738416812340510, 0.232102326775050};
315 static const double G_TRI_W19[] = {
316 9.71357962827961025E-002, 3.13347002271398278E-002, 3.13347002271398278E-002, 3.13347002271398278E-002, 7.78275410047754301E-002,
317 7.78275410047754301E-002, 7.78275410047754301E-002, 7.96477389272090969E-002, 7.96477389272090969E-002, 7.96477389272090969E-002,
318 2.55776756586981006E-002, 2.55776756586981006E-002, 2.55776756586981006E-002, 4.32835393772893970E-002, 4.32835393772893970E-002,
319 4.32835393772893970E-002, 4.32835393772893970E-002, 4.32835393772893970E-002, 4.32835393772893970E-002};
320 
321 static const double G_TRI_X28[] = {
322 0.333333333333333, 0.948021718143423, 0.025989140928288, 0.025989140928288, 0.811424994704155,
323 0.094287502647923, 0.094287502647923, 0.010726449965571, 0.494636775017215, 0.494636775017215,
324 0.585313234770972, 0.207343382614514, 0.207343382614514, 0.122184388599019, 0.438907805700491,
325 0.438907805700491, 0.677937654882590, 0.677937654882590, 0.044841677589131, 0.044841677589131,
326 0.277220667528279, 0.277220667528279, 0.858870281282636, 0.858870281282636, 0.0000000000000000,
327 0.0000000000000000, 0.141129718717364, 0.141129718717364};
328 static const double G_TRI_Y28[] = {
329 0.333333333333333, 0.025989140928288, 0.948021718143423, 0.025989140928288, 0.094287502647923,
330 0.811424994704155, 0.094287502647923, 0.494636775017215, 0.010726449965571, 0.494636775017215,
331 0.207343382614514, 0.585313234770972, 0.207343382614514, 0.438907805700491, 0.122184388599019,
332 0.438907805700491, 0.044841677589131, 0.277220667528279, 0.677937654882590 , 0.277220667528279,
333 0.677937654882590, 0.044841677589131, 0.000000000000000, 0.141129718717364, 0.858870281282636,
334 0.141129718717364, 0.858870281282636, 0.000000000000000};
335 static const double G_TRI_W28[] = {
336 0.08797730116222190, 0.008744311553736190, 0.008744311553736190, 0.008744311553736190, 0.03808157199393533,
337 0.03808157199393533, 0.03808157199393533, 0.01885544805613125, 0.01885544805613125, 0.01885544805613125,
338 0.07215969754474100, 0.07215969754474100, 0.07215969754474100, 0.06932913870553720, 0.06932913870553720,
339 0.06932913870553720, 0.04105631542928860, 0.04105631542928860, 0.04105631542928860, 0.04105631542928860,
340 0.04105631542928860, 0.04105631542928860, 0.007362383783300573, 0.007362383783300573, 0.007362383783300573,
341 0.007362383783300573, 0.007362383783300573, 0.007362383783300573};
342 
343 static const double G_TRI_X37[] = {
344 0.333333333333333, 0.950275662924106, 0.024862168537947, 0.024862168537947, 0.171614914923835,
345 0.414192542538082, 0.414192542538082, 0.539412243677190, 0.230293878161405, 0.230293878161405,
346 0.772160036676533, 0.113919981661734, 0.113919981661734, 0.009085399949835, 0.495457300025082,
347 0.495457300025082, 0.062277290305887, 0.468861354847056, 0.468861354847056, 0.022076289653624,
348 0.022076289653624, 0.851306504174348, 0.851306504174348, 0.126617206172027, 0.126617206172027,
349 0.018620522802521, 0.018620522802521, 0.689441970728591, 0.689441970728591, 0.291937506468888,
350 0.291937506468888, 0.096506481292159, 0.096506481292159, 0.635867859433873, 0.635867859433873,
351 0.267625659273968, 0.267625659273968};
352 static const double G_TRI_Y37[] = {
353 0.333333333333333, 0.024862168537947, 0.950275662924106, 0.024862168537947, 0.414192542538082,
354 0.171614914923835, 0.414192542538082, 0.230293878161405, 0.539412243677190, 0.230293878161405,
355 0.113919981661734, 0.772160036676533, 0.113919981661734, 0.495457300025082, 0.009085399949835,
356 0.495457300025082, 0.468861354847056, 0.062277290305887, 0.468861354847056, 0.851306504174348,
357 0.126617206172027, 0.022076289653624, 0.126617206172027, 0.022076289653624, 0.851306504174348,
358 0.689441970728591, 0.291937506468888, 0.018620522802521, 0.291937506468888, 0.018620522802521,
359 0.689441970728591, 0.635867859433873, 0.267625659273968, 0.096506481292159, 0.267625659273968,
360 0.096506481292159, 0.635867859433873};
361 static const double G_TRI_W37[] = {
362 0.051739766065744, 0.008007799555565, 0.008007799555565, 0.008007799555565, 0.046868898981822,
363 0.046868898981822, 0.046868898981822, 0.046590940183976, 0.046590940183976, 0.046590940183976,
364 0.031016943313796, 0.031016943313796, 0.031016943313796, 0.010791612736631, 0.010791612736631,
365 0.010791612736631, 0.032195534242432, 0.032195534242432, 0.032195534242432, 0.015445834210702,
366 0.015445834210702, 0.015445834210702, 0.015445834210702, 0.015445834210702, 0.015445834210702,
367 0.017822989923179, 0.017822989923179, 0.017822989923179, 0.017822989923179, 0.017822989923179,
368 0.017822989923179, 0.037038683681385, 0.037038683681385, 0.037038683681385, 0.037038683681385,
369 0.037038683681385, 0.037038683681385};
370 static const double G_TRI_X286[] = {
371  0.04347826086956522, 0.1304347826086956, 0.2173913043478261, 0.3043478260869565, 0.391304347826087, 0.4782608695652174, 0.5652173913043478, 0.6521739130434783,
372  0.7391304347826086, 0.8260869565217391, 0.9130434782608695, 0.04347826086956522, 0.1304347826086956, 0.2173913043478261, 0.3043478260869565, 0.391304347826087,
373  0.4782608695652174, 0.5652173913043478, 0.6521739130434783, 0.7391304347826086, 0.8260869565217391, 0.04347826086956522, 0.1304347826086956, 0.2173913043478261,
374  0.3043478260869565, 0.391304347826087, 0.4782608695652174, 0.5652173913043478, 0.6521739130434783, 0.7391304347826086, 0.04347826086956522, 0.1304347826086956,
375  0.2173913043478261, 0.3043478260869565, 0.391304347826087, 0.4782608695652174, 0.5652173913043478, 0.6521739130434783, 0.04347826086956522, 0.1304347826086956,
376  0.2173913043478261, 0.3043478260869565, 0.391304347826087, 0.4782608695652174, 0.5652173913043478, 0.04347826086956522, 0.1304347826086956, 0.2173913043478261,
377  0.3043478260869565, 0.391304347826087, 0.4782608695652174, 0.04347826086956522, 0.1304347826086956, 0.2173913043478261, 0.3043478260869565, 0.391304347826087,
378  0.04347826086956522, 0.1304347826086956, 0.2173913043478261, 0.3043478260869565, 0.04347826086956522, 0.1304347826086956, 0.2173913043478261, 0.04347826086956522,
379  0.1304347826086956, 0.04347826086956522, 0.04761904761904762, 0.1428571428571428, 0.2380952380952381, 0.3333333333333333, 0.4285714285714285, 0.5238095238095238,
380  0.6190476190476191, 0.7142857142857143, 0.8095238095238095, 0.9047619047619048, 0.04761904761904762, 0.1428571428571428, 0.2380952380952381, 0.3333333333333333,
381  0.4285714285714285, 0.5238095238095238, 0.6190476190476191, 0.7142857142857143, 0.8095238095238095, 0.04761904761904762, 0.1428571428571428, 0.2380952380952381,
382  0.3333333333333333, 0.4285714285714285, 0.5238095238095238, 0.6190476190476191, 0.7142857142857143, 0.04761904761904762, 0.1428571428571428, 0.2380952380952381,
383  0.3333333333333333, 0.4285714285714285, 0.5238095238095238, 0.6190476190476191, 0.04761904761904762, 0.1428571428571428, 0.2380952380952381, 0.3333333333333333,
384  0.4285714285714285, 0.5238095238095238, 0.04761904761904762, 0.1428571428571428, 0.2380952380952381, 0.3333333333333333, 0.4285714285714285, 0.04761904761904762,
385  0.1428571428571428, 0.2380952380952381, 0.3333333333333333, 0.04761904761904762, 0.1428571428571428, 0.2380952380952381, 0.04761904761904762, 0.1428571428571428,
386  0.04761904761904762, 0.05263157894736842, 0.1578947368421053, 0.2631578947368421, 0.3684210526315789, 0.4736842105263158, 0.5789473684210527, 0.6842105263157895,
387  0.7894736842105263, 0.8947368421052632, 0.05263157894736842, 0.1578947368421053, 0.2631578947368421, 0.3684210526315789, 0.4736842105263158, 0.5789473684210527,
388  0.6842105263157895, 0.7894736842105263, 0.05263157894736842, 0.1578947368421053, 0.2631578947368421, 0.3684210526315789, 0.4736842105263158, 0.5789473684210527,
389  0.6842105263157895, 0.05263157894736842, 0.1578947368421053, 0.2631578947368421, 0.3684210526315789, 0.4736842105263158, 0.5789473684210527, 0.05263157894736842,
390  0.1578947368421053, 0.2631578947368421, 0.3684210526315789, 0.4736842105263158, 0.05263157894736842, 0.1578947368421053, 0.2631578947368421, 0.3684210526315789,
391  0.05263157894736842, 0.1578947368421053, 0.2631578947368421, 0.05263157894736842, 0.1578947368421053, 0.05263157894736842, 0.05882352941176471, 0.1764705882352941,
392  0.2941176470588235, 0.4117647058823529, 0.5294117647058824, 0.6470588235294118, 0.7647058823529411, 0.8823529411764706, 0.05882352941176471, 0.1764705882352941,
393  0.2941176470588235, 0.4117647058823529, 0.5294117647058824, 0.6470588235294118, 0.7647058823529411, 0.05882352941176471, 0.1764705882352941, 0.2941176470588235,
394  0.4117647058823529, 0.5294117647058824, 0.6470588235294118, 0.05882352941176471, 0.1764705882352941, 0.2941176470588235, 0.4117647058823529, 0.5294117647058824,
395  0.05882352941176471, 0.1764705882352941, 0.2941176470588235, 0.4117647058823529, 0.05882352941176471, 0.1764705882352941, 0.2941176470588235, 0.05882352941176471,
396  0.1764705882352941, 0.05882352941176471, 0.06666666666666667, 0.2, 0.3333333333333333, 0.4666666666666667, 0.6, 0.7333333333333333,
397  0.8666666666666667, 0.06666666666666667, 0.2, 0.3333333333333333, 0.4666666666666667, 0.6, 0.7333333333333333, 0.06666666666666667,
398  0.2, 0.3333333333333333, 0.4666666666666667, 0.6, 0.06666666666666667, 0.2, 0.3333333333333333, 0.4666666666666667,
399  0.06666666666666667, 0.2, 0.3333333333333333, 0.06666666666666667, 0.2, 0.06666666666666667, 0.07692307692307693, 0.2307692307692308,
400  0.3846153846153846, 0.5384615384615384, 0.6923076923076923, 0.8461538461538461, 0.07692307692307693, 0.2307692307692308, 0.3846153846153846, 0.5384615384615384,
401  0.6923076923076923, 0.07692307692307693, 0.2307692307692308, 0.3846153846153846, 0.5384615384615384, 0.07692307692307693, 0.2307692307692308, 0.3846153846153846,
402  0.07692307692307693, 0.2307692307692308, 0.07692307692307693, 0.09090909090909091, 0.2727272727272727, 0.4545454545454545, 0.6363636363636364, 0.8181818181818182,
403  0.09090909090909091, 0.2727272727272727, 0.4545454545454545, 0.6363636363636364, 0.09090909090909091, 0.2727272727272727, 0.4545454545454545, 0.09090909090909091,
404  0.2727272727272727, 0.09090909090909091, 0.1111111111111111, 0.3333333333333333, 0.5555555555555556, 0.7777777777777778, 0.1111111111111111, 0.3333333333333333,
405  0.5555555555555556, 0.1111111111111111, 0.3333333333333333, 0.1111111111111111, 0.1428571428571428, 0.4285714285714285, 0.7142857142857143, 0.1428571428571428,
406  0.4285714285714285, 0.1428571428571428, 0.2, 0.6, 0.2, 0.3333333333333333
407 };
408 static const double G_TRI_Y286[] = {
409  0.04347826086956522, 0.04347826086956522, 0.04347826086956522, 0.04347826086956522, 0.04347826086956522, 0.04347826086956522, 0.04347826086956522, 0.04347826086956522,
410  0.04347826086956522, 0.04347826086956522, 0.04347826086956522, 0.1304347826086956, 0.1304347826086956, 0.1304347826086956, 0.1304347826086956, 0.1304347826086956,
411  0.1304347826086956, 0.1304347826086956, 0.1304347826086956, 0.1304347826086956, 0.1304347826086956, 0.2173913043478261, 0.2173913043478261, 0.2173913043478261,
412  0.2173913043478261, 0.2173913043478261, 0.2173913043478261, 0.2173913043478261, 0.2173913043478261, 0.2173913043478261, 0.3043478260869565, 0.3043478260869565,
413  0.3043478260869565, 0.3043478260869565, 0.3043478260869565, 0.3043478260869565, 0.3043478260869565, 0.3043478260869565, 0.391304347826087, 0.391304347826087,
414  0.391304347826087, 0.391304347826087, 0.391304347826087, 0.391304347826087, 0.391304347826087, 0.4782608695652174, 0.4782608695652174, 0.4782608695652174,
415  0.4782608695652174, 0.4782608695652174, 0.4782608695652174, 0.5652173913043478, 0.5652173913043478, 0.5652173913043478, 0.5652173913043478, 0.5652173913043478,
416  0.6521739130434783, 0.6521739130434783, 0.6521739130434783, 0.6521739130434783, 0.7391304347826086, 0.7391304347826086, 0.7391304347826086, 0.8260869565217391,
417  0.8260869565217391, 0.9130434782608695, 0.04761904761904762, 0.04761904761904762, 0.04761904761904762, 0.04761904761904762, 0.04761904761904762, 0.04761904761904762,
418  0.04761904761904762, 0.04761904761904762, 0.04761904761904762, 0.04761904761904762, 0.1428571428571428, 0.1428571428571428, 0.1428571428571428, 0.1428571428571428,
419  0.1428571428571428, 0.1428571428571428, 0.1428571428571428, 0.1428571428571428, 0.1428571428571428, 0.2380952380952381, 0.2380952380952381, 0.2380952380952381,
420  0.2380952380952381, 0.2380952380952381, 0.2380952380952381, 0.2380952380952381, 0.2380952380952381, 0.3333333333333333, 0.3333333333333333, 0.3333333333333333,
421  0.3333333333333333, 0.3333333333333333, 0.3333333333333333, 0.3333333333333333, 0.4285714285714285, 0.4285714285714285, 0.4285714285714285, 0.4285714285714285,
422  0.4285714285714285, 0.4285714285714285, 0.5238095238095238, 0.5238095238095238, 0.5238095238095238, 0.5238095238095238, 0.5238095238095238, 0.6190476190476191,
423  0.6190476190476191, 0.6190476190476191, 0.6190476190476191, 0.7142857142857143, 0.7142857142857143, 0.7142857142857143, 0.8095238095238095, 0.8095238095238095,
424  0.9047619047619048, 0.05263157894736842, 0.05263157894736842, 0.05263157894736842, 0.05263157894736842, 0.05263157894736842, 0.05263157894736842, 0.05263157894736842,
425  0.05263157894736842, 0.05263157894736842, 0.1578947368421053, 0.1578947368421053, 0.1578947368421053, 0.1578947368421053, 0.1578947368421053, 0.1578947368421053,
426  0.1578947368421053, 0.1578947368421053, 0.2631578947368421, 0.2631578947368421, 0.2631578947368421, 0.2631578947368421, 0.2631578947368421, 0.2631578947368421,
427  0.2631578947368421, 0.3684210526315789, 0.3684210526315789, 0.3684210526315789, 0.3684210526315789, 0.3684210526315789, 0.3684210526315789, 0.4736842105263158,
428  0.4736842105263158, 0.4736842105263158, 0.4736842105263158, 0.4736842105263158, 0.5789473684210527, 0.5789473684210527, 0.5789473684210527, 0.5789473684210527,
429  0.6842105263157895, 0.6842105263157895, 0.6842105263157895, 0.7894736842105263, 0.7894736842105263, 0.8947368421052632, 0.05882352941176471, 0.05882352941176471,
430  0.05882352941176471, 0.05882352941176471, 0.05882352941176471, 0.05882352941176471, 0.05882352941176471, 0.05882352941176471, 0.1764705882352941, 0.1764705882352941,
431  0.1764705882352941, 0.1764705882352941, 0.1764705882352941, 0.1764705882352941, 0.1764705882352941, 0.2941176470588235, 0.2941176470588235, 0.2941176470588235,
432  0.2941176470588235, 0.2941176470588235, 0.2941176470588235, 0.4117647058823529, 0.4117647058823529, 0.4117647058823529, 0.4117647058823529, 0.4117647058823529,
433  0.5294117647058824, 0.5294117647058824, 0.5294117647058824, 0.5294117647058824, 0.6470588235294118, 0.6470588235294118, 0.6470588235294118, 0.7647058823529411,
434  0.7647058823529411, 0.8823529411764706, 0.06666666666666667, 0.06666666666666667, 0.06666666666666667, 0.06666666666666667, 0.06666666666666667, 0.06666666666666667,
435  0.06666666666666667, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.3333333333333333,
436  0.3333333333333333, 0.3333333333333333, 0.3333333333333333, 0.3333333333333333, 0.4666666666666667, 0.4666666666666667, 0.4666666666666667, 0.4666666666666667,
437  0.6, 0.6, 0.6, 0.7333333333333333, 0.7333333333333333, 0.8666666666666667, 0.07692307692307693, 0.07692307692307693,
438  0.07692307692307693, 0.07692307692307693, 0.07692307692307693, 0.07692307692307693, 0.2307692307692308, 0.2307692307692308, 0.2307692307692308, 0.2307692307692308,
439  0.2307692307692308, 0.3846153846153846, 0.3846153846153846, 0.3846153846153846, 0.3846153846153846, 0.5384615384615384, 0.5384615384615384, 0.5384615384615384,
440  0.6923076923076923, 0.6923076923076923, 0.8461538461538461, 0.09090909090909091, 0.09090909090909091, 0.09090909090909091, 0.09090909090909091, 0.09090909090909091,
441  0.2727272727272727, 0.2727272727272727, 0.2727272727272727, 0.2727272727272727, 0.4545454545454545, 0.4545454545454545, 0.4545454545454545, 0.6363636363636364,
442  0.6363636363636364, 0.8181818181818182, 0.1111111111111111, 0.1111111111111111, 0.1111111111111111, 0.1111111111111111, 0.3333333333333333, 0.3333333333333333,
443  0.3333333333333333, 0.5555555555555556, 0.5555555555555556, 0.7777777777777778, 0.1428571428571428, 0.1428571428571428, 0.1428571428571428, 0.4285714285714285,
444  0.4285714285714285, 0.7142857142857143, 0.2, 0.2, 0.6, 0.3333333333333333
445 };
446  static const double G_TRI_W286[] = {
447  2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668,
448  2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668,
449  2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668,
450  2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668,
451  2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668,
452  2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668,
453  2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668,
454  2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668, 2.912193380035668,
455  2.912193380035668, 2.912193380035668, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852,
456  -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852,
457  -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852,
458  -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852,
459  -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852,
460  -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852,
461  -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852, -9.914451197589852,
462  -9.914451197589852, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992,
463  13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992,
464  13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992,
465  13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992,
466  13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992,
467  13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, 13.33158527957992, -9.027792408986382, -9.027792408986382,
468  -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382,
469  -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382,
470  -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382,
471  -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382, -9.027792408986382,
472  -9.027792408986382, -9.027792408986382, 3.258672079964582, 3.258672079964582, 3.258672079964582, 3.258672079964582, 3.258672079964582, 3.258672079964582,
473  3.258672079964582, 3.258672079964582, 3.258672079964582, 3.258672079964582, 3.258672079964582, 3.258672079964582, 3.258672079964582, 3.258672079964582,
474  3.258672079964582, 3.258672079964582, 3.258672079964582, 3.258672079964582, 3.258672079964582, 3.258672079964582, 3.258672079964582, 3.258672079964582,
475  3.258672079964582, 3.258672079964582, 3.258672079964582, 3.258672079964582, 3.258672079964582, 3.258672079964582, -0.6133639040302452, -0.6133639040302452,
476  -0.6133639040302452, -0.6133639040302452, -0.6133639040302452, -0.6133639040302452, -0.6133639040302452, -0.6133639040302452, -0.6133639040302452, -0.6133639040302452,
477  -0.6133639040302452, -0.6133639040302452, -0.6133639040302452, -0.6133639040302452, -0.6133639040302452, -0.6133639040302452, -0.6133639040302452, -0.6133639040302452,
478  -0.6133639040302452, -0.6133639040302452, -0.6133639040302452, 0.05511571669513555, 0.05511571669513555, 0.05511571669513555, 0.05511571669513555, 0.05511571669513555,
479  0.05511571669513555, 0.05511571669513555, 0.05511571669513555, 0.05511571669513555, 0.05511571669513555, 0.05511571669513555, 0.05511571669513555, 0.05511571669513555,
480  0.05511571669513555, 0.05511571669513555, -0.001979122382447095, -0.001979122382447095, -0.001979122382447095, -0.001979122382447095, -0.001979122382447095, -0.001979122382447095,
481  -0.001979122382447095, -0.001979122382447095, -0.001979122382447095, -0.001979122382447095, 2.02054621415273e-05, 2.02054621415273e-05, 2.02054621415273e-05, 2.02054621415273e-05,
482  2.02054621415273e-05, 2.02054621415273e-05, -2.874940020535803e-08, -2.874940020535803e-08, -2.874940020535803e-08, 8.829438425435718e-13
483 };
484 
485 //TET
486 static const double G_TET_X1[] = {
487 0.25
488 };
489 static const double G_TET_Y1[] = {
490 0.25
491 };
492 static const double G_TET_Z1[] = {
493 0.25
494 };
495 static const double G_TET_W1[] = {
496 1.
497 };
498 static const double G_TET_X4[] = {
499 0.1757281246520584, 0.2445310270213291,
500 0.5556470949048655, 0.0240937534217468
501 };
502 static const double G_TET_Y4[] = {
503 0.5656137776620919, 0.0501800797762026,
504 0.1487681308666864, 0.2354380116950194
505 };
506 static const double G_TET_Z4[] = {
507 0.2180665126782654, 0.5635595064952189,
508 0.0350112499848832, 0.1833627308416330
509 };
510 static const double G_TET_W4[] = {
511 0.25,0.25,0.25,0.25
512 };
513 static const double G_TET_X5[] = {
514 0.25000000000000000, 0.50000000000000000,
515 0.16666666666666667, 0.16666666666666667,
516 0.16666666666666667
517 };
518 static const double G_TET_Y5[] = {
519 0.25000000000000000, 0.16666666666666667,
520 0.50000000000000000, 0.16666666666666667,
521 0.16666666666666667
522 };
523 static const double G_TET_Z5[] = {
524 0.25000000000000000, 0.16666666666666667,
525 0.16666666666666667, 0.50000000000000000,
526 0.16666666666666667
527 };
528 static const double G_TET_W5[] = {
529 -0.80000000000000000, 0.45000000000000000,
530 0.45000000000000000, 0.45000000000000000,
531 0.45000000000000000
532 };
533 static const double G_TET_X10[] = {
534  0.5684305841968444, 0.1438564719343852, 0.1438564719343852, 0.1438564719343852,
535  0.0000000000000000, 0.5000000000000000, 0.5000000000000000, 0.5000000000000000,
536  0.0000000000000000, 0.0000000000000000
537 };
538 static const double G_TET_Y10[] = {
539  0.1438564719343852, 0.1438564719343852, 0.1438564719343852, 0.5684305841968444,
540  0.5000000000000000, 0.0000000000000000, 0.5000000000000000, 0.0000000000000000,
541  0.5000000000000000, 0.0000000000000000
542 };
543 static const double G_TET_Z10[] = {
544  0.1438564719343852, 0.1438564719343852, 0.5684305841968444, 0.1438564719343852,
545  0.5000000000000000, 0.5000000000000000, 0.0000000000000000, 0.0000000000000000,
546  0.0000000000000000, 0.5000000000000000
547 };
548 static const double G_TET_W10[] = {
549  0.2177650698804054, 0.2177650698804054, 0.2177650698804054, 0.2177650698804054,
550  0.0214899534130631, 0.0214899534130631, 0.0214899534130631, 0.0214899534130631,
551  0.0214899534130631, 0.0214899534130631
552 };
553 
554 static const double G_TET_X45[] = {
555 0.2500000000000000, 0.6175871903000830, 0.1274709365666390, 0.1274709365666390, 0.1274709365666390,
556 0.9037635088221031, 0.0320788303926323, 0.0320788303926323, 0.0320788303926323, 0.4502229043567190,
557 0.0497770956432810, 0.0497770956432810, 0.0497770956432810, 0.4502229043567190, 0.4502229043567190,
558 0.3162695526014501, 0.1837304473985499, 0.1837304473985499, 0.1837304473985499, 0.3162695526014501,
559 0.3162695526014501, 0.0229177878448171, 0.2319010893971509, 0.2319010893971509, 0.5132800333608811,
560 0.2319010893971509, 0.2319010893971509, 0.2319010893971509, 0.0229177878448171, 0.5132800333608811,
561 0.2319010893971509, 0.0229177878448171, 0.5132800333608811, 0.7303134278075384, 0.0379700484718286,
562 0.0379700484718286, 0.1937464752488044, 0.0379700484718286, 0.0379700484718286, 0.0379700484718286,
563 0.7303134278075384, 0.1937464752488044, 0.0379700484718286, 0.7303134278075384, 0.1937464752488044
564 };
565 static const double G_TET_Y45[] = {
566 0.2500000000000000, 0.1274709365666390, 0.1274709365666390, 0.1274709365666390, 0.6175871903000830,
567 0.0320788303926323, 0.0320788303926323, 0.0320788303926323, 0.9037635088221031, 0.0497770956432810,
568 0.4502229043567190, 0.0497770956432810, 0.4502229043567190, 0.0497770956432810, 0.4502229043567190,
569 0.1837304473985499, 0.3162695526014501, 0.1837304473985499, 0.3162695526014501, 0.1837304473985499,
570 0.3162695526014501, 0.2319010893971509, 0.0229177878448171, 0.2319010893971509, 0.2319010893971509,
571 0.5132800333608811, 0.2319010893971509, 0.0229177878448171, 0.5132800333608811, 0.2319010893971509,
572 0.5132800333608811, 0.2319010893971509, 0.0229177878448171, 0.0379700484718286, 0.7303134278075384,
573 0.0379700484718286, 0.0379700484718286, 0.1937464752488044, 0.0379700484718286, 0.7303134278075384,
574 0.1937464752488044, 0.0379700484718286, 0.1937464752488044, 0.0379700484718286, 0.7303134278075384
575 };
576 static const double G_TET_Z45[] = {
577 0.2500000000000000, 0.1274709365666390, 0.1274709365666390, 0.6175871903000830, 0.1274709365666390,
578 0.0320788303926323, 0.0320788303926323, 0.9037635088221031, 0.0320788303926323, 0.0497770956432810,
579 0.0497770956432810, 0.4502229043567190, 0.4502229043567190, 0.4502229043567190, 0.0497770956432810,
580 0.1837304473985499, 0.1837304473985499, 0.3162695526014501, 0.3162695526014501, 0.3162695526014501,
581 0.1837304473985499, 0.2319010893971509, 0.2319010893971509, 0.0229177878448171, 0.2319010893971509,
582 0.2319010893971509, 0.5132800333608811, 0.5132800333608811, 0.2319010893971509, 0.0229177878448171,
583 0.0229177878448171, 0.5132800333608811, 0.2319010893971509, 0.0379700484718286, 0.0379700484718286,
584 0.7303134278075384, 0.0379700484718286, 0.0379700484718286, 0.1937464752488044, 0.1937464752488044,
585 0.0379700484718286, 0.7303134278075384, 0.7303134278075384, 0.1937464752488044, 0.0379700484718286
586 };
587 static const double G_TET_W45[] = {
588  -0.2359620398477557, 0.0244878963560562, 0.0244878963560562, 0.0244878963560562, 0.0244878963560562,
589  0.0039485206398261, 0.0039485206398261, 0.0039485206398261, 0.0039485206398261, 0.0263055529507371,
590  0.0263055529507371, 0.0263055529507371, 0.0263055529507371, 0.0263055529507371, 0.0263055529507371,
591  0.0829803830550589, 0.0829803830550589, 0.0829803830550589, 0.0829803830550589, 0.0829803830550589,
592  0.0829803830550589, 0.0254426245481023, 0.0254426245481023, 0.0254426245481023, 0.0254426245481023,
593  0.0254426245481023, 0.0254426245481023, 0.0254426245481023, 0.0254426245481023, 0.0254426245481023,
594  0.0254426245481023, 0.0254426245481023, 0.0254426245481023, 0.0134324384376852, 0.0134324384376852,
595  0.0134324384376852, 0.0134324384376852, 0.0134324384376852, 0.0134324384376852, 0.0134324384376852,
596  0.0134324384376852, 0.0134324384376852, 0.0134324384376852, 0.0134324384376852, 0.0134324384376852
597 };
598 static const double NC_TET_X84[] = {
599  0.1000000000000000, 0.1000000000000000, 0.1000000000000000, 0.7000000000000000, 0.1000000000000000, 0.1000000000000000, 0.1000000000000000,
600  0.1000000000000000, 0.1000000000000000, 0.1000000000000000, 0.6000000000000000, 0.6000000000000000, 0.6000000000000000, 0.2000000000000000,
601  0.2000000000000000, 0.2000000000000000, 0.1000000000000000, 0.1000000000000000, 0.1000000000000000, 0.1000000000000000, 0.1000000000000000,
602  0.1000000000000000, 0.5000000000000000, 0.5000000000000000, 0.5000000000000000, 0.3000000000000000, 0.3000000000000000, 0.3000000000000000,
603  0.2000000000000000, 0.2000000000000000, 0.2000000000000000, 0.2000000000000000, 0.2000000000000000, 0.2000000000000000, 0.1000000000000000,
604  0.1000000000000000, 0.1000000000000000, 0.5000000000000000, 0.5000000000000000, 0.5000000000000000, 0.4000000000000000, 0.4000000000000000,
605  0.4000000000000000, 0.1000000000000000, 0.1000000000000000, 0.1000000000000000, 0.4000000000000000, 0.4000000000000000, 0.4000000000000000,
606  0.4000000000000000, 0.4000000000000000, 0.4000000000000000, 0.3000000000000000, 0.3000000000000000, 0.3000000000000000, 0.3000000000000000,
607  0.3000000000000000, 0.3000000000000000, 0.2000000000000000, 0.2000000000000000, 0.2000000000000000, 0.2000000000000000, 0.2000000000000000,
608  0.2000000000000000, 0.1000000000000000, 0.1000000000000000, 0.1000000000000000, 0.1000000000000000, 0.1000000000000000, 0.1000000000000000,
609  0.2000000000000000, 0.2000000000000000, 0.2000000000000000, 0.4000000000000000, 0.3000000000000000, 0.3000000000000000, 0.3000000000000000,
610  0.1000000000000000, 0.3000000000000000, 0.3000000000000000, 0.3000000000000000, 0.2000000000000000, 0.2000000000000000, 0.2000000000000000
611 };
612 static const double NC_TET_Y84[] = {
613  0.1000000000000000, 0.1000000000000000, 0.7000000000000000, 0.1000000000000000, 0.1000000000000000, 0.1000000000000000, 0.6000000000000000,
614  0.6000000000000000, 0.2000000000000000, 0.2000000000000000, 0.1000000000000000, 0.1000000000000000, 0.2000000000000000, 0.1000000000000000,
615  0.1000000000000000, 0.6000000000000000, 0.1000000000000000, 0.1000000000000000, 0.5000000000000000, 0.5000000000000000, 0.3000000000000000,
616  0.3000000000000000, 0.1000000000000000, 0.1000000000000000, 0.3000000000000000, 0.1000000000000000, 0.1000000000000000, 0.5000000000000000,
617  0.2000000000000000, 0.2000000000000000, 0.1000000000000000, 0.1000000000000000, 0.5000000000000000, 0.5000000000000000, 0.2000000000000000,
618  0.2000000000000000, 0.5000000000000000, 0.2000000000000000, 0.2000000000000000, 0.1000000000000000, 0.4000000000000000, 0.1000000000000000,
619  0.1000000000000000, 0.4000000000000000, 0.4000000000000000, 0.1000000000000000, 0.3000000000000000, 0.3000000000000000, 0.2000000000000000,
620  0.2000000000000000, 0.1000000000000000, 0.1000000000000000, 0.4000000000000000, 0.4000000000000000, 0.2000000000000000, 0.2000000000000000,
621  0.1000000000000000, 0.1000000000000000, 0.4000000000000000, 0.4000000000000000, 0.3000000000000000, 0.3000000000000000, 0.1000000000000000,
622  0.1000000000000000, 0.4000000000000000, 0.4000000000000000, 0.3000000000000000, 0.3000000000000000, 0.2000000000000000, 0.2000000000000000,
623  0.2000000000000000, 0.2000000000000000, 0.4000000000000000, 0.2000000000000000, 0.3000000000000000, 0.3000000000000000, 0.1000000000000000,
624  0.3000000000000000, 0.3000000000000000, 0.2000000000000000, 0.2000000000000000, 0.3000000000000000, 0.3000000000000000, 0.2000000000000000
625 };
626 static const double NC_TET_Z84[] = {
627  0.1000000000000000, 0.7000000000000000, 0.1000000000000000, 0.1000000000000000, 0.6000000000000000, 0.2000000000000000, 0.1000000000000000,
628  0.2000000000000000, 0.1000000000000000, 0.6000000000000000, 0.1000000000000000, 0.2000000000000000, 0.1000000000000000, 0.1000000000000000,
629  0.6000000000000000, 0.1000000000000000, 0.5000000000000000, 0.3000000000000000, 0.1000000000000000, 0.3000000000000000, 0.1000000000000000,
630  0.5000000000000000, 0.1000000000000000, 0.3000000000000000, 0.1000000000000000, 0.1000000000000000, 0.5000000000000000, 0.1000000000000000,
631  0.1000000000000000, 0.5000000000000000, 0.2000000000000000, 0.5000000000000000, 0.2000000000000000, 0.1000000000000000, 0.2000000000000000,
632  0.5000000000000000, 0.2000000000000000, 0.2000000000000000, 0.1000000000000000, 0.2000000000000000, 0.1000000000000000, 0.4000000000000000,
633  0.1000000000000000, 0.4000000000000000, 0.1000000000000000, 0.4000000000000000, 0.2000000000000000, 0.1000000000000000, 0.3000000000000000,
634  0.1000000000000000, 0.3000000000000000, 0.2000000000000000, 0.2000000000000000, 0.1000000000000000, 0.4000000000000000, 0.1000000000000000,
635  0.4000000000000000, 0.2000000000000000, 0.3000000000000000, 0.1000000000000000, 0.4000000000000000, 0.1000000000000000, 0.4000000000000000,
636  0.3000000000000000, 0.3000000000000000, 0.2000000000000000, 0.4000000000000000, 0.2000000000000000, 0.4000000000000000, 0.3000000000000000,
637  0.2000000000000000, 0.4000000000000000, 0.2000000000000000, 0.2000000000000000, 0.3000000000000000, 0.1000000000000000, 0.3000000000000000,
638  0.3000000000000000, 0.2000000000000000, 0.3000000000000000, 0.2000000000000000, 0.3000000000000000, 0.2000000000000000, 0.3000000000000000
639 };
640 static const double NC_TET_W84[] = {
641  0.2843915343915344, 0.2843915343915344, 0.2843915343915344, 0.2843915343915344, -0.3882275132275133, -0.3882275132275133, -0.3882275132275133,
642  -0.3882275132275133, -0.3882275132275133, -0.3882275132275133, -0.3882275132275133, -0.3882275132275133, -0.3882275132275133, -0.3882275132275133,
643  -0.3882275132275133, -0.3882275132275133, 0.8776455026455027, 0.8776455026455027, 0.8776455026455027, 0.8776455026455027, 0.8776455026455027,
644  0.8776455026455027, 0.8776455026455027, 0.8776455026455027, 0.8776455026455027, 0.8776455026455027, 0.8776455026455027, 0.8776455026455027,
645  0.1236772486772487, 0.1236772486772487, 0.1236772486772487, 0.1236772486772487, 0.1236772486772487, 0.1236772486772487, 0.1236772486772487,
646  0.1236772486772487, 0.1236772486772487, 0.1236772486772487, 0.1236772486772487, 0.1236772486772487, -0.8584656084656085, -0.8584656084656085,
647  -0.8584656084656085, -0.8584656084656085, -0.8584656084656085, -0.8584656084656085, -0.2632275132275133, -0.2632275132275133, -0.2632275132275133,
648  -0.2632275132275133, -0.2632275132275133, -0.2632275132275133, -0.2632275132275133, -0.2632275132275133, -0.2632275132275133, -0.2632275132275133,
649  -0.2632275132275133, -0.2632275132275133, -0.2632275132275133, -0.2632275132275133, -0.2632275132275133, -0.2632275132275133, -0.2632275132275133,
650  -0.2632275132275133, -0.2632275132275133, -0.2632275132275133, -0.2632275132275133, -0.2632275132275133, -0.2632275132275133, -0.2632275132275133,
651  0.0145502645502645, 0.0145502645502645, 0.0145502645502645, 0.0145502645502645, 1.0165343915343916, 1.0165343915343916, 1.0165343915343916,
652  1.0165343915343916, -0.0251322751322751, -0.0251322751322751, -0.0251322751322751, -0.0251322751322751, -0.0251322751322751, -0.0251322751322751
653 };
654 
655 #ifdef __cplusplus
656 }
657 #endif
658 
659 #endif
static const double G_TET_Y10[]
Definition: fem_tools.h:538
static const double G_TRI_Y7[]
Definition: fem_tools.h:282
PetscErrorCode ShapeMBTETQ_inverse(double *N, double *diffN, const double *elem_coords, const double *glob_coords, double *loc_coords, const double eps)
Definition: fem_tools.c:1019
PetscErrorCode ShapeInvJacVolume(double *jac)
Definition: fem_tools.c:51
static const double NC_TET_W84[]
Definition: fem_tools.h:640
static const double G_TET_X4[]
Definition: fem_tools.h:498
static const double G_TRI_W37[]
Definition: fem_tools.h:361
static const double G_TET_Z45[]
Definition: fem_tools.h:576
PetscErrorCode ShapeFaceBaseMBTRI(double *diffN, const double *coords, double *normal, double *s1, double *s2)
Definition: fem_tools.c:216
static const double G_TET_Y45[]
Definition: fem_tools.h:565
double ShapeVolumeMBTETQ(const double *diffN, const double *coords, int G_DIM, double *G_TET_W)
Definition: fem_tools.c:1005
static const double G_TRI_W28[]
Definition: fem_tools.h:335
PetscErrorCode Spin(double *spinOmega, double *vecOmega)
calculate spin matrix from vector
Definition: fem_tools.c:558
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:318
static const double G_TRI_X286[]
Definition: fem_tools.h:370
static const double G_TRI_Y286[]
Definition: fem_tools.h:408
static const double G_TRI_X37[]
Definition: fem_tools.h:343
static const double G_TRI_W286[]
Definition: fem_tools.h:446
static const double G_TRI_X1[]
Definition: fem_tools.h:247
PetscErrorCode ShapeMBTRIQ(double *N, const double *X, const double *Y, const int G_DIM)
Definition: fem_tools.c:792
PetscErrorCode ShapeFaceNormalMBTRI(double *diffN, const double *coords, double *normal)
Definition: fem_tools.c:241
static const double G_TET_Z10[]
Definition: fem_tools.h:543
void print_mat_complex(__CLPK_doublecomplex *M, int m, int n)
priint complex matrix
double ShapeVolumeMBTET(double *diffN, const double *coords)
calculate TET volume
Definition: fem_tools.c:310
static const double G_TET_X1[]
Definition: fem_tools.h:486
PetscErrorCode make_complex_matrix(double *reA, double *imA, __CLPK_doublecomplex *xA)
Compose complex matrix (3x3) from two real matrices.
Definition: fem_tools.c:569
static const double G_TRI_X28[]
Definition: fem_tools.h:321
static const double G_TRI_W3[]
Definition: fem_tools.h:262
CBLAS_TRANSPOSE
Definition: cblas.h:11
PetscErrorCode Grundmann_Moeller_integration_points_2D_TRI(int rule, double *G_TRI_X, double *G_TRI_Y, double *G_TRI_W)
Compute weights and integration points for 2D Triangle using Grundmann_Moeller rule.
Definition: fem_tools.c:110
PetscErrorCode MakeComplexTensor(double *reA, double *imA, __CLPK_doublecomplex *xA)
Definition: fem_tools.c:499
static const double G_TRI_Y3[]
Definition: fem_tools.h:259
static const double G_TET_Y1[]
Definition: fem_tools.h:489
void ShapeJacMBTRI(double *diffN, const double *coords, double *Jac)
calculate jacobioan
static const double G_TRI_W1[]
Definition: fem_tools.h:253
PetscErrorCode ShapeDiffMBTETinvJ(double *diffN, double *invJac, double *diffNinvJac)
calculate shape functions derivatives in space
Definition: fem_tools.c:427
static const double NC_TET_Z84[]
Definition: fem_tools.h:626
static const double G_TRI_W13[]
Definition: fem_tools.h:300
static const double G_TET_Z4[]
Definition: fem_tools.h:506
PetscErrorCode ShapeMBTET_inverse(double *N, double *diffN, const double *elem_coords, const double *glob_coords, double *loc_coords)
calculate local coordinates for given global coordinates
Definition: fem_tools.c:347
static const double NC_TET_X84[]
Definition: fem_tools.h:598
PetscErrorCode ShapeJacMBTETQ(const double *diffN, const double *coords, double *Jac)
Definition: fem_tools.c:979
PetscErrorCode ShapeMBTETQ_GAUSS(double *N, const double *X, const double *Y, const double *Z, const int G_DIM)
Definition: fem_tools.c:920
PetscErrorCode ShapeMBTETQ_detJac_at_Gauss_Points(double *detJac_at_Gauss_Points, const double *diffN, const double *coords, int G_DIM)
Definition: fem_tools.c:991
PetscErrorCode Base_scale(__CLPK_doublecomplex *xnormal, __CLPK_doublecomplex *xs1, __CLPK_doublecomplex *xs2)
Definition: fem_tools.c:753
PetscErrorCode GradientOfDeformation(double *diffN, double *dofs, double *F)
calculate gradient of deformation
Definition: fem_tools.c:437
static const double G_TRI_Y28[]
Definition: fem_tools.h:328
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:331
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:206
PetscErrorCode ShapeMBEDGE(double *N, const double *G_X, int DIM)
Definition: fem_tools.c:773
static const double G_TET_X10[]
Definition: fem_tools.h:533
static const double NC_TET_Y84[]
Definition: fem_tools.h:612
static const double G_TET_Z1[]
Definition: fem_tools.h:492
static const double G_TET_W1[]
Definition: fem_tools.h:495
PetscErrorCode Normal_hierarchical(int order_approx, int *order_edge_approx, int order, int *order_edge, double *diffN, double *diffN_face, double *diffN_edge[], double *dofs, double *dofs_edge[], double *dofs_face, double *idofs, double *idofs_edge[], double *idofs_face, __CLPK_doublecomplex *xnormal, __CLPK_doublecomplex *s1, __CLPK_doublecomplex *s2, int gg)
Complex normal.
Definition: fem_tools.c:581
static const double G_TRI_Y37[]
Definition: fem_tools.h:352
static const double G_TRI_X13[]
Definition: fem_tools.h:290
static const double G_TRI_W4[]
Definition: fem_tools.h:273
static const double G_TRI_X7[]
Definition: fem_tools.h:277
static const double G_TET_W5[]
Definition: fem_tools.h:528
static const double G_TET_X45[]
Definition: fem_tools.h:554
PetscErrorCode ShapeMBTETQ(double *N, const double x, const double y, const double z)
Definition: fem_tools.c:870
PetscErrorCode Grundmann_Moeller_integration_points_3D_TET(int rule, double *G_TET_X, double *G_TET_Y, double *G_TET_Z, double *G_TET_W)
Compute weights and integration points for 3D Tet using Grundmann_Moeller rule.
Definition: fem_tools.c:155
PetscErrorCode DeterminantComplexGradient(__CLPK_doublecomplex *xF, __CLPK_doublecomplex *det_xF)
Definition: fem_tools.c:538
PetscErrorCode ShapeDiffMBEDGE(double *diffN)
Definition: fem_tools.c:783
static const double G_TET_Z5[]
Definition: fem_tools.h:523
void ShapeDiffMBTRIinvJ(double *diffN, double *invJac, double *diffNinvJac)
calculate derivatives of shape functions in space
static const double G_TRI_X3[]
Definition: fem_tools.h:256
PetscErrorCode ShapeFaceDiffNormalMBTRI(double *diffN, const double *coords, double *diff_normal)
calculate derivative of normal in respect to nodal positions
Definition: fem_tools.c:249
void print_mat(double *M, int m, int n)
print matric M
PetscErrorCode ShapeDiffMBTRIQ(double *diffN, const double *X, const double *Y, const int G_DIM)
Definition: fem_tools.c:807
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition: fem_tools.c:194
void ShapeDiffMBTETinvJ_complex(double *diffN, __CLPK_doublecomplex *invJac, __CLPK_doublecomplex *diffNinvJac, const enum CBLAS_TRANSPOSE Trans)
Definition: fem_tools.c:449
static const double G_TRI_W19[]
Definition: fem_tools.h:315
static const double G_TET_Y5[]
Definition: fem_tools.h:518
static const double G_TRI_Y1[]
Definition: fem_tools.h:250
static const double G_TET_W4[]
Definition: fem_tools.h:510
PetscErrorCode ShapeDiffMBTETQ(double *diffN, const double x, const double y, const double z)
Definition: fem_tools.c:885
static const double G_TET_X5[]
Definition: fem_tools.h:513
void print_mat_sym_upper(double *M, int m, int n)
print upper part of the symmetric matrix
static const double G_TRI_W7[]
Definition: fem_tools.h:287
static const double eps
double ShapeDetJacVolume(double *jac)
determined of jacobian
Definition: fem_tools.c:34
PetscErrorCode ShapeJacMBTET(double *diffN, const double *coords, double *jac)
calculate jacobian
Definition: fem_tools.c:300
static const double G_TRI_Y19[]
Definition: fem_tools.h:310
static const double G_TRI_X19[]
Definition: fem_tools.h:305
static const double G_TET_W10[]
Definition: fem_tools.h:548
PetscErrorCode ShapeMBTRI_inverse(double *N, double *diffN, const double *elem_coords, const double *glob_coords, double *loc_coords)
calculate local coordinates of triangle element for given global coordinates in 2D (Assume e....
Definition: fem_tools.c:392
PetscErrorCode InvertComplexGradient(__CLPK_doublecomplex *xF)
Definition: fem_tools.c:511
static const double G_TRI_X4[]
Definition: fem_tools.h:265
static const double G_TET_Y4[]
Definition: fem_tools.h:502
const int N
Definition: speed_test.cpp:3
static const double G_TET_W45[]
Definition: fem_tools.h:587
static const double G_TRI_Y4[]
Definition: fem_tools.h:269
PetscErrorCode Grundmann_Moeller_integration_points_1D_EDGE(int rule, double *G_TRI_X, double *G_TRI_W)
Compute weights and integration points for edge using Grundmann_Moeller rule.
Definition: fem_tools.c:67
PetscErrorCode ShapeDiffMBTETQ_GAUSS(double *diffN, const double *X, const double *Y, const double *Z, const int G_DIM)
Definition: fem_tools.c:939
static const double G_TRI_Y13[]
Definition: fem_tools.h:295
PetscErrorCode ShapeFaceNormalMBTRI_complex(double *diffN, __CLPK_doublecomplex *xcoords, __CLPK_doublecomplex *xnormal)
Definition: fem_tools.c:464