v0.10.0
fem_tools.c
Go to the documentation of this file.
1 /** \file fem_tools.c
2  * \brief Loose implementation of some useful functions
3  */
4 
5 /* This file is part of MoFEM.
6  * MoFEM is free software: you can redistribute it and/or modify it under
7  * the terms of the GNU Lesser General Public License as published by the
8  * Free Software Foundation, either version 3 of the License, or (at your
9  * option) any later version.
10  *
11  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
12  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
14  * License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
18 
19 #include <definitions.h>
20 #include <stdlib.h>
21 #include <stdio.h>
22 #include <math.h>
23 #include <complex.h>
24 #include <assert.h>
25 #include <string.h>
26 
27 #include <fem_tools.h>
28 #include <gm_rule.h>
29 
30 #include <h1_hdiv_hcurl_l2.h>
31 
32 static PetscErrorCode ierr;
33 
34 double ShapeDetJacVolume(double *jac) {
35  double det_jac;
36  __CLPK_integer ipiv[4];
37  __CLPK_integer info = lapack_dgetrf(3, 3, jac, 3, ipiv);
38  if (info != 0)
39  return -1;
40  int i = 0, j = 0;
41  det_jac = 1.;
42  for (; i < 3; i++) {
43  det_jac *= jac[3 * i + i];
44  if (ipiv[i] != i + 1)
45  j++;
46  }
47  if ((j - (j / 2) * 2) != 0)
48  det_jac = -det_jac;
49  return det_jac;
50 }
51 PetscErrorCode ShapeInvJacVolume(double *jac) {
52  MoFEMFunctionBeginHot;
53  __CLPK_integer ipiv[4];
54  __CLPK_doublereal work[3];
55  __CLPK_integer lwork = 3;
56  __CLPK_integer info;
57  info = lapack_dgetrf(3, 3, jac, 3, ipiv);
58  if (info != 0)
59  SETERRQ1(PETSC_COMM_SELF, 1, "info = %d", info);
60  info = lapack_dgetri(3, jac, 3, ipiv, work, lwork);
61  if (info != 0)
62  SETERRQ1(PETSC_COMM_SELF, 1, "info = %d", info);
63  MoFEMFunctionReturnHot(0);
64 }
65 
66 // MBTRI
67 PetscErrorCode Grundmann_Moeller_integration_points_1D_EDGE(int rule,
68  double *G_TRI_X,
69  double *G_TRI_W) {
70  MoFEMFunctionBeginHot;
71 
72  int dim_num = 1;
73  int point;
74  int point_num;
75  double *w;
76  double *x;
77 
78  // GM_RULE_SET determines the weights and abscissas
79  // pof a Grundmann-Moeller quadrature rule for
80  // the DIM_NUM dimensional simplex,
81  // using a rule of in index RULE,
82  // which will have degree of exactness 2*RULE+1.
83 
84  // printf ( " Here we use DIM_NUM = %d\n", dim_num );
85  // printf ( " RULE = %d\n", rule );
86  // printf ( " DEGREE = %d\n", 2 * rule + 1 );
87 
88  point_num = gm_rule_size(rule, dim_num);
89 
90  ierr = PetscMalloc(point_num * sizeof(double), &w);
91  CHKERRQ(ierr);
92  ierr = PetscMalloc(dim_num * point_num * sizeof(double), &x);
93  CHKERRQ(ierr);
94 
95  gm_rule_set(rule, dim_num, point_num, w, x);
96 
97  for (point = 0; point < point_num; point++) {
98  G_TRI_X[point] = x[0 + point * dim_num];
99  G_TRI_W[point] = w[point];
100  }
101 
102  ierr = PetscFree(w);
103  CHKERRQ(ierr);
104  ierr = PetscFree(x);
105  CHKERRQ(ierr);
106 
107  MoFEMFunctionReturnHot(0);
108 }
109 
110 PetscErrorCode Grundmann_Moeller_integration_points_2D_TRI(int rule,
111  double *G_TRI_X,
112  double *G_TRI_Y,
113  double *G_TRI_W) {
114  MoFEMFunctionBeginHot;
115 
116  int dim_num = 2;
117  int point;
118  int point_num;
119  double *w;
120  double *x;
121 
122  // GM_RULE_SET determines the weights and abscissas
123  // pof a Grundmann-Moeller quadrature rule for
124  // the DIM_NUM dimensional simplex,
125  // using a rule of in index RULE,
126  // which will have degree of exactness 2*RULE+1.
127 
128  // printf ( " Here we use DIM_NUM = %d\n", dim_num );
129  // printf ( " RULE = %d\n", rule );
130  // printf ( " DEGREE = %d\n", 2 * rule + 1 );
131 
132  point_num = gm_rule_size(rule, dim_num);
133 
134  ierr = PetscMalloc(point_num * sizeof(double), &w);
135  CHKERRQ(ierr);
136  ierr = PetscMalloc(dim_num * point_num * sizeof(double), &x);
137  CHKERRQ(ierr);
138 
139  gm_rule_set(rule, dim_num, point_num, w, x);
140 
141  for (point = 0; point < point_num; point++) {
142  G_TRI_X[point] = x[0 + point * dim_num];
143  G_TRI_Y[point] = x[1 + point * dim_num];
144  G_TRI_W[point] = w[point];
145  }
146 
147  ierr = PetscFree(w);
148  CHKERRQ(ierr);
149  ierr = PetscFree(x);
150  CHKERRQ(ierr);
151 
152  MoFEMFunctionReturnHot(0);
153 }
154 
155 PetscErrorCode Grundmann_Moeller_integration_points_3D_TET(int rule,
156  double *G_TET_X,
157  double *G_TET_Y,
158  double *G_TET_Z,
159  double *G_TET_W) {
160  MoFEMFunctionBeginHot;
161 
162  int dim_num = 3;
163  int point;
164  int point_num;
165  double *w;
166  double *x;
167 
168  // printf ( " Here we use DIM_NUM = %d\n", dim_num );
169  // printf ( " RULE = %d\n", rule );
170  // printf ( " DEGREE = %d\n", 2 * rule + 1 );
171 
172  point_num = gm_rule_size(rule, dim_num);
173 
174  ierr = PetscMalloc(point_num * sizeof(double), &w);
175  CHKERRQ(ierr);
176  ierr = PetscMalloc(dim_num * point_num * sizeof(double), &x);
177  CHKERRQ(ierr);
178 
179  gm_rule_set(rule, dim_num, point_num, w, x);
180  for (point = 0; point < point_num; point++) {
181  G_TET_X[point] = x[0 + point * dim_num];
182  G_TET_Y[point] = x[1 + point * dim_num];
183  G_TET_Z[point] = x[2 + point * dim_num];
184  G_TET_W[point] = w[point];
185  }
186 
187  ierr = PetscFree(w);
188  CHKERRQ(ierr);
189  ierr = PetscFree(x);
190  CHKERRQ(ierr);
191 
192  MoFEMFunctionReturnHot(0);
193 }
194 PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y,
195  const int G_DIM) {
196  MoFEMFunctionBeginHot;
197  int ii = 0;
198  for (; ii < G_DIM; ii++) {
199  double x = X[ii], y = Y[ii];
200  N[3 * ii + 0] = N_MBTRI0(x, y);
201  N[3 * ii + 1] = N_MBTRI1(x, y);
202  N[3 * ii + 2] = N_MBTRI2(x, y);
203  }
204  MoFEMFunctionReturnHot(0);
205 }
206 PetscErrorCode ShapeDiffMBTRI(double *diffN) {
207  MoFEMFunctionBeginHot;
208  diffN[0] = diffN_MBTRI0x;
209  diffN[1] = diffN_MBTRI0y;
210  diffN[2] = diffN_MBTRI1x;
211  diffN[3] = diffN_MBTRI1y;
212  diffN[4] = diffN_MBTRI2x;
213  diffN[5] = diffN_MBTRI2y;
214  MoFEMFunctionReturnHot(0);
215 }
216 PetscErrorCode ShapeFaceBaseMBTRI(double *diffN, const double *coords,
217  double *normal, double *s1, double *s2) {
218  MoFEMFunctionBeginHot;
219 
220  double diffX_ksi[3];
221  double diffX_eta[3];
222  int ii = 0;
223  for (; ii < 3; ii++) {
224  diffX_ksi[ii] = cblas_ddot(3, &coords[ii], 3, &diffN[0], 2);
225  diffX_eta[ii] = cblas_ddot(3, &coords[ii], 3, &diffN[1], 2);
226  }
227  if (s1 != NULL) {
228  cblas_dcopy(3, diffX_ksi, 1, s1, 1);
229  }
230  if (s2 != NULL) {
231  cblas_dcopy(3, diffX_eta, 1, s2, 1);
232  }
233  double Spin_diffX_ksi[9];
234  ierr = Spin(Spin_diffX_ksi, diffX_ksi);
235  CHKERRQ(ierr);
236  cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1., Spin_diffX_ksi, 3,
237  diffX_eta, 1, 0., normal, 1);
238  MoFEMFunctionReturnHot(0);
239 }
240 
241 PetscErrorCode ShapeFaceNormalMBTRI(double *diffN, const double *coords,
242  double *normal) {
243  MoFEMFunctionBeginHot;
244  ierr = ShapeFaceBaseMBTRI(diffN, coords, normal, NULL, NULL);
245  CHKERRQ(ierr);
246  MoFEMFunctionReturnHot(0);
247 }
248 
249 PetscErrorCode ShapeFaceDiffNormalMBTRI(double *diffN, const double *coords,
250  double *diff_normal) {
251  MoFEMFunctionBeginHot;
252  // N = Spin(dX/dksi)*dX/deta = -Spin(dX/deta)*dX/dksi
253 
254  double diffX_ksi[3];
255  double diffX_eta[3];
256  int ii = 0;
257  for (; ii < 3; ii++) {
258  diffX_ksi[ii] = cblas_ddot(3, &coords[ii], 3, &diffN[0], 2);
259  diffX_eta[ii] = cblas_ddot(3, &coords[ii], 3, &diffN[1], 2);
260  }
261  double Spin_diffX_ksi[9];
262  ierr = Spin(Spin_diffX_ksi, diffX_ksi);
263  CHKERRQ(ierr);
264  double Spin_diffX_eta[9];
265  ierr = Spin(Spin_diffX_eta, diffX_eta);
266  CHKERRQ(ierr);
267  double B_ksi[3 * 9];
268  bzero(B_ksi, 3 * 9 * sizeof(double));
269  double B_eta[3 * 9];
270  bzero(B_eta, 3 * 9 * sizeof(double));
271  // B_ksi[] = [
272  // diffN[2*0+0], 0, 0, diffN[2*1+0], 0, 0, diffN[2*2+0], 0,
273  // 0
274  // 0, diffN[2*0+0], 0, 0, diffN[2*1+0], 0, 0, diffN[2*2+0],
275  // 0
276  // 0, 0, diffM[2*0+0], 0, 0, diffN[2*1+0], 0, 0,
277  // diffN[2*2+0]
278  //]
279  // B_eta[] = [
280  // diffN[2*0+1], 0, 0, diffN[2*1+1], 0, 0, diffN[2*2+1], 0,
281  // 0
282  // 0, diffN[2*0+1], 0, 0, diffN[2*1+1], 0, 0, diffN[2*2+1],
283  // 0
284  // 0, 0, diffM[2*0+1], 0, 0, diffN[2*1+1], 0, 0,
285  // diffN[2*2+1]
286  //]
287  ii = 0;
288  for (; ii < 3; ii++) {
289  cblas_dcopy(3, &diffN[0], 2, &B_ksi[ii * 9 + ii], 3);
290  cblas_dcopy(3, &diffN[1], 2, &B_eta[ii * 9 + ii], 3);
291  }
292  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 3, 9, 3, +1.,
293  Spin_diffX_ksi, 3, B_eta, 9, 0., diff_normal, 9);
294  cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 3, 9, 3, -1.,
295  Spin_diffX_eta, 3, B_ksi, 9, 1., diff_normal, 9);
296  MoFEMFunctionReturnHot(0);
297 }
298 
299 // MBTET
300 PetscErrorCode ShapeJacMBTET(double *diffN, const double *coords, double *jac) {
301  MoFEMFunctionBeginHot;
302  int ii, jj, kk;
303  bzero(jac, sizeof(double) * 9);
304  for (ii = 0; ii < 4; ii++) // shape func.
305  for (jj = 0; jj < 3; jj++) // space
306  for (kk = 0; kk < 3; kk++) // derivative of shape func.
307  jac[jj * 3 + kk] += diffN[ii * 3 + kk] * coords[ii * 3 + jj];
308  MoFEMFunctionReturnHot(0);
309 }
310 double ShapeVolumeMBTET(double *diffN, const double *coords) {
311  double Jac[9];
312  ShapeJacMBTET(diffN, coords, Jac);
313  double detJac = ShapeDetJacVolume(Jac);
314  // printf("detJac = +%6.4e\n",detJac);
315  // print_mat(Jac,3,3);
316  return detJac * G_TET_W1[0] / 6.;
317 }
318 PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y,
319  const double *G_Z, int DIM) {
320  MoFEMFunctionBeginHot;
321  int ii = 0;
322  for (; ii < DIM; ii++) {
323  double x = G_X[ii], y = G_Y[ii], z = G_Z[ii];
324  N[4 * ii + 0] = N_MBTET0(x, y, z);
325  N[4 * ii + 1] = N_MBTET1(x, y, z);
326  N[4 * ii + 2] = N_MBTET2(x, y, z);
327  N[4 * ii + 3] = N_MBTET3(x, y, z);
328  }
329  MoFEMFunctionReturnHot(0);
330 }
331 PetscErrorCode ShapeDiffMBTET(double *diffN) {
332  MoFEMFunctionBeginHot;
333  diffN[0] = diffN_MBTET0x;
334  diffN[1] = diffN_MBTET0y;
335  diffN[2] = diffN_MBTET0z;
336  diffN[3] = diffN_MBTET1x;
337  diffN[4] = diffN_MBTET1y;
338  diffN[5] = diffN_MBTET1z;
339  diffN[6] = diffN_MBTET2x;
340  diffN[7] = diffN_MBTET2y;
341  diffN[8] = diffN_MBTET2z;
342  diffN[9] = diffN_MBTET3x;
343  diffN[10] = diffN_MBTET3y;
344  diffN[11] = diffN_MBTET3z;
345  MoFEMFunctionReturnHot(0);
346 }
347 PetscErrorCode ShapeMBTET_inverse(double *N, double *diffN,
348  const double *elem_coords,
349  const double *glob_coords,
350  double *loc_coords) {
351  MoFEMFunctionBeginHot;
352  double A[3 * 3];
353  int IPIV[3];
354  // COL MAJOR
355  // X
356  A[0 + 3 * 0] =
357  cblas_ddot(4, &diffN[0 * 3 + 0], 3, &elem_coords[0 * 3 + 0], 3);
358  A[0 + 3 * 1] =
359  cblas_ddot(4, &diffN[0 * 3 + 1], 3, &elem_coords[0 * 3 + 0], 3);
360  A[0 + 3 * 2] =
361  cblas_ddot(4, &diffN[0 * 3 + 2], 3, &elem_coords[0 * 3 + 0], 3);
362  loc_coords[0] =
363  glob_coords[0] - cblas_ddot(4, &N[0], 1, &elem_coords[0 * 3 + 0], 3);
364  // printf("A\n[ %3.2f %3.2f %3.2f ] %3.2f \n",A[0*3],A[1*3],A[2*3],R[0]);
365  // Y
366  A[1 + 3 * 0] =
367  cblas_ddot(4, &diffN[0 * 3 + 0], 3, &elem_coords[0 * 3 + 1], 3);
368  A[1 + 3 * 1] =
369  cblas_ddot(4, &diffN[0 * 3 + 1], 3, &elem_coords[0 * 3 + 1], 3);
370  A[1 + 3 * 2] =
371  cblas_ddot(4, &diffN[0 * 3 + 2], 3, &elem_coords[0 * 3 + 1], 3);
372  loc_coords[1] =
373  glob_coords[1] - cblas_ddot(4, &N[0], 1, &elem_coords[0 * 3 + 1], 3);
374  // printf("[ %3.2f %3.2f %3.2f ] %3.2f \n",A[1+3*0],A[1+3*1],A[1+3*2],R[1]);
375  // Z
376  A[2 + 3 * 0] =
377  cblas_ddot(4, &diffN[0 * 3 + 0], 3, &elem_coords[0 * 3 + 2], 3);
378  A[2 + 3 * 1] =
379  cblas_ddot(4, &diffN[0 * 3 + 1], 3, &elem_coords[0 * 3 + 2], 3);
380  A[2 + 3 * 2] =
381  cblas_ddot(4, &diffN[0 * 3 + 2], 3, &elem_coords[0 * 3 + 2], 3);
382  loc_coords[2] =
383  glob_coords[2] - cblas_ddot(4, &N[0], 1, &elem_coords[0 * 3 + 2], 3);
384  // printf("[ %3.2f %3.2f %3.2f ] %3.2f \n",A[2+3*0],A[2+3*1],A[2+3*2],R[1]);
385  int info =
386  lapack_dgesv(3, 1, &A[0], 3, (__CLPK_integer *)IPIV, loc_coords, 3);
387  if (info != 0)
388  SETERRQ1(PETSC_COMM_SELF, 1, "info == %d", info);
389  MoFEMFunctionReturnHot(0);
390 }
391 
392 PetscErrorCode ShapeMBTRI_inverse(double *N, double *diffN,
393  const double *elem_coords,
394  const double *glob_coords,
395  double *loc_coords) {
396  MoFEMFunctionBeginHot;
397  double A[2 * 2];
398 
399  // 1st and 2nd element of matrix A
400  A[0] = cblas_ddot(3, &diffN[0], 2, &elem_coords[0], 2); // dot product
401  A[1] = cblas_ddot(3, &diffN[1], 2, &elem_coords[0], 2);
402  loc_coords[0] = glob_coords[0] - cblas_ddot(3, &N[0], 1, &elem_coords[0], 2);
403 
404  // 3rd and 4th element of matrix A
405  A[2] = cblas_ddot(3, &diffN[0], 2, &elem_coords[1], 2);
406  A[3] = cblas_ddot(3, &diffN[1], 2, &elem_coords[1], 2);
407  loc_coords[1] = glob_coords[1] - cblas_ddot(3, &N[0], 1, &elem_coords[1], 2);
408 
409  // calculate directly the solution (as the size of matrix is only 2x2)
410  double invA[2 * 2], detA;
411  detA = A[0] * A[3] - A[1] * A[2];
412  detA = 1.0 / detA;
413  invA[0] = A[3] * detA;
414  invA[1] = -1.0 * A[1] * detA;
415  invA[2] = -1.0 * A[2] * detA;
416  invA[3] = A[0] * detA;
417 
418  double loc_coords_new[2];
419  loc_coords_new[0] = invA[0] * loc_coords[0] + invA[1] * loc_coords[1];
420  loc_coords_new[1] = invA[2] * loc_coords[0] + invA[3] * loc_coords[1];
421 
422  loc_coords[0] = loc_coords_new[0];
423  loc_coords[1] = loc_coords_new[1];
424  MoFEMFunctionReturnHot(0);
425 }
426 
427 PetscErrorCode ShapeDiffMBTETinvJ(double *diffN, double *invJac,
428  double *diffNinvJac) {
429  MoFEMFunctionBeginHot;
430  int ii = 0;
431  for (; ii < 4; ii++) {
432  cblas_dgemv(CblasRowMajor, CblasTrans, 3, 3, 1., invJac, 3, &diffN[ii * 3],
433  1, 0., &diffNinvJac[ii * 3], 1);
434  }
435  MoFEMFunctionReturnHot(0);
436 }
437 PetscErrorCode GradientOfDeformation(double *diffN, double *dofs, double *F) {
438  MoFEMFunctionBeginHot;
439  int col, row = 0;
440  for (; row < 3; row++)
441  for (col = 0; col < 3; col++) {
442  F[3 * row + col] = cblas_ddot(4, &diffN[col], 3, &dofs[row], 3);
443  }
444  MoFEMFunctionReturnHot(0);
445 }
446 
447 // Come functions with complex variables if one like to calculate derivative
448 // using complex variable
449 void ShapeDiffMBTETinvJ_complex(double *diffN, __CLPK_doublecomplex *invJac,
450  __CLPK_doublecomplex *diffNinvJac,
451  const enum CBLAS_TRANSPOSE Trans) {
452  __CLPK_doublecomplex tmp1 = {1., 0.}, tmp2 = {0., 0.};
453  int ii = 0, jj;
454  for (; ii < 4; ii++) {
455  __CLPK_doublecomplex tmp3[3];
456  for (jj = 0; jj < 3; jj++) {
457  tmp3[jj].r = diffN[ii * 3 + jj];
458  tmp3[jj].i = 0;
459  }
460  cblas_zgemv(CblasRowMajor, Trans, 3, 3, &tmp1, invJac, 3, tmp3, 1, &tmp2,
461  &diffNinvJac[ii * 3], 1);
462  }
463 }
464 PetscErrorCode ShapeFaceNormalMBTRI_complex(double *diffN,
465  __CLPK_doublecomplex *xcoords,
466  __CLPK_doublecomplex *xnormal) {
467  MoFEMFunctionBeginHot;
468  double complex diffX_x, diffX_y, diffX_z;
469  double complex diffY_x, diffY_y, diffY_z;
470  diffX_x = diffX_y = diffX_z = 0.;
471  diffY_x = diffY_y = diffY_z = 0.;
472  int ii;
473  for (ii = 0; ii < 3; ii++) {
474  diffX_x +=
475  (xcoords[3 * ii + 0].r + I * xcoords[3 * ii + 0].i) * diffN[2 * ii + 0];
476  diffX_y +=
477  (xcoords[3 * ii + 1].r + I * xcoords[3 * ii + 1].i) * diffN[2 * ii + 0];
478  diffX_z +=
479  (xcoords[3 * ii + 2].r + I * xcoords[3 * ii + 2].i) * diffN[2 * ii + 0];
480  diffY_x +=
481  (xcoords[3 * ii + 0].r + I * xcoords[3 * ii + 0].i) * diffN[2 * ii + 1];
482  diffY_y +=
483  (xcoords[3 * ii + 1].r + I * xcoords[3 * ii + 1].i) * diffN[2 * ii + 1];
484  diffY_z +=
485  (xcoords[3 * ii + 2].r + I * xcoords[3 * ii + 2].i) * diffN[2 * ii + 1];
486  }
487  double complex tmp;
488  tmp = diffX_y * diffY_z - diffX_z * diffY_y;
489  xnormal[0].r = creal(tmp);
490  xnormal[0].i = cimag(tmp);
491  tmp = diffX_z * diffY_x - diffX_x * diffY_z;
492  xnormal[1].r = creal(tmp);
493  xnormal[1].i = cimag(tmp);
494  tmp = diffX_x * diffY_y - diffX_y * diffY_x;
495  xnormal[2].r = creal(tmp);
496  xnormal[2].i = cimag(tmp);
497  MoFEMFunctionReturnHot(0);
498 }
499 PetscErrorCode MakeComplexTensor(double *reA, double *imA,
500  __CLPK_doublecomplex *xA) {
501  MoFEMFunctionBeginHot;
502  int ii = 0, jj;
503  for (; ii < 3; ii++) {
504  for (jj = 0; jj < 3; jj++) {
505  xA[3 * ii + jj].r = reA[3 * ii + jj];
506  xA[3 * ii + jj].i = imA[3 * ii + jj];
507  }
508  }
509  MoFEMFunctionReturnHot(0);
510 }
511 PetscErrorCode InvertComplexGradient(__CLPK_doublecomplex *xF) {
512  MoFEMFunctionBeginHot;
513  __CLPK_integer IPIV[4];
514  __CLPK_doublecomplex WORK[4];
515  __CLPK_integer LWORK = 4;
516  __CLPK_integer info;
517  info = lapack_zgetrf(3, 3, xF, 3, IPIV);
518  if (info != 0)
519  SETERRQ(PETSC_COMM_SELF, 1, "info == 0");
520  info = lapack_zgetri(3, xF, 3, IPIV, WORK, LWORK);
521  if (info != 0)
522  SETERRQ(PETSC_COMM_SELF, 1, "info == 0");
523  MoFEMFunctionReturnHot(0);
524 }
525 PetscErrorCode InvertComplexSymmMatrix3by3(__CLPK_doublecomplex *xC) {
526  MoFEMFunctionBeginHot;
527  __CLPK_integer info;
528  info = lapack_zpotrf('L', 3, xC, 3);
529  if (info == 0)
530  SETERRQ(PETSC_COMM_SELF, 1, "info == 0");
531  // assert(info == 0);
532  info = lapack_zpotri('L', 3, xC, 3);
533  if (info == 0)
534  SETERRQ(PETSC_COMM_SELF, 1, "info == 0");
535  // assert(info == 0);
536  MoFEMFunctionReturnHot(0);
537 }
538 PetscErrorCode DeterminantComplexGradient(__CLPK_doublecomplex *xF,
539  __CLPK_doublecomplex *det_xF) {
540  MoFEMFunctionBeginHot;
541  __CLPK_integer IPIV[4];
542  if (lapack_zgetrf(3, 3, xF, 3, IPIV) != 0) {
543  SETERRQ(PETSC_COMM_SELF, 1, "lapack_zgetrf(3,3,xF,3,IPIV) != 0");
544  }
545  double complex det = 1;
546  int i = 0, j = 0;
547  for (; i < 3; i++) {
548  det *= xF[3 * i + i].r + I * xF[3 * i + i].i;
549  if (IPIV[i] != i + 1)
550  j++;
551  }
552  if ((j - (j / 2) * 2) != 0)
553  det = -det;
554  (*det_xF).r = creal(det);
555  (*det_xF).i = cimag(det);
556  MoFEMFunctionReturnHot(0);
557 }
558 PetscErrorCode Spin(double *spinOmega, double *vecOmega) {
559  MoFEMFunctionBeginHot;
560  bzero(spinOmega, 9 * sizeof(double));
561  spinOmega[0 * 3 + 1] = -vecOmega[2];
562  spinOmega[0 * 3 + 2] = +vecOmega[1];
563  spinOmega[1 * 3 + 0] = +vecOmega[2];
564  spinOmega[1 * 3 + 2] = -vecOmega[0];
565  spinOmega[2 * 3 + 0] = -vecOmega[1];
566  spinOmega[2 * 3 + 1] = +vecOmega[0];
567  MoFEMFunctionReturnHot(0);
568 }
569 PetscErrorCode make_complex_matrix(double *reA, double *imA,
570  __CLPK_doublecomplex *xA) {
571  MoFEMFunctionBeginHot;
572  int ii = 0, jj;
573  for (; ii < 3; ii++) {
574  for (jj = 0; jj < 3; jj++) {
575  xA[3 * ii + jj].r = reA[3 * ii + jj];
576  xA[3 * ii + jj].i = imA[3 * ii + jj];
577  }
578  }
579  MoFEMFunctionReturnHot(0);
580 }
581 PetscErrorCode Normal_hierarchical(
582  int order_approx, int *order_edge_approx, int order, int *order_edge,
583  double *diffN, double *diffN_face, double *diffN_edge[], double *dofs,
584  double *dofs_edge[], double *dofs_face, double *idofs, double *idofs_edge[],
585  double *idofs_face, __CLPK_doublecomplex *xnormal,
586  __CLPK_doublecomplex *xs1, __CLPK_doublecomplex *xs2, int gg) {
587  MoFEMFunctionBeginHot;
588  int nn, ee, dd;
589  // node
590  double complex diffX_x_node, diffX_y_node, diffX_z_node;
591  double complex diffY_x_node, diffY_y_node, diffY_z_node;
592  diffX_x_node = 0.;
593  diffX_y_node = 0.;
594  diffX_z_node = 0.;
595  diffY_x_node = 0.;
596  diffY_y_node = 0.;
597  diffY_z_node = 0.;
598  if (dofs != NULL || idofs != NULL) {
599  nn = 0;
600  for (; nn < 3; nn++) {
601  if (dofs != NULL) {
602  diffX_x_node += dofs[3 * nn + 0] * diffN[2 * nn + 0];
603  diffX_y_node += dofs[3 * nn + 1] * diffN[2 * nn + 0];
604  diffX_z_node += dofs[3 * nn + 2] * diffN[2 * nn + 0];
605  diffY_x_node += dofs[3 * nn + 0] * diffN[2 * nn + 1];
606  diffY_y_node += dofs[3 * nn + 1] * diffN[2 * nn + 1];
607  diffY_z_node += dofs[3 * nn + 2] * diffN[2 * nn + 1];
608  }
609  if (idofs != NULL) {
610  diffX_x_node += I * idofs[3 * nn + 0] * diffN[2 * nn + 0];
611  diffX_y_node += I * idofs[3 * nn + 1] * diffN[2 * nn + 0];
612  diffX_z_node += I * idofs[3 * nn + 2] * diffN[2 * nn + 0];
613  diffY_x_node += I * idofs[3 * nn + 0] * diffN[2 * nn + 1];
614  diffY_y_node += I * idofs[3 * nn + 1] * diffN[2 * nn + 1];
615  diffY_z_node += I * idofs[3 * nn + 2] * diffN[2 * nn + 1];
616  }
617  }
618  }
619  double complex diffX_x, diffX_y, diffX_z;
620  double complex diffY_x, diffY_y, diffY_z;
621  diffX_x = diffX_x_node;
622  diffX_y = diffX_y_node;
623  diffX_z = diffX_z_node;
624  diffY_x = diffY_x_node;
625  diffY_y = diffY_y_node;
626  diffY_z = diffY_z_node;
627  if (dofs_face != NULL || idofs_face != NULL) {
628  int nb_dofs_face = NBFACETRI_H1(order);
629  int nb_dofs_approx_face = NBFACETRI_H1(order_approx);
630  if (nb_dofs_face > 0) {
631  if (dofs_face != NULL) {
632  diffX_x += cblas_ddot(nb_dofs_face, &dofs_face[0], 3,
633  &diffN_face[gg * 2 * nb_dofs_approx_face + 0], 2);
634  diffX_y += cblas_ddot(nb_dofs_face, &dofs_face[1], 3,
635  &diffN_face[gg * 2 * nb_dofs_approx_face + 0], 2);
636  diffX_z += cblas_ddot(nb_dofs_face, &dofs_face[2], 3,
637  &diffN_face[gg * 2 * nb_dofs_approx_face + 0], 2);
638  diffY_x += cblas_ddot(nb_dofs_face, &dofs_face[0], 3,
639  &diffN_face[gg * 2 * nb_dofs_approx_face + 1], 2);
640  diffY_y += cblas_ddot(nb_dofs_face, &dofs_face[1], 3,
641  &diffN_face[gg * 2 * nb_dofs_approx_face + 1], 2);
642  diffY_z += cblas_ddot(nb_dofs_face, &dofs_face[2], 3,
643  &diffN_face[gg * 2 * nb_dofs_approx_face + 1], 2);
644  }
645  if (idofs_face != NULL) {
646  diffX_x +=
647  I * cblas_ddot(nb_dofs_face, &idofs_face[0], 3,
648  &diffN_face[gg * 2 * nb_dofs_approx_face + 0], 2);
649  diffX_y +=
650  I * cblas_ddot(nb_dofs_face, &idofs_face[1], 3,
651  &diffN_face[gg * 2 * nb_dofs_approx_face + 0], 2);
652  diffX_z +=
653  I * cblas_ddot(nb_dofs_face, &idofs_face[2], 3,
654  &diffN_face[gg * 2 * nb_dofs_approx_face + 0], 2);
655  diffY_x +=
656  I * cblas_ddot(nb_dofs_face, &idofs_face[0], 3,
657  &diffN_face[gg * 2 * nb_dofs_approx_face + 1], 2);
658  diffY_y +=
659  I * cblas_ddot(nb_dofs_face, &idofs_face[1], 3,
660  &diffN_face[gg * 2 * nb_dofs_approx_face + 1], 2);
661  diffY_z +=
662  I * cblas_ddot(nb_dofs_face, &idofs_face[2], 3,
663  &diffN_face[gg * 2 * nb_dofs_approx_face + 1], 2);
664  }
665  }
666  }
667  ee = 0;
668  if (dofs_edge != NULL || idofs_edge != NULL) {
669  for (; ee < 3; ee++) {
670  int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
671  int nb_dofs_approx_edge = NBEDGE_H1(order_edge_approx[ee]);
672  if (nb_dofs_edge > 0) {
673  if (dofs_edge != NULL) {
674  if (dofs_edge[ee] != NULL) {
675  diffX_x += cblas_ddot(
676  nb_dofs_edge, &(dofs_edge[ee])[0], 3,
677  &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 0], 2);
678  diffX_y += cblas_ddot(
679  nb_dofs_edge, &(dofs_edge[ee])[1], 3,
680  &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 0], 2);
681  diffX_z += cblas_ddot(
682  nb_dofs_edge, &(dofs_edge[ee])[2], 3,
683  &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 0], 2);
684  diffY_x += cblas_ddot(
685  nb_dofs_edge, &(dofs_edge[ee])[0], 3,
686  &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 1], 2);
687  diffY_y += cblas_ddot(
688  nb_dofs_edge, &(dofs_edge[ee])[1], 3,
689  &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 1], 2);
690  diffY_z += cblas_ddot(
691  nb_dofs_edge, &(dofs_edge[ee])[2], 3,
692  &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 1], 2);
693  }
694  }
695  if (idofs_edge != NULL) {
696  if (idofs_edge[ee] == NULL)
697  continue;
698  diffX_x +=
699  I * cblas_ddot(
700  nb_dofs_edge, &(idofs_edge[ee])[0], 3,
701  &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 0], 2);
702  diffX_y +=
703  I * cblas_ddot(
704  nb_dofs_edge, &(idofs_edge[ee])[1], 3,
705  &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 0], 2);
706  diffX_z +=
707  I * cblas_ddot(
708  nb_dofs_edge, &(idofs_edge[ee])[2], 3,
709  &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 0], 2);
710  diffY_x +=
711  I * cblas_ddot(
712  nb_dofs_edge, &(idofs_edge[ee])[0], 3,
713  &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 1], 2);
714  diffY_y +=
715  I * cblas_ddot(
716  nb_dofs_edge, &(idofs_edge[ee])[1], 3,
717  &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 1], 2);
718  diffY_z +=
719  I * cblas_ddot(
720  nb_dofs_edge, &(idofs_edge[ee])[2], 3,
721  &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 1], 2);
722  }
723  }
724  }
725  }
726  double complex normal[3];
727  normal[0] = diffX_y * diffY_z - diffX_z * diffY_y;
728  normal[1] = diffX_z * diffY_x - diffX_x * diffY_z;
729  normal[2] = diffX_x * diffY_y - diffX_y * diffY_x;
730  dd = 0;
731  for (; dd < 3; dd++) {
732  xnormal[dd].r = creal(normal[dd]);
733  xnormal[dd].i = cimag(normal[dd]);
734  }
735  if (xs1 != NULL) {
736  xs1[0].r = creal(diffX_x);
737  xs1[0].i = cimag(diffX_x);
738  xs1[1].r = creal(diffX_y);
739  xs1[1].i = cimag(diffX_y);
740  xs1[2].r = creal(diffX_z);
741  xs1[2].i = cimag(diffX_z);
742  }
743  if (xs2 != NULL) {
744  xs2[0].r = creal(diffY_x);
745  xs2[0].i = cimag(diffY_x);
746  xs2[1].r = creal(diffY_y);
747  xs2[1].i = cimag(diffY_y);
748  xs2[2].r = creal(diffY_z);
749  xs2[2].i = cimag(diffY_z);
750  }
751  MoFEMFunctionReturnHot(0);
752 }
753 PetscErrorCode Base_scale(__CLPK_doublecomplex *xnormal,
754  __CLPK_doublecomplex *xs1,
755  __CLPK_doublecomplex *xs2) {
756  MoFEMFunctionBeginHot;
757  complex double xnrm2_normal = csqrt(cpow(xnormal[0].r + I * xnormal[0].i, 2) +
758  cpow(xnormal[1].r + I * xnormal[1].i, 2) +
759  cpow(xnormal[2].r + I * xnormal[2].i, 2));
760  int dd = 0;
761  for (; dd < 3; dd++) {
762  complex double s1 = (xs1[dd].r + I * xs1[dd].i) * xnrm2_normal;
763  complex double s2 = (xs2[dd].r + I * xs2[dd].i) * xnrm2_normal;
764  xs1[dd].r = creal(s1);
765  xs1[dd].i = cimag(s1);
766  xs2[dd].r = creal(s2);
767  xs2[dd].i = cimag(s2);
768  }
769  MoFEMFunctionReturnHot(0);
770 }
771 
772 // MBEDGE
773 PetscErrorCode ShapeMBEDGE(double *N, const double *G_X, int DIM) {
774  MoFEMFunctionBeginHot;
775  int ii = 0;
776  for (; ii < DIM; ii++) {
777  double x = G_X[ii];
778  N[2 * ii + 0] = N_MBEDGE0(x);
779  N[2 * ii + 1] = N_MBEDGE1(x);
780  }
781  MoFEMFunctionReturnHot(0);
782 }
783 PetscErrorCode ShapeDiffMBEDGE(double *diffN) {
784  MoFEMFunctionBeginHot;
785  diffN[0] = diffN_MBEDGE0;
786  diffN[1] = diffN_MBEDGE1;
787  MoFEMFunctionReturnHot(0);
788 }
789 
790 // FIXME: NOT PROPERLY TESTED YET
791 // HO
792 PetscErrorCode ShapeMBTRIQ(double *N, const double *X, const double *Y,
793  const int G_DIM) {
794  MoFEMFunctionBeginHot;
795  int ii = 0;
796  for (; ii < G_DIM; ii++) {
797  double x = X[ii], y = Y[ii];
798  N[6 * ii + 0] = N_MBTRIQ0(x, y);
799  N[6 * ii + 1] = N_MBTRIQ1(x, y);
800  N[6 * ii + 2] = N_MBTRIQ2(x, y);
801  N[6 * ii + 3] = N_MBTRIQ3(x, y);
802  N[6 * ii + 4] = N_MBTRIQ4(x, y);
803  N[6 * ii + 5] = N_MBTRIQ5(x, y);
804  }
805  MoFEMFunctionReturnHot(0);
806 }
807 PetscErrorCode ShapeDiffMBTRIQ(double *diffN, const double *X, const double *Y,
808  const int G_DIM) {
809  MoFEMFunctionBeginHot;
810  int ii = 0;
811  for (; ii < G_DIM; ii++) {
812  double x = X[ii], y = Y[ii];
813  diffN[12 * ii + 0] = diffN_MBTRIQ0x(x, y);
814  diffN[12 * ii + 1] = diffN_MBTRIQ0y(x, y);
815  diffN[12 * ii + 2] = diffN_MBTRIQ1x(x, y);
816  diffN[12 * ii + 3] = diffN_MBTRIQ1y(x, y);
817  diffN[12 * ii + 4] = diffN_MBTRIQ2x(x, y);
818  diffN[12 * ii + 5] = diffN_MBTRIQ2y(x, y);
819  diffN[12 * ii + 6] = diffN_MBTRIQ3x(x, y);
820  diffN[12 * ii + 7] = diffN_MBTRIQ3y(x, y);
821  diffN[12 * ii + 8] = diffN_MBTRIQ4x(x, y);
822  diffN[12 * ii + 9] = diffN_MBTRIQ4y(x, y);
823  diffN[12 * ii + 10] = diffN_MBTRIQ5x(x, y);
824  diffN[12 * ii + 11] = diffN_MBTRIQ5y(x, y);
825  }
826  MoFEMFunctionReturnHot(0);
827 }
828 
829 // MBTETQ (JULIEN WORK)
830 #define N_MBTETQ0(x, y, z) ((2. * (1. - x - y - z) - 1.) * (1. - x - y - z))
831 #define N_MBTETQ1(x, y, z) ((2. * x - 1.) * x)
832 #define N_MBTETQ2(x, y, z) ((2. * y - 1.) * y)
833 #define N_MBTETQ3(x, y, z) ((2. * z - 1.) * z)
834 #define N_MBTETQ4(x, y, z) (4. * (1. - x - y - z) * x)
835 #define N_MBTETQ5(x, y, z) (4. * x * y)
836 #define N_MBTETQ6(x, y, z) (4. * (1. - x - y - z) * y)
837 #define N_MBTETQ7(x, y, z) (4. * (1. - x - y - z) * z)
838 #define N_MBTETQ8(x, y, z) (4. * x * z)
839 #define N_MBTETQ9(x, y, z) (4. * y * z)
840 #define diffN_MBTETQ0x(x, y, z) (-3. + 4. * x + 4. * y + 4. * z)
841 #define diffN_MBTETQ0y(x, y, z) (-3. + 4. * x + 4. * y + 4. * z)
842 #define diffN_MBTETQ0z(x, y, z) (-3. + 4. * x + 4. * y + 4. * z)
843 #define diffN_MBTETQ1x(x, y, z) (4. * x - 1.)
844 #define diffN_MBTETQ1y(x, y, z) (0.)
845 #define diffN_MBTETQ1z(x, y, z) (0.)
846 #define diffN_MBTETQ2x(x, y, z) (0.)
847 #define diffN_MBTETQ2y(x, y, z) (4. * y - 1.)
848 #define diffN_MBTETQ2z(x, y, z) (0.)
849 #define diffN_MBTETQ3x(x, y, z) (0.)
850 #define diffN_MBTETQ3y(x, y, z) (0.)
851 #define diffN_MBTETQ3z(x, y, z) (4. * z - 1.)
852 #define diffN_MBTETQ4x(x, y, z) (-8. * x + 4. - 4. * y - 4. * z)
853 #define diffN_MBTETQ4y(x, y, z) (-4. * x)
854 #define diffN_MBTETQ4z(x, y, z) (-4. * x)
855 #define diffN_MBTETQ5x(x, y, z) (4. * y)
856 #define diffN_MBTETQ5y(x, y, z) (4. * x)
857 #define diffN_MBTETQ5z(x, y, z) (0.)
858 #define diffN_MBTETQ6x(x, y, z) (-4. * y)
859 #define diffN_MBTETQ6y(x, y, z) (-8. * y + 4. - 4. * x - 4. * z)
860 #define diffN_MBTETQ6z(x, y, z) (-4. * y)
861 #define diffN_MBTETQ7x(x, y, z) (-4. * z)
862 #define diffN_MBTETQ7y(x, y, z) (-4. * z)
863 #define diffN_MBTETQ7z(x, y, z) (-8. * z + 4. - 4. * x - 4. * y)
864 #define diffN_MBTETQ8x(x, y, z) (4. * z)
865 #define diffN_MBTETQ8y(x, y, z) (0.)
866 #define diffN_MBTETQ8z(x, y, z) (4. * x)
867 #define diffN_MBTETQ9x(x, y, z) (0.)
868 #define diffN_MBTETQ9y(x, y, z) (4. * z)
869 #define diffN_MBTETQ9z(x, y, z) (4. * y)
870 PetscErrorCode ShapeMBTETQ(double *N, const double x, const double y,
871  const double z) {
872  MoFEMFunctionBeginHot;
873  N[0] = N_MBTETQ0(x, y, z);
874  N[1] = N_MBTETQ1(x, y, z);
875  N[2] = N_MBTETQ2(x, y, z);
876  N[3] = N_MBTETQ3(x, y, z);
877  N[4] = N_MBTETQ4(x, y, z);
878  N[5] = N_MBTETQ5(x, y, z);
879  N[6] = N_MBTETQ6(x, y, z);
880  N[7] = N_MBTETQ7(x, y, z);
881  N[8] = N_MBTETQ8(x, y, z);
882  N[9] = N_MBTETQ9(x, y, z);
883  MoFEMFunctionReturnHot(0);
884 }
885 PetscErrorCode ShapeDiffMBTETQ(double *diffN, const double x, const double y,
886  const double z) {
887  MoFEMFunctionBeginHot;
888  diffN[0] = diffN_MBTETQ0x(x, y, z);
889  diffN[1] = diffN_MBTETQ0y(x, y, z);
890  diffN[2] = diffN_MBTETQ0z(x, y, z);
891  diffN[3] = diffN_MBTETQ1x(x, y, z);
892  diffN[4] = diffN_MBTETQ1y(x, y, z);
893  diffN[5] = diffN_MBTETQ1z(x, y, z);
894  diffN[6] = diffN_MBTETQ2x(x, y, z);
895  diffN[7] = diffN_MBTETQ2y(x, y, z);
896  diffN[8] = diffN_MBTETQ2z(x, y, z);
897  diffN[9] = diffN_MBTETQ3x(x, y, z);
898  diffN[10] = diffN_MBTETQ3y(x, y, z);
899  diffN[11] = diffN_MBTETQ3z(x, y, z);
900  diffN[12] = diffN_MBTETQ4x(x, y, z);
901  diffN[13] = diffN_MBTETQ4y(x, y, z);
902  diffN[14] = diffN_MBTETQ4z(x, y, z);
903  diffN[15] = diffN_MBTETQ5x(x, y, z);
904  diffN[16] = diffN_MBTETQ5y(x, y, z);
905  diffN[17] = diffN_MBTETQ5z(x, y, z);
906  diffN[18] = diffN_MBTETQ6x(x, y, z);
907  diffN[19] = diffN_MBTETQ6y(x, y, z);
908  diffN[20] = diffN_MBTETQ6z(x, y, z);
909  diffN[21] = diffN_MBTETQ7x(x, y, z);
910  diffN[22] = diffN_MBTETQ7y(x, y, z);
911  diffN[23] = diffN_MBTETQ7z(x, y, z);
912  diffN[24] = diffN_MBTETQ8x(x, y, z);
913  diffN[25] = diffN_MBTETQ8y(x, y, z);
914  diffN[26] = diffN_MBTETQ8z(x, y, z);
915  diffN[27] = diffN_MBTETQ9x(x, y, z);
916  diffN[28] = diffN_MBTETQ9y(x, y, z);
917  diffN[29] = diffN_MBTETQ9z(x, y, z);
918  MoFEMFunctionReturnHot(0);
919 }
920 PetscErrorCode ShapeMBTETQ_GAUSS(double *N, const double *X, const double *Y,
921  const double *Z, const int G_DIM) {
922  MoFEMFunctionBeginHot;
923  int ii = 0;
924  for (; ii < G_DIM; ii++) {
925  double x = X[ii], y = Y[ii], z = Z[ii];
926  N[10 * ii + 0] = N_MBTETQ0(x, y, z);
927  N[10 * ii + 1] = N_MBTETQ1(x, y, z);
928  N[10 * ii + 2] = N_MBTETQ2(x, y, z);
929  N[10 * ii + 3] = N_MBTETQ3(x, y, z);
930  N[10 * ii + 4] = N_MBTETQ4(x, y, z);
931  N[10 * ii + 5] = N_MBTETQ5(x, y, z);
932  N[10 * ii + 6] = N_MBTETQ6(x, y, z);
933  N[10 * ii + 7] = N_MBTETQ7(x, y, z);
934  N[10 * ii + 8] = N_MBTETQ8(x, y, z);
935  N[10 * ii + 9] = N_MBTETQ9(x, y, z);
936  }
937  MoFEMFunctionReturnHot(0);
938 }
939 PetscErrorCode ShapeDiffMBTETQ_GAUSS(double *diffN, const double *X,
940  const double *Y, const double *Z,
941  const int G_DIM) {
942  MoFEMFunctionBeginHot;
943  int ii = 0;
944  for (; ii < G_DIM; ii++) {
945  double x = X[ii], y = Y[ii], z = Z[ii];
946  diffN[30 * ii + 0] = diffN_MBTETQ0x(x, y, z);
947  diffN[30 * ii + 1] = diffN_MBTETQ0y(x, y, z);
948  diffN[30 * ii + 2] = diffN_MBTETQ0z(x, y, z);
949  diffN[30 * ii + 3] = diffN_MBTETQ1x(x, y, z);
950  diffN[30 * ii + 4] = diffN_MBTETQ1y(x, y, z);
951  diffN[30 * ii + 5] = diffN_MBTETQ1z(x, y, z);
952  diffN[30 * ii + 6] = diffN_MBTETQ2x(x, y, z);
953  diffN[30 * ii + 7] = diffN_MBTETQ2y(x, y, z);
954  diffN[30 * ii + 8] = diffN_MBTETQ2z(x, y, z);
955  diffN[30 * ii + 9] = diffN_MBTETQ3x(x, y, z);
956  diffN[30 * ii + 10] = diffN_MBTETQ3y(x, y, z);
957  diffN[30 * ii + 11] = diffN_MBTETQ3z(x, y, z);
958  diffN[30 * ii + 12] = diffN_MBTETQ4x(x, y, z);
959  diffN[30 * ii + 13] = diffN_MBTETQ4y(x, y, z);
960  diffN[30 * ii + 14] = diffN_MBTETQ4z(x, y, z);
961  diffN[30 * ii + 15] = diffN_MBTETQ5x(x, y, z);
962  diffN[30 * ii + 16] = diffN_MBTETQ5y(x, y, z);
963  diffN[30 * ii + 17] = diffN_MBTETQ5z(x, y, z);
964  diffN[30 * ii + 18] = diffN_MBTETQ6x(x, y, z);
965  diffN[30 * ii + 19] = diffN_MBTETQ6y(x, y, z);
966  diffN[30 * ii + 20] = diffN_MBTETQ6z(x, y, z);
967  diffN[30 * ii + 21] = diffN_MBTETQ7x(x, y, z);
968  diffN[30 * ii + 22] = diffN_MBTETQ7y(x, y, z);
969  diffN[30 * ii + 23] = diffN_MBTETQ7z(x, y, z);
970  diffN[30 * ii + 24] = diffN_MBTETQ8x(x, y, z);
971  diffN[30 * ii + 25] = diffN_MBTETQ8y(x, y, z);
972  diffN[30 * ii + 26] = diffN_MBTETQ8z(x, y, z);
973  diffN[30 * ii + 27] = diffN_MBTETQ9x(x, y, z);
974  diffN[30 * ii + 28] = diffN_MBTETQ9y(x, y, z);
975  diffN[30 * ii + 29] = diffN_MBTETQ9z(x, y, z);
976  }
977  MoFEMFunctionReturnHot(0);
978 }
979 PetscErrorCode ShapeJacMBTETQ(const double *diffN, const double *coords,
980  double *Jac) {
981  MoFEMFunctionBeginHot;
982  int ii, jj, kk;
983  bzero(Jac, sizeof(double) * 9);
984  for (ii = 0; ii < 10; ii++) // shape func.
985  for (jj = 0; jj < 3; jj++) // space
986  for (kk = 0; kk < 3; kk++) // direvative of shape func.
987  Jac[jj * 3 + kk] += diffN[ii * 3 + kk] * coords[ii * 3 + jj];
988  MoFEMFunctionReturnHot(0);
989 }
990 PetscErrorCode
991 ShapeMBTETQ_detJac_at_Gauss_Points(double *detJac_at_Gauss_Points,
992  const double *diffN, const double *coords,
993  int G_DIM) {
994  MoFEMFunctionBeginHot;
995 
996  double Jac[9];
997  int ii = 0;
998  for (; ii < G_DIM; ii++) {
999  ierr = ShapeJacMBTETQ(&diffN[30 * ii], coords, Jac);
1000  CHKERRQ(ierr);
1001  detJac_at_Gauss_Points[ii] = ShapeDetJacVolume(Jac);
1002  }
1003  MoFEMFunctionReturnHot(0);
1004 }
1005 double ShapeVolumeMBTETQ(const double *diffN, const double *coords, int G_DIM,
1006  double *G_TET_W) {
1007 
1008  int ii = 0;
1009  double vol = 0;
1010  double detJac_at_Gauss_Points[G_DIM];
1011  ierr = ShapeMBTETQ_detJac_at_Gauss_Points(detJac_at_Gauss_Points, diffN,
1012  coords, G_DIM);
1013  CHKERRQ(ierr);
1014  for (; ii < G_DIM; ii++) {
1015  vol += G_TET_W[ii] * (detJac_at_Gauss_Points[ii]) / 6;
1016  }
1017  return vol;
1018 }
1019 PetscErrorCode ShapeMBTETQ_inverse(double *N, double *diffN,
1020  const double *elem_coords,
1021  const double *glob_coords,
1022  double *loc_coords, const double eps) {
1023  MoFEMFunctionBeginHot;
1024  double A[3 * 3];
1025  double R[3];
1026  int IPIV[3];
1027  float NORM_dR = 1000.;
1028  float NORM_R0;
1029  ShapeMBTETQ(N, 0.1, 0.1, 0.1);
1030  ShapeDiffMBTETQ(diffN, 0.1, 0.1, 0.1);
1031  R[0] = glob_coords[0] - cblas_ddot(10, &N[0], 1, &elem_coords[0], 3);
1032  R[1] = glob_coords[1] - cblas_ddot(10, &N[0], 1, &elem_coords[1], 3);
1033  R[2] = glob_coords[2] - cblas_ddot(10, &N[0], 1, &elem_coords[2], 3);
1034  NORM_R0 = cblas_dnrm2(3, &R[0], 1);
1035  while ((NORM_dR / NORM_R0) > eps) {
1036  // COL MAJOR
1037  // X
1038  A[0 + 3 * 0] = cblas_ddot(10, &diffN[0 * 3 + 0], 3, &elem_coords[0], 3);
1039  A[0 + 3 * 1] = cblas_ddot(10, &diffN[0 * 3 + 1], 3, &elem_coords[0], 3);
1040  A[0 + 3 * 2] = cblas_ddot(10, &diffN[0 * 3 + 2], 3, &elem_coords[0], 3);
1041  R[0] = glob_coords[0] - cblas_ddot(10, &N[0], 1, &elem_coords[0], 3);
1042  // Y
1043  A[1 + 3 * 0] = cblas_ddot(10, &diffN[0 * 3 + 0], 3, &elem_coords[1], 3);
1044  A[1 + 3 * 1] = cblas_ddot(10, &diffN[0 * 3 + 1], 3, &elem_coords[1], 3);
1045  A[1 + 3 * 2] = cblas_ddot(10, &diffN[0 * 3 + 2], 3, &elem_coords[1], 3);
1046  R[1] = glob_coords[1] - cblas_ddot(10, &N[0], 1, &elem_coords[1], 3);
1047  // Z
1048  A[2 + 3 * 0] =
1049  cblas_ddot(10, &diffN[0 * 3 + 0], 3, &elem_coords[0 * 3 + 2], 3);
1050  A[2 + 3 * 1] =
1051  cblas_ddot(10, &diffN[0 * 3 + 1], 3, &elem_coords[0 * 3 + 2], 3);
1052  A[2 + 3 * 2] =
1053  cblas_ddot(10, &diffN[0 * 3 + 2], 3, &elem_coords[0 * 3 + 2], 3);
1054  R[2] = glob_coords[2] - cblas_ddot(10, &N[0], 1, &elem_coords[2], 3);
1055  int info = lapack_dgesv(3, 1, &A[0], 3, (__CLPK_integer *)IPIV, R, 3);
1056  assert(info == 0);
1057  NOT_USED(info);
1058  cblas_daxpy(3, 1., R, 1, loc_coords, 1);
1059  NORM_dR = cblas_dnrm2(3, &R[0], 1);
1060  ShapeMBTETQ(N, loc_coords[0], loc_coords[1], loc_coords[2]);
1061  ShapeDiffMBTETQ(diffN, loc_coords[0], loc_coords[1], loc_coords[2]);
1062  }
1063  MoFEMFunctionReturnHot(0);
1064 }
I
static Index< 'I', 3 > I
Definition: BasicFeTools.hpp:88
invJac
MatrixDouble invJac
Definition: HookeElement.hpp:698
CblasRowMajor
@ CblasRowMajor
Definition: cblas.h:10
cblas_ddot
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: cblas_ddot.c:12
cblas_dnrm2
double cblas_dnrm2(const int N, const double *X, const int incX)
Definition: cblas_dnrm2.c:12
CBLAS_TRANSPOSE
CBLAS_TRANSPOSE
Definition: cblas.h:11
CblasNoTrans
@ CblasNoTrans
Definition: cblas.h:11
CblasTrans
@ CblasTrans
Definition: cblas.h:11
cblas_dgemv
void cblas_dgemv(const enum CBLAS_ORDER order, const enum CBLAS_TRANSPOSE TransA, const int M, const int N, const double alpha, const double *A, const int lda, const double *X, const int incX, const double beta, double *Y, const int incY)
Definition: cblas_dgemv.c:11
cblas_dcopy
void cblas_dcopy(const int N, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_dcopy.c:11
cblas_daxpy
void cblas_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_daxpy.c:11
cblas_dgemm
void cblas_dgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K, const double alpha, const double *A, const int lda, const double *B, const int ldb, const double beta, double *C, const int ldc)
Definition: cblas_dgemm.c:12
cblas_zgemv
void cblas_zgemv(const enum CBLAS_ORDER order, const enum CBLAS_TRANSPOSE TransA, const int M, const int N, const void *alpha, const void *A, const int lda, const void *X, const int incX, const void *beta, void *Y, const int incY)
Definition: cblas_zgemv.c:12
definitions.h
useful compiler directives and definitions
NOT_USED
#define NOT_USED(x)
Definition: definitions.h:311
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
ShapeVolumeMBTETQ
double ShapeVolumeMBTETQ(const double *diffN, const double *coords, int G_DIM, double *G_TET_W)
Definition: fem_tools.c:1005
MakeComplexTensor
PetscErrorCode MakeComplexTensor(double *reA, double *imA, __CLPK_doublecomplex *xA)
Definition: fem_tools.c:499
Normal_hierarchical
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 *xs1, __CLPK_doublecomplex *xs2, int gg)
Complex normal.
Definition: fem_tools.c:581
diffN_MBTETQ5x
#define diffN_MBTETQ5x(x, y, z)
Definition: fem_tools.c:855
diffN_MBTETQ1y
#define diffN_MBTETQ1y(x, y, z)
Definition: fem_tools.c:844
InvertComplexSymmMatrix3by3
PetscErrorCode InvertComplexSymmMatrix3by3(__CLPK_doublecomplex *xC)
Definition: fem_tools.c:525
diffN_MBTETQ1z
#define diffN_MBTETQ1z(x, y, z)
Definition: fem_tools.c:845
ShapeDiffMBTET
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:331
N_MBTETQ2
#define N_MBTETQ2(x, y, z)
Definition: fem_tools.c:832
ShapeJacMBTET
PetscErrorCode ShapeJacMBTET(double *diffN, const double *coords, double *jac)
calculate jacobian
Definition: fem_tools.c:300
N_MBTETQ9
#define N_MBTETQ9(x, y, z)
Definition: fem_tools.c:839
N_MBTETQ4
#define N_MBTETQ4(x, y, z)
Definition: fem_tools.c:834
GradientOfDeformation
PetscErrorCode GradientOfDeformation(double *diffN, double *dofs, double *F)
calculate gradient of deformation
Definition: fem_tools.c:437
N_MBTETQ8
#define N_MBTETQ8(x, y, z)
Definition: fem_tools.c:838
Grundmann_Moeller_integration_points_2D_TRI
PetscErrorCode Grundmann_Moeller_integration_points_2D_TRI(int rule, double *G_TRI_X, double *G_TRI_Y, double *G_TRI_W)
Definition: fem_tools.c:110
diffN_MBTETQ8z
#define diffN_MBTETQ8z(x, y, z)
Definition: fem_tools.c:866
diffN_MBTETQ7y
#define diffN_MBTETQ7y(x, y, z)
Definition: fem_tools.c:862
ShapeMBTET_inverse
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
diffN_MBTETQ3y
#define diffN_MBTETQ3y(x, y, z)
Definition: fem_tools.c:850
ShapeDiffMBTRIQ
PetscErrorCode ShapeDiffMBTRIQ(double *diffN, const double *X, const double *Y, const int G_DIM)
Definition: fem_tools.c:807
ShapeFaceDiffNormalMBTRI
PetscErrorCode ShapeFaceDiffNormalMBTRI(double *diffN, const double *coords, double *diff_normal)
calculate derivative of normal in respect to nodal positions
Definition: fem_tools.c:249
ShapeFaceBaseMBTRI
PetscErrorCode ShapeFaceBaseMBTRI(double *diffN, const double *coords, double *normal, double *s1, double *s2)
Definition: fem_tools.c:216
ShapeFaceNormalMBTRI_complex
PetscErrorCode ShapeFaceNormalMBTRI_complex(double *diffN, __CLPK_doublecomplex *xcoords, __CLPK_doublecomplex *xnormal)
Definition: fem_tools.c:464
N_MBTETQ6
#define N_MBTETQ6(x, y, z)
Definition: fem_tools.c:836
ShapeFaceNormalMBTRI
PetscErrorCode ShapeFaceNormalMBTRI(double *diffN, const double *coords, double *normal)
Definition: fem_tools.c:241
ShapeDiffMBTETQ
PetscErrorCode ShapeDiffMBTETQ(double *diffN, const double x, const double y, const double z)
Definition: fem_tools.c:885
diffN_MBTETQ0y
#define diffN_MBTETQ0y(x, y, z)
Definition: fem_tools.c:841
ierr
static PetscErrorCode ierr
Definition: fem_tools.c:32
diffN_MBTETQ9y
#define diffN_MBTETQ9y(x, y, z)
Definition: fem_tools.c:868
Grundmann_Moeller_integration_points_3D_TET
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)
Definition: fem_tools.c:155
diffN_MBTETQ9x
#define diffN_MBTETQ9x(x, y, z)
Definition: fem_tools.c:867
ShapeInvJacVolume
PetscErrorCode ShapeInvJacVolume(double *jac)
Definition: fem_tools.c:51
ShapeDiffMBTETinvJ
PetscErrorCode ShapeDiffMBTETinvJ(double *diffN, double *invJac, double *diffNinvJac)
calculate shape functions derivatives in space
Definition: fem_tools.c:427
diffN_MBTETQ0x
#define diffN_MBTETQ0x(x, y, z)
Definition: fem_tools.c:840
ShapeDiffMBTRI
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:206
diffN_MBTETQ6y
#define diffN_MBTETQ6y(x, y, z)
Definition: fem_tools.c:859
ShapeDiffMBTETinvJ_complex
void ShapeDiffMBTETinvJ_complex(double *diffN, __CLPK_doublecomplex *invJac, __CLPK_doublecomplex *diffNinvJac, const enum CBLAS_TRANSPOSE Trans)
Definition: fem_tools.c:449
ShapeMBTET
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
diffN_MBTETQ6x
#define diffN_MBTETQ6x(x, y, z)
Definition: fem_tools.c:858
ShapeDiffMBEDGE
PetscErrorCode ShapeDiffMBEDGE(double *diffN)
Definition: fem_tools.c:783
ShapeMBTETQ_GAUSS
PetscErrorCode ShapeMBTETQ_GAUSS(double *N, const double *X, const double *Y, const double *Z, const int G_DIM)
Definition: fem_tools.c:920
DeterminantComplexGradient
PetscErrorCode DeterminantComplexGradient(__CLPK_doublecomplex *xF, __CLPK_doublecomplex *det_xF)
Definition: fem_tools.c:538
Spin
PetscErrorCode Spin(double *spinOmega, double *vecOmega)
calculate spin matrix from vector
Definition: fem_tools.c:558
diffN_MBTETQ2z
#define diffN_MBTETQ2z(x, y, z)
Definition: fem_tools.c:848
diffN_MBTETQ6z
#define diffN_MBTETQ6z(x, y, z)
Definition: fem_tools.c:860
diffN_MBTETQ4z
#define diffN_MBTETQ4z(x, y, z)
Definition: fem_tools.c:854
diffN_MBTETQ0z
#define diffN_MBTETQ0z(x, y, z)
Definition: fem_tools.c:842
Grundmann_Moeller_integration_points_1D_EDGE
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
N_MBTETQ7
#define N_MBTETQ7(x, y, z)
Definition: fem_tools.c:837
InvertComplexGradient
PetscErrorCode InvertComplexGradient(__CLPK_doublecomplex *xF)
Definition: fem_tools.c:511
diffN_MBTETQ3z
#define diffN_MBTETQ3z(x, y, z)
Definition: fem_tools.c:851
diffN_MBTETQ4x
#define diffN_MBTETQ4x(x, y, z)
Definition: fem_tools.c:852
diffN_MBTETQ2y
#define diffN_MBTETQ2y(x, y, z)
Definition: fem_tools.c:847
N_MBTETQ5
#define N_MBTETQ5(x, y, z)
Definition: fem_tools.c:835
diffN_MBTETQ9z
#define diffN_MBTETQ9z(x, y, z)
Definition: fem_tools.c:869
diffN_MBTETQ5y
#define diffN_MBTETQ5y(x, y, z)
Definition: fem_tools.c:856
diffN_MBTETQ4y
#define diffN_MBTETQ4y(x, y, z)
Definition: fem_tools.c:853
ShapeMBTRIQ
PetscErrorCode ShapeMBTRIQ(double *N, const double *X, const double *Y, const int G_DIM)
Definition: fem_tools.c:792
ShapeMBTETQ
PetscErrorCode ShapeMBTETQ(double *N, const double x, const double y, const double z)
Definition: fem_tools.c:870
diffN_MBTETQ3x
#define diffN_MBTETQ3x(x, y, z)
Definition: fem_tools.c:849
ShapeMBEDGE
PetscErrorCode ShapeMBEDGE(double *N, const double *G_X, int DIM)
Definition: fem_tools.c:773
ShapeVolumeMBTET
double ShapeVolumeMBTET(double *diffN, const double *coords)
calculate TET volume
Definition: fem_tools.c:310
diffN_MBTETQ2x
#define diffN_MBTETQ2x(x, y, z)
Definition: fem_tools.c:846
N_MBTETQ1
#define N_MBTETQ1(x, y, z)
Definition: fem_tools.c:831
ShapeMBTETQ_inverse
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
make_complex_matrix
PetscErrorCode make_complex_matrix(double *reA, double *imA, __CLPK_doublecomplex *xA)
Compose complex matrix (3x3) from two real matrices.
Definition: fem_tools.c:569
ShapeDiffMBTETQ_GAUSS
PetscErrorCode ShapeDiffMBTETQ_GAUSS(double *diffN, const double *X, const double *Y, const double *Z, const int G_DIM)
Definition: fem_tools.c:939
N_MBTETQ3
#define N_MBTETQ3(x, y, z)
Definition: fem_tools.c:833
ShapeMBTETQ_detJac_at_Gauss_Points
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
ShapeDetJacVolume
double ShapeDetJacVolume(double *jac)
determined of jacobian
Definition: fem_tools.c:34
diffN_MBTETQ7x
#define diffN_MBTETQ7x(x, y, z)
Definition: fem_tools.c:861
Base_scale
PetscErrorCode Base_scale(__CLPK_doublecomplex *xnormal, __CLPK_doublecomplex *xs1, __CLPK_doublecomplex *xs2)
Definition: fem_tools.c:753
diffN_MBTETQ1x
#define diffN_MBTETQ1x(x, y, z)
Definition: fem_tools.c:843
diffN_MBTETQ8y
#define diffN_MBTETQ8y(x, y, z)
Definition: fem_tools.c:865
diffN_MBTETQ7z
#define diffN_MBTETQ7z(x, y, z)
Definition: fem_tools.c:863
ShapeJacMBTETQ
PetscErrorCode ShapeJacMBTETQ(const double *diffN, const double *coords, double *Jac)
Definition: fem_tools.c:979
ShapeMBTRI_inverse
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
N_MBTETQ0
#define N_MBTETQ0(x, y, z)
Definition: fem_tools.c:830
diffN_MBTETQ5z
#define diffN_MBTETQ5z(x, y, z)
Definition: fem_tools.c:857
diffN_MBTETQ8x
#define diffN_MBTETQ8x(x, y, z)
Definition: fem_tools.c:864
ShapeMBTRI
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition: fem_tools.c:194
fem_tools.h
Loose implementation of some useful functions.
diffN_MBTET0z
#define diffN_MBTET0z
derivative of tetrahedral shape function
Definition: fem_tools.h:44
N_MBTET0
#define N_MBTET0(x, y, z)
tetrahedral shape function
Definition: fem_tools.h:38
N_MBTRIQ2
#define N_MBTRIQ2(x, y)
Definition: fem_tools.h:89
diffN_MBTRIQ0x
#define diffN_MBTRIQ0x(x, y)
Definition: fem_tools.h:93
diffN_MBTRIQ1x
#define diffN_MBTRIQ1x(x, y)
Definition: fem_tools.h:95
diffN_MBTRI2x
#define diffN_MBTRI2x
derivative of triangle shape function
Definition: fem_tools.h:63
N_MBTRIQ3
#define N_MBTRIQ3(x, y)
Definition: fem_tools.h:90
N_MBTET1
#define N_MBTET1(x, y, z)
tetrahedral shape function
Definition: fem_tools.h:39
diffN_MBTRI0y
#define diffN_MBTRI0y
derivative of triangle shape function
Definition: fem_tools.h:60
diffN_MBTRI1y
#define diffN_MBTRI1y
derivative of triangle shape function
Definition: fem_tools.h:62
N_MBTRI1
#define N_MBTRI1(x, y)
triangle shape function
Definition: fem_tools.h:57
diffN_MBTRIQ5x
#define diffN_MBTRIQ5x(x, y)
Definition: fem_tools.h:103
diffN_MBTRIQ2x
#define diffN_MBTRIQ2x(x, y)
Definition: fem_tools.h:97
diffN_MBEDGE0
#define diffN_MBEDGE0
derivative of edge shape function
Definition: fem_tools.h:83
diffN_MBTRIQ2y
#define diffN_MBTRIQ2y(x, y)
Definition: fem_tools.h:98
diffN_MBTET2z
#define diffN_MBTET2z
derivative of tetrahedral shape function
Definition: fem_tools.h:50
diffN_MBEDGE1
#define diffN_MBEDGE1
derivative of edge shape function
Definition: fem_tools.h:84
N_MBEDGE0
#define N_MBEDGE0(x)
edge shape function
Definition: fem_tools.h:81
diffN_MBTET3y
#define diffN_MBTET3y
derivative of tetrahedral shape function
Definition: fem_tools.h:52
diffN_MBTET0y
#define diffN_MBTET0y
derivative of tetrahedral shape function
Definition: fem_tools.h:43
G_TET_W1
static const double G_TET_W1[]
Definition: fem_tools.h:1091
diffN_MBTRIQ1y
#define diffN_MBTRIQ1y(x, y)
Definition: fem_tools.h:96
diffN_MBTET1z
#define diffN_MBTET1z
derivative of tetrahedral shape function
Definition: fem_tools.h:47
diffN_MBTRIQ4y
#define diffN_MBTRIQ4y(x, y)
Definition: fem_tools.h:102
diffN_MBTET3z
#define diffN_MBTET3z
derivative of tetrahedral shape function
Definition: fem_tools.h:53
diffN_MBTRIQ5y
#define diffN_MBTRIQ5y(x, y)
Definition: fem_tools.h:104
diffN_MBTRIQ3x
#define diffN_MBTRIQ3x(x, y)
Definition: fem_tools.h:99
N_MBTRIQ0
#define N_MBTRIQ0(x, y)
Definition: fem_tools.h:87
N_MBTRIQ4
#define N_MBTRIQ4(x, y)
Definition: fem_tools.h:91
diffN_MBTET1y
#define diffN_MBTET1y
derivative of tetrahedral shape function
Definition: fem_tools.h:46
diffN_MBTRI0x
#define diffN_MBTRI0x
derivative of triangle shape function
Definition: fem_tools.h:59
N_MBTRIQ5
#define N_MBTRIQ5(x, y)
Definition: fem_tools.h:92
diffN_MBTET2y
#define diffN_MBTET2y
derivative of tetrahedral shape function
Definition: fem_tools.h:49
diffN_MBTET2x
#define diffN_MBTET2x
derivative of tetrahedral shape function
Definition: fem_tools.h:48
N_MBTRI0
#define N_MBTRI0(x, y)
triangle shape function
Definition: fem_tools.h:56
diffN_MBTRI2y
#define diffN_MBTRI2y
derivative of triangle shape function
Definition: fem_tools.h:64
N_MBTRIQ1
#define N_MBTRIQ1(x, y)
Definition: fem_tools.h:88
diffN_MBTRIQ3y
#define diffN_MBTRIQ3y(x, y)
Definition: fem_tools.h:100
N_MBTRI2
#define N_MBTRI2(x, y)
triangle shape function
Definition: fem_tools.h:58
N_MBTET2
#define N_MBTET2(x, y, z)
tetrahedral shape function
Definition: fem_tools.h:40
diffN_MBTRIQ4x
#define diffN_MBTRIQ4x(x, y)
Definition: fem_tools.h:101
N_MBEDGE1
#define N_MBEDGE1(x)
edge shape function
Definition: fem_tools.h:82
diffN_MBTRI1x
#define diffN_MBTRI1x
derivative of triangle shape function
Definition: fem_tools.h:61
diffN_MBTET3x
#define diffN_MBTET3x
derivative of tetrahedral shape function
Definition: fem_tools.h:51
diffN_MBTRIQ0y
#define diffN_MBTRIQ0y(x, y)
Definition: fem_tools.h:94
diffN_MBTET1x
#define diffN_MBTET1x
derivative of tetrahedral shape function
Definition: fem_tools.h:45
diffN_MBTET0x
#define diffN_MBTET0x
derivative of tetrahedral shape function
Definition: fem_tools.h:42
N_MBTET3
#define N_MBTET3(x, y, z)
tetrahedral shape function
Definition: fem_tools.h:41
eps
constexpr double eps
Definition: forces_and_sources_getting_mult_H1_H1_atom_test.cpp:30
gm_rule_set
void gm_rule_set(int rule, int dim_num, int point_num, double w[], double x[])
Definition: gm_rule.c:152
gm_rule_size
int gm_rule_size(int rule, int dim_num)
Definition: gm_rule.c:294
gm_rule.h
h1_hdiv_hcurl_l2.h
Functions to approximate hierarchical spaces.
NBEDGE_H1
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
Definition: h1_hdiv_hcurl_l2.h:55
NBFACETRI_H1
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
Definition: h1_hdiv_hcurl_l2.h:60
lapack_zgetri
static __CLPK_integer lapack_zgetri(__CLPK_integer n, __CLPK_doublecomplex *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublecomplex *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:345
lapack_zpotri
static __CLPK_integer lapack_zpotri(char uplo, __CLPK_integer n, __CLPK_doublecomplex *a, __CLPK_integer lda)
Definition: lapack_wrap.h:353
lapack_zgetrf
static __CLPK_integer lapack_zgetrf(__CLPK_integer m, __CLPK_integer n, __CLPK_doublecomplex *a, __CLPK_integer lda, __CLPK_integer *ipiv)
Definition: lapack_wrap.h:335
lapack_dgesv
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition: lapack_wrap.h:188
__CLPK_integer
long int __CLPK_integer
Definition: lapack_wrap.h:35
lapack_dgetrf
static __CLPK_integer lapack_dgetrf(__CLPK_integer m, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv)
Definition: lapack_wrap.h:169
lapack_dgetri
static __CLPK_integer lapack_dgetri(__CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:197
lapack_zpotrf
static __CLPK_integer lapack_zpotrf(char uplo, __CLPK_integer n, __CLPK_doublecomplex *a, __CLPK_integer lda)
Definition: lapack_wrap.h:360
__CLPK_doublereal
double __CLPK_doublereal
Definition: lapack_wrap.h:38
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
i
FTensor::Index< 'i', 3 > i
Definition: matrix_function.cpp:18
ContactOps::normal
FTensor::Tensor1< double, SPACE_DIM > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_disp)
Definition: ContactOps.hpp:80
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
OpContactTools::w
double w(const double g, const double t)
Definition: ContactOperators.hpp:144
ReactionDiffusionEquation::r
const double r
rate factor
Definition: reaction_diffusion.cpp:33
N
const int N
Definition: speed_test.cpp:3
__CLPK_doublecomplex
Definition: lapack_wrap.h:46
__CLPK_doublecomplex::i
__CLPK_doublereal i
Definition: lapack_wrap.h:47
__CLPK_doublecomplex::r
__CLPK_doublereal r
Definition: lapack_wrap.h:47
Generated on Tue Feb 16 2021 20:34:35 for MoFEM by Doxygen 1.9.1 and hosted at University of Glasgow