v0.14.0
Functions | Variables
h1.c File Reference
#include <petscsys.h>
#include <definitions.h>
#include <cblas.h>
#include <h1_hdiv_hcurl_l2.h>

Go to the source code of this file.

Functions

PetscErrorCode H1_EdgeShapeFunctions_MBTRI (int *sense, int *p, double *N, double *diffN, double *edgeN[3], double *diff_edgeN[3], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 H1_EdgeShapeFunctions_MBTRI. More...
 
PetscErrorCode H1_FaceShapeFunctions_MBTRI (const int *face_nodes, int p, double *N, double *diffN, double *faceN, double *diff_faceN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 
PetscErrorCode H1_EdgeShapeFunctions_MBTET (int *sense, int *p, double *N, double *diffN, double *edgeN[], double *diff_edgeN[], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 
PetscErrorCode H1_FaceShapeFunctions_MBTET (int *faces_nodes, int *p, double *N, double *diffN, double *faceN[], double *diff_faceN[], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 
PetscErrorCode H1_VolumeShapeFunctions_MBTET (int p, double *N, double *diffN, double *volumeN, double *diff_volumeN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 
PetscErrorCode H1_EdgeShapeDiffMBTETinvJ (int *base_p, int *p, double *edge_diffN[], double *invJac, double *edge_diffNinvJac[], int GDIM)
 
PetscErrorCode H1_FaceShapeDiffMBTETinvJ (int *base_p, int *p, double *face_diffN[], double *invJac, double *face_diffNinvJac[], int GDIM)
 
PetscErrorCode H1_VolumeShapeDiffMBTETinvJ (int base_p, int p, double *volume_diffN, double *invJac, double *volume_diffNinvJac, int GDIM)
 
PetscErrorCode H1_EdgeGradientOfDeformation_hierarchical (int p, double *diffN, double *dofs, double *F)
 
PetscErrorCode H1_FaceGradientOfDeformation_hierarchical (int p, double *diffN, double *dofs, double *F)
 
PetscErrorCode H1_VolumeGradientOfDeformation_hierarchical (int p, double *diffN, double *dofs, double *F)
 
PetscErrorCode H1_QuadShapeFunctions_MBPRISM (int *faces_nodes, int *p, double *N, double *diffN, double *faceN[], double *diff_faceN[], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 
PetscErrorCode H1_VolumeShapeFunctions_MBPRISM (int p, double *N, double *diffN, double *volumeN, double *diff_volumeN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 
PetscErrorCode H1_QuadShapeFunctions_MBQUAD (int *faces_nodes, int p, double *N, double *diffN, double *faceN, double *diff_faceN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 
PetscErrorCode H1_EdgeShapeFunctions_MBQUAD (int *sense, int *p, double *N, double *diffN, double *edgeN[4], double *diff_edgeN[4], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
 

Variables

static PetscErrorCode ierr
 

Detailed Description

Based on Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes, by Mark Ainsworth and Joe Coyle Shape functions for MBTRI and H1 approximation

Definition in file h1.c.

Function Documentation

◆ H1_EdgeGradientOfDeformation_hierarchical()

PetscErrorCode H1_EdgeGradientOfDeformation_hierarchical ( int  p,
double *  diffN,
double *  dofs,
double *  F 
)

Definition at line 609 of file h1.c.

611  {
612  MoFEMFunctionBeginHot;
613  int col, row = 0;
614  for (; row < 3; row++)
615  for (col = 0; col < 3; col++)
616  F[3 * row + col] =
617  cblas_ddot(NBEDGE_H1(p), &diffN[col], 3, &dofs[row], 3);
618  MoFEMFunctionReturnHot(0);
619 }

◆ H1_EdgeShapeDiffMBTETinvJ()

PetscErrorCode H1_EdgeShapeDiffMBTETinvJ ( int *  base_p,
int *  p,
double *  edge_diffN[],
double *  invJac,
double *  edge_diffNinvJac[],
int  GDIM 
)

Definition at line 556 of file h1.c.

558  {
559  MoFEMFunctionBeginHot;
560  int ee = 0, ii, gg;
561  for (; ee < 6; ee++) {
562  for (ii = 0; ii < NBEDGE_H1(p[ee]); ii++) {
563  for (gg = 0; gg < GDIM; gg++) {
564  int shift1 = NBEDGE_H1(base_p[ee]) * gg;
565  int shift2 = NBEDGE_H1(p[ee]) * gg;
566  cblas_dgemv(CblasRowMajor, CblasTrans, 3, 3, 1., invJac, 3,
567  &(edge_diffN[ee])[3 * shift1 + 3 * ii], 1, 0.,
568  &(edge_diffNinvJac[ee])[3 * shift2 + 3 * ii], 1);
569  }
570  }
571  }
572  MoFEMFunctionReturnHot(0);
573 }

◆ H1_EdgeShapeFunctions_MBQUAD()

PetscErrorCode H1_EdgeShapeFunctions_MBQUAD ( int *  sense,
int *  p,
double *  N,
double *  diffN,
double *  edgeN[4],
double *  diff_edgeN[4],
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)  base_polynomials 
)

Definition at line 1091 of file h1.c.

1096  {
1097  MoFEMFunctionBeginHot;
1098 
1099  double *edgeN01 = NULL, *edgeN12 = NULL, *edgeN23 = NULL, *edgeN30 = NULL;
1100  if (edgeN != NULL) {
1101  edgeN01 = edgeN[0];
1102  edgeN12 = edgeN[1];
1103  edgeN23 = edgeN[2];
1104  edgeN30 = edgeN[3];
1105  }
1106  double *diff_edgeN01 = NULL, *diff_edgeN12 = NULL, *diff_edgeN23 = NULL,
1107  *diff_edgeN30 = NULL;
1108  if (diff_edgeN != NULL) {
1109  diff_edgeN01 = diff_edgeN[0];
1110  diff_edgeN12 = diff_edgeN[1];
1111  diff_edgeN23 = diff_edgeN[2];
1112  diff_edgeN30 = diff_edgeN[3];
1113  }
1114  int P[4];
1115  int ee = 0;
1116  for (; ee < 4; ee++)
1117  P[ee] = NBEDGE_H1(p[ee]);
1118 
1119  int n0 = 0;
1120  int n1 = 1;
1121  int n2 = 2;
1122  int n3 = 3;
1123 
1124  int ii = 0;
1125  for (; ii < GDIM; ii++) {
1126  int node_shift = ii * 4;
1127  int node_diff_shift = 2 * node_shift;
1128 
1129  double shape0 = N[node_shift + n0];
1130  double shape1 = N[node_shift + n1];
1131  double shape2 = N[node_shift + n2];
1132  double shape3 = N[node_shift + n3];
1133 
1134  double ksi01 = (shape1 + shape2 - shape0 - shape3) * sense[n0];
1135  double ksi12 = (shape2 + shape3 - shape1 - shape0) * sense[n1];
1136  double ksi23 = (shape3 + shape0 - shape2 - shape1) * sense[n2];
1137  double ksi30 = (shape0 + shape1 - shape3 - shape2) * sense[n3];
1138 
1139  double extrude_zeta01 = shape0 + shape1;
1140  double extrude_ksi12 = shape1 + shape2;
1141  double extrude_zeta23 = shape2 + shape3;
1142  double extrude_ksi30 = shape0 + shape3;
1143 
1144  double bubble_ksi = extrude_ksi12 * extrude_ksi30;
1145  double bubble_zeta = extrude_zeta01 * extrude_zeta23;
1146 
1147  double diff_ksi01[2], diff_ksi12[2], diff_ksi23[2], diff_ksi30[2];
1148  double diff_extrude_zeta01[2];
1149  double diff_extrude_ksi12[2];
1150  double diff_extrude_zeta23[2];
1151  double diff_extrude_ksi30[2];
1152  double diff_bubble_ksi[2];
1153  double diff_bubble_zeta[2];
1154 
1155  int d = 0;
1156  for (; d < 2; d++) {
1157  double diff_shape0 = diffN[node_diff_shift + 2 * n0 + d];
1158  double diff_shape1 = diffN[node_diff_shift + 2 * n1 + d];
1159  double diff_shape2 = diffN[node_diff_shift + 2 * n2 + d];
1160  double diff_shape3 = diffN[node_diff_shift + 2 * n3 + d];
1161  diff_ksi01[d] =
1162  (diff_shape1 + diff_shape2 - diff_shape0 - diff_shape3) * sense[n0];
1163  diff_ksi12[d] =
1164  (diff_shape2 + diff_shape3 - diff_shape1 - diff_shape0) * sense[n1];
1165  diff_ksi23[d] =
1166  (diff_shape3 + diff_shape0 - diff_shape2 - diff_shape1) * sense[n2];
1167  diff_ksi30[d] =
1168  (diff_shape0 + diff_shape1 - diff_shape3 - diff_shape2) * sense[n3];
1169  diff_extrude_zeta01[d] = diff_shape0 + diff_shape1;
1170  diff_extrude_ksi12[d] = diff_shape1 + diff_shape2;
1171  diff_extrude_zeta23[d] = diff_shape2 + diff_shape3;
1172  diff_extrude_ksi30[d] = diff_shape0 + diff_shape3;
1173  diff_bubble_ksi[d] = diff_extrude_ksi12[d] * extrude_ksi30 +
1174  extrude_ksi12 * diff_extrude_ksi30[d];
1175  diff_bubble_zeta[d] = diff_extrude_zeta01[d] * extrude_zeta23 +
1176  extrude_zeta01 * diff_extrude_zeta23[d];
1177  }
1178 
1179  double L01[p[0] + 1], L12[p[1] + 1], L23[p[2] + 1], L30[p[3] + 1];
1180  double diffL01[2 * (p[0] + 1)], diffL12[2 * (p[1] + 1)],
1181  diffL23[2 * (p[2] + 1)], diffL30[2 * (p[3] + 1)];
1182  ierr = base_polynomials(p[0], ksi01, diff_ksi01, L01, diffL01, 2);
1183  CHKERRQ(ierr);
1184  ierr = base_polynomials(p[1], ksi12, diff_ksi12, L12, diffL12, 2);
1185  CHKERRQ(ierr);
1186  ierr = base_polynomials(p[2], ksi23, diff_ksi23, L23, diffL23, 2);
1187  CHKERRQ(ierr);
1188  ierr = base_polynomials(p[3], ksi30, diff_ksi30, L30, diffL30, 2);
1189  CHKERRQ(ierr);
1190 
1191  int shift;
1192  if (edgeN != NULL) {
1193  // edge01
1194  shift = ii * (P[0]);
1195  cblas_dcopy(P[0], L01, 1, &edgeN01[shift], 1);
1196  cblas_dscal(P[0], bubble_ksi * extrude_zeta01, &edgeN01[shift], 1);
1197  // edge12
1198  shift = ii * (P[1]);
1199  cblas_dcopy(P[1], L12, 1, &edgeN12[shift], 1);
1200  cblas_dscal(P[1], bubble_zeta * extrude_ksi12, &edgeN12[shift], 1);
1201  // edge23
1202  shift = ii * (P[2]);
1203  cblas_dcopy(P[2], L23, 1, &edgeN23[shift], 1);
1204  cblas_dscal(P[2], bubble_ksi * extrude_zeta23, &edgeN23[shift], 1);
1205  // edge30
1206  shift = ii * (P[3]);
1207  cblas_dcopy(P[3], L30, 1, &edgeN30[shift], 1);
1208  cblas_dscal(P[3], bubble_zeta * extrude_ksi30, &edgeN30[shift], 1);
1209  }
1210  if (diff_edgeN != NULL) {
1211  if (P[0] > 0) {
1212  // edge01
1213  shift = ii * (P[0]);
1214  bzero(&diff_edgeN01[2 * shift], sizeof(double) * 2 * (P[0]));
1215  int d = 0;
1216  for (; d != 2; ++d) {
1217  cblas_daxpy(P[0], bubble_ksi * extrude_zeta01,
1218  &diffL01[d * (p[0] + 1)], 1, &diff_edgeN01[2 * shift + d],
1219  2);
1220  cblas_daxpy(P[0],
1221  diff_bubble_ksi[d] * extrude_zeta01 +
1222  bubble_ksi * diff_extrude_zeta01[d],
1223  L01, 1, &diff_edgeN01[2 * shift + d], 2);
1224  }
1225  }
1226  if (P[1] > 0) {
1227  // edge12
1228  shift = ii * (P[1]);
1229  bzero(&diff_edgeN12[2 * shift], sizeof(double) * 2 * (P[1]));
1230  int d = 0;
1231  for (; d != 2; ++d) {
1232  cblas_daxpy(P[1], bubble_zeta * extrude_ksi12,
1233  &diffL12[d * (p[1] + 1)], 1, &diff_edgeN12[2 * shift + d],
1234  2);
1235  cblas_daxpy(P[1],
1236  diff_bubble_zeta[d] * extrude_ksi12 +
1237  bubble_zeta * diff_extrude_ksi12[d],
1238  L12, 1, &diff_edgeN12[2 * shift + d], 2);
1239  }
1240  }
1241  if (P[2] > 0) {
1242  // edge23
1243  shift = ii * (P[2]);
1244  bzero(&diff_edgeN23[2 * shift], sizeof(double) * 2 * (P[2]));
1245  int d = 0;
1246  for (; d != 2; ++d) {
1247  cblas_daxpy(P[2], bubble_ksi * extrude_zeta23,
1248  &diffL23[d * (p[2] + 1)], 1, &diff_edgeN23[2 * shift + d],
1249  2);
1250  cblas_daxpy(P[2],
1251  diff_bubble_ksi[d] * extrude_zeta23 +
1252  bubble_ksi * diff_extrude_zeta23[d],
1253  L23, 1, &diff_edgeN23[2 * shift + d], 2);
1254  }
1255  }
1256  if (P[3] > 0) {
1257  // edge30
1258  shift = ii * (P[3]);
1259  bzero(&diff_edgeN30[2 * shift], sizeof(double) * 2 * (P[3]));
1260  int d = 0;
1261  for (; d != 2; ++d) {
1262  cblas_daxpy(P[3], bubble_zeta * extrude_ksi30,
1263  &diffL30[d * (p[3] + 1)], 1, &diff_edgeN30[2 * shift + d],
1264  2);
1265  cblas_daxpy(P[3],
1266  diff_bubble_zeta[d] * extrude_ksi30 +
1267  bubble_zeta * diff_extrude_ksi30[d],
1268  L30, 1, &diff_edgeN30[2 * shift + d], 2);
1269  }
1270  }
1271  }
1272  }
1273  MoFEMFunctionReturnHot(0);
1274 }

◆ H1_EdgeShapeFunctions_MBTET()

PetscErrorCode H1_EdgeShapeFunctions_MBTET ( int *  sense,
int *  p,
double *  N,
double *  diffN,
double *  edgeN[],
double *  diff_edgeN[],
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)  base_polynomials 
)

Definition at line 274 of file h1.c.

279  {
280  MoFEMFunctionBeginHot;
281 
282  int P[6];
283  int ee = 0;
284  for (; ee < 6; ee++)
285  P[ee] = NBEDGE_H1(p[ee]);
286  int edges_nodes[2 * 6] = {0, 1, 1, 2, 2, 0, 0, 3, 1, 3, 2, 3};
287  double diff_ksi01[3], diff_ksi12[3], diff_ksi20[3];
288  double diff_ksi03[3], diff_ksi13[3], diff_ksi23[3];
289  double *edges_diff_ksi[6] = {diff_ksi01, diff_ksi12, diff_ksi20,
290  diff_ksi03, diff_ksi13, diff_ksi23};
291  for (ee = 0; ee < 6; ee++) {
292  int dd = 0;
293  for (; dd < 3; dd++) {
294  edges_diff_ksi[ee][dd] = (diffN[edges_nodes[2 * ee + 1] * 3 + dd] -
295  diffN[edges_nodes[2 * ee + 0] * 3 + dd]) *
296  sense[ee];
297  }
298  }
299  int ii = 0;
300  for (; ii < GDIM; ii++) {
301  int node_shift = ii * 4;
302  double edges_ksi[6];
303  for (ee = 0; ee < 6; ee++) {
304  edges_ksi[ee] = (N[node_shift + edges_nodes[2 * ee + 1]] -
305  N[node_shift + edges_nodes[2 * ee + 0]]) *
306  sense[ee];
307  }
308  for (ee = 0; ee < 6; ee++) {
309  if (P[ee] == 0)
310  continue;
311  double L[p[ee] + 1], diffL[3 * (p[ee] + 1)];
312  ierr = base_polynomials(p[ee], edges_ksi[ee], edges_diff_ksi[ee], L,
313  diffL, 3);
314  CHKERRQ(ierr);
315  double v = N[node_shift + edges_nodes[2 * ee + 0]] *
316  N[node_shift + edges_nodes[2 * ee + 1]];
317  if (edgeN != NULL)
318  if (edgeN[ee] != NULL) {
319  int shift = ii * P[ee];
320  {
321  double *edge_n_ptr = &edgeN[ee][shift];
322  double *l_ptr = L;
323  double scalar = N[node_shift + edges_nodes[2 * ee + 0]] *
324  N[node_shift + edges_nodes[2 * ee + 1]];
325  int size = P[ee];
326  for (size_t jj = 0; jj != size; ++jj, ++l_ptr) {
327  *edge_n_ptr = (*l_ptr) * scalar;
328  ++edge_n_ptr;
329  }
330  }
331  }
332  if (diff_edgeN != NULL)
333  if (diff_edgeN[ee] != NULL) {
334  int shift = ii * P[ee];
335  {
336  double *diff_edge_n_ptr = &diff_edgeN[ee][3 * shift];
337  double *diff_l_x = &diffL[0 * (p[ee] + 1)];
338  double *diff_l_y = &diffL[1 * (p[ee] + 1)];
339  double *diff_l_z = &diffL[2 * (p[ee] + 1)];
340  double *l_ptr = L;
341  double scalar_x = diffN[3 * edges_nodes[2 * ee + 0] + 0] *
342  N[node_shift + edges_nodes[2 * ee + 1]] +
343  N[node_shift + edges_nodes[2 * ee + 0]] *
344  diffN[3 * edges_nodes[2 * ee + 1] + 0];
345  double scalar_y = diffN[3 * edges_nodes[2 * ee + 0] + 1] *
346  N[node_shift + edges_nodes[2 * ee + 1]] +
347  N[node_shift + edges_nodes[2 * ee + 0]] *
348  diffN[3 * edges_nodes[2 * ee + 1] + 1];
349  double scalar_z = diffN[3 * edges_nodes[2 * ee + 0] + 2] *
350  N[node_shift + edges_nodes[2 * ee + 1]] +
351  N[node_shift + edges_nodes[2 * ee + 0]] *
352  diffN[3 * edges_nodes[2 * ee + 1] + 2];
353 
354  int size = P[ee];
355  for (size_t jj = 0; jj != size;
356  ++jj, ++diff_l_x, ++diff_l_y, ++diff_l_z, ++l_ptr) {
357 
358  *diff_edge_n_ptr = v * (*diff_l_x) + scalar_x * (*l_ptr);
359  ++diff_edge_n_ptr;
360  *diff_edge_n_ptr = v * (*diff_l_y) + scalar_y * (*l_ptr);
361  ++diff_edge_n_ptr;
362  *diff_edge_n_ptr = v * (*diff_l_z) + scalar_z * (*l_ptr);
363  ++diff_edge_n_ptr;
364 
365  }
366  }
367 
368  }
369  }
370  }
371  MoFEMFunctionReturnHot(0);
372 }

◆ H1_EdgeShapeFunctions_MBTRI()

PetscErrorCode H1_EdgeShapeFunctions_MBTRI ( int *  sense,
int *  p,
double *  N,
double *  diffN,
double *  edgeN[3],
double *  diff_edgeN[3],
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)  base_polynomials 
)

H1_EdgeShapeFunctions_MBTRI.

Parameters
senseof edges, it is array of integers dim 3 (3-edges of triangle)
pof edges

Definition at line 17 of file h1.c.

22  {
23  MoFEMFunctionBeginHot;
24 
25  double *edgeN01 = NULL, *edgeN12 = NULL, *edgeN20 = NULL;
26  if (edgeN != NULL) {
27  edgeN01 = edgeN[0];
28  edgeN12 = edgeN[1];
29  edgeN20 = edgeN[2];
30  }
31  double *diff_edgeN01 = NULL, *diff_edgeN12 = NULL, *diff_edgeN20 = NULL;
32  if (diff_edgeN != NULL) {
33  diff_edgeN01 = diff_edgeN[0];
34  diff_edgeN12 = diff_edgeN[1];
35  diff_edgeN20 = diff_edgeN[2];
36  }
37  int P[3];
38  int ee = 0;
39  for (; ee < 3; ee++) {
40  P[ee] = NBEDGE_H1(p[ee]);
41  }
42  int dd = 0;
43  double diff_ksi01[2], diff_ksi12[2], diff_ksi20[2];
44  if (diff_edgeN != NULL) {
45  for (; dd < 2; dd++) {
46  diff_ksi01[dd] = (diffN[1 * 2 + dd] - diffN[0 * 2 + dd]) * sense[0];
47  diff_ksi12[dd] = (diffN[2 * 2 + dd] - diffN[1 * 2 + dd]) * sense[1];
48  diff_ksi20[dd] = (diffN[0 * 2 + dd] - diffN[2 * 2 + dd]) * sense[2];
49  }
50  }
51  int ii = 0;
52  for (; ii < GDIM; ii++) {
53  int node_shift = ii * 3;
54  double ksi01 = (N[node_shift + 1] - N[node_shift + 0]) * sense[0];
55  double ksi12 = (N[node_shift + 2] - N[node_shift + 1]) * sense[1];
56  double ksi20 = (N[node_shift + 0] - N[node_shift + 2]) * sense[2];
57  double L01[p[0] + 1], L12[p[1] + 1], L20[p[2] + 1];
58  double diffL01[2 * (p[0] + 1)], diffL12[2 * (p[1] + 1)],
59  diffL20[2 * (p[2] + 1)];
60  ierr = base_polynomials(p[0], ksi01, diff_ksi01, L01, diffL01, 2);
61  CHKERRQ(ierr);
62  ierr = base_polynomials(p[1], ksi12, diff_ksi12, L12, diffL12, 2);
63  CHKERRQ(ierr);
64  ierr = base_polynomials(p[2], ksi20, diff_ksi20, L20, diffL20, 2);
65  CHKERRQ(ierr);
66  int shift;
67  if (edgeN != NULL) {
68  // edge 01
69  shift = ii * (P[0]);
70  {
71  double *edge_n_ptr = &edgeN01[shift];
72  double *l_ptr = L01;
73  double scalar = N[node_shift + 0] * N[node_shift + 1];
74  int size = P[0];
75  for (size_t jj = 0; jj != size; ++jj, ++l_ptr) {
76  *edge_n_ptr = (*l_ptr) * scalar;
77  ++edge_n_ptr;
78  }
79  }
80 
81  // edge 12
82  shift = ii * (P[1]);
83  {
84  double *edge_n_ptr = &edgeN12[shift];
85  double *l_ptr = L12;
86  double scalar = N[node_shift + 1] * N[node_shift + 2];
87  int size = P[1];
88  for (size_t jj = 0; jj != size; ++jj, ++l_ptr) {
89  *edge_n_ptr = (*l_ptr) * scalar;
90  ++edge_n_ptr;
91  }
92  }
93 
94  // edge 20
95  shift = ii * (P[2]);
96  {
97  double *edge_n_ptr = &edgeN20[shift];
98  double *l_ptr = L20;
99  double scalar = N[node_shift + 2] * N[node_shift + 0];
100  int size = P[2];
101  for (size_t jj = 0; jj != size; ++jj, ++l_ptr) {
102  *edge_n_ptr = (*l_ptr) * scalar;
103  ++edge_n_ptr;
104  }
105  }
106  }
107  if (diff_edgeN != NULL) {
108  if (P[0] > 0) {
109  // edge01
110  shift = ii * (P[0]);
111  {
112  double *diff_edge_n_ptr = &diff_edgeN01[2 * shift];
113  double *diff_l_x = &diffL01[0 * (p[0] + 1)];
114  double *diff_l_y = &diffL01[1 * (p[0] + 1)];
115  double scalar_x = diffN[2 * 0 + 0] * N[node_shift + 1] +
116  N[node_shift + 0] * diffN[2 * 1 + 0];
117  double scalar_y = diffN[2 * 0 + 1] * N[node_shift + 1] +
118  N[node_shift + 0] * diffN[2 * 1 + 1];
119  double v = N[node_shift + 0] * N[node_shift + 1];
120  double *l_ptr = L01;
121 
122  int size = P[0];
123  for (size_t jj = 0; jj != size;
124  ++jj, ++diff_l_x, ++diff_l_y, ++l_ptr) {
125 
126  *diff_edge_n_ptr = v * (*diff_l_x) + scalar_x * (*l_ptr);
127  ++diff_edge_n_ptr;
128  *diff_edge_n_ptr = v * (*diff_l_y) + scalar_y * (*l_ptr);
129  ++diff_edge_n_ptr;
130  }
131  }
132 
133  }
134  if (P[1] > 0) {
135  // edge12
136  shift = ii * (P[1]);
137  {
138  double *diff_edge_n_ptr = &diff_edgeN12[2 * shift];
139  double *diff_l_x = &diffL12[0 * (p[1] + 1)];
140  double *diff_l_y = &diffL12[1 * (p[1] + 1)];
141  double scalar_x = diffN[2 * 1 + 0] * N[node_shift + 2] +
142  N[node_shift + 1] * diffN[2 * 2 + 0];
143  double scalar_y = diffN[2 * 1 + 1] * N[node_shift + 2] +
144  N[node_shift + 1] * diffN[2 * 2 + 1];
145  double v = N[node_shift + 1] * N[node_shift + 2];
146  double *l_ptr = L12;
147 
148  int size = P[1];
149  for (size_t jj = 0; jj != size;
150  ++jj, ++diff_l_x, ++diff_l_y, ++l_ptr) {
151 
152  *diff_edge_n_ptr = v * (*diff_l_x) + scalar_x * (*l_ptr);
153  ++diff_edge_n_ptr;
154  *diff_edge_n_ptr = v * (*diff_l_y) + scalar_y * (*l_ptr);
155  ++diff_edge_n_ptr;
156  }
157  }
158 
159  }
160  if (P[2] > 0) {
161  // edge20
162  shift = ii * (P[2]);
163 
164  {
165  double *diff_edge_n_ptr = &diff_edgeN20[2 * shift];
166  double *diff_l_x = &diffL20[0 * (p[2] + 1)];
167  double *diff_l_y = &diffL20[1 * (p[2] + 1)];
168  double scalar_x = diffN[2 * 2 + 0] * N[node_shift + 0] +
169  N[node_shift + 2] * diffN[2 * 0 + 0];
170  double scalar_y = diffN[2 * 2 + 1] * N[node_shift + 0] +
171  N[node_shift + 2] * diffN[2 * 0 + 1];
172  double v = N[node_shift + 2] * N[node_shift + 0];
173  double *l_ptr = L20;
174 
175  int size = P[2];
176  for (size_t jj = 0; jj != size;
177  ++jj, ++diff_l_x, ++diff_l_y, ++l_ptr) {
178 
179  *diff_edge_n_ptr = v * (*diff_l_x) + scalar_x * (*l_ptr);
180  ++diff_edge_n_ptr;
181  *diff_edge_n_ptr = v * (*diff_l_y) + scalar_y * (*l_ptr);
182  ++diff_edge_n_ptr;
183  }
184  }
185 
186  }
187  }
188  }
189  MoFEMFunctionReturnHot(0);
190 }

◆ H1_FaceGradientOfDeformation_hierarchical()

PetscErrorCode H1_FaceGradientOfDeformation_hierarchical ( int  p,
double *  diffN,
double *  dofs,
double *  F 
)

Definition at line 620 of file h1.c.

622  {
623  MoFEMFunctionBeginHot;
624  int col, row = 0;
625  for (; row < 3; row++)
626  for (col = 0; col < 3; col++)
627  F[3 * row + col] =
628  cblas_ddot(NBFACETRI_H1(p), &diffN[col], 3, &dofs[row], 3);
629  MoFEMFunctionReturnHot(0);
630 }

◆ H1_FaceShapeDiffMBTETinvJ()

PetscErrorCode H1_FaceShapeDiffMBTETinvJ ( int *  base_p,
int *  p,
double *  face_diffN[],
double *  invJac,
double *  face_diffNinvJac[],
int  GDIM 
)

Definition at line 574 of file h1.c.

576  {
577  MoFEMFunctionBeginHot;
578  int ff = 0, ii, gg;
579  for (; ff < 4; ff++) {
580  for (ii = 0; ii < NBFACETRI_H1(p[ff]); ii++) {
581  for (gg = 0; gg < GDIM; gg++) {
582  int shift1 = NBFACETRI_H1(base_p[ff]) * gg;
583  int shift2 = NBFACETRI_H1(p[ff]) * gg;
584  cblas_dgemv(CblasRowMajor, CblasTrans, 3, 3, 1., invJac, 3,
585  &(face_diffN[ff])[3 * shift1 + 3 * ii], 1, 0.,
586  &(face_diffNinvJac[ff])[3 * shift2 + 3 * ii], 1);
587  }
588  }
589  }
590  MoFEMFunctionReturnHot(0);
591 }

◆ H1_FaceShapeFunctions_MBTET()

PetscErrorCode H1_FaceShapeFunctions_MBTET ( int *  faces_nodes,
int *  p,
double *  N,
double *  diffN,
double *  faceN[],
double *  diff_faceN[],
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)  base_polynomials 
)

Definition at line 373 of file h1.c.

378  {
379  MoFEMFunctionBeginHot;
380 
381  int P[4];
382  int ff = 0;
383  for (; ff < 4; ff++)
384  P[ff] = NBFACETRI_H1(p[ff]);
385  double diff_ksiL0F0[3], diff_ksiL1F0[3];
386  double diff_ksiL0F1[3], diff_ksiL1F1[3];
387  double diff_ksiL0F2[3], diff_ksiL1F2[3];
388  double diff_ksiL0F3[3], diff_ksiL1F3[3];
389  double *diff_ksi_faces[] = {diff_ksiL0F0, diff_ksiL1F0, diff_ksiL0F1,
390  diff_ksiL1F1, diff_ksiL0F2, diff_ksiL1F2,
391  diff_ksiL0F3, diff_ksiL1F3};
392  int dd = 0;
393  for (; dd < 3; dd++) {
394  for (ff = 0; ff < 4; ff++) {
395  diff_ksi_faces[ff * 2 + 0][dd] =
396  (diffN[faces_nodes[3 * ff + 1] * 3 + dd] -
397  diffN[faces_nodes[3 * ff + 0] * 3 + dd]);
398  diff_ksi_faces[ff * 2 + 1][dd] =
399  (diffN[faces_nodes[3 * ff + 2] * 3 + dd] -
400  diffN[faces_nodes[3 * ff + 0] * 3 + dd]);
401  }
402  }
403  int ii = 0;
404  for (; ii < GDIM; ii++) {
405  int node_shift = ii * 4;
406  double ksi_faces[8];
407  for (ff = 0; ff < 4; ff++) {
408  ksi_faces[2 * ff + 0] = N[node_shift + faces_nodes[3 * ff + 1]] -
409  N[node_shift + faces_nodes[3 * ff + 0]];
410  ksi_faces[2 * ff + 1] = N[node_shift + faces_nodes[3 * ff + 2]] -
411  N[node_shift + faces_nodes[3 * ff + 0]];
412  }
413  int shift;
414  for (ff = 0; ff < 4; ff++) {
415  if (P[ff] == 0)
416  continue;
417  double L0[p[ff] + 1], L1[p[ff] + 1];
418  double diffL0[3 * (p[ff] + 1)], diffL1[3 * (p[ff] + 1)];
419  ierr = base_polynomials(p[ff], ksi_faces[ff * 2 + 0],
420  diff_ksi_faces[ff * 2 + 0], L0, diffL0, 3);
421  CHKERRQ(ierr);
422  ierr = base_polynomials(p[ff], ksi_faces[ff * 2 + 1],
423  diff_ksi_faces[ff * 2 + 1], L1, diffL1, 3);
424  CHKERRQ(ierr);
425  double v = N[node_shift + faces_nodes[3 * ff + 0]] *
426  N[node_shift + faces_nodes[3 * ff + 1]] *
427  N[node_shift + faces_nodes[3 * ff + 2]];
428  double v2[3] = {0, 0, 0};
429  dd = 0;
430  double n1n2 = N[node_shift + faces_nodes[3 * ff + 1]] *
431  N[node_shift + faces_nodes[3 * ff + 2]];
432  double n0n2 = N[node_shift + faces_nodes[3 * ff + 0]] *
433  N[node_shift + faces_nodes[3 * ff + 2]];
434  double n0n1 = N[node_shift + faces_nodes[3 * ff + 0]] *
435  N[node_shift + faces_nodes[3 * ff + 1]];
436  for (; dd < 3; dd++) {
437  v2[dd] = diffN[3 * faces_nodes[3 * ff + 0] + dd] * n1n2 +
438  diffN[3 * faces_nodes[3 * ff + 1] + dd] * n0n2 +
439  diffN[3 * faces_nodes[3 * ff + 2] + dd] * n0n1;
440  }
441  shift = ii * P[ff];
442  int jj = 0;
443  int oo = 0;
444  for (; oo <= (p[ff] - 3); oo++) {
445  int pp0 = 0;
446  for (; pp0 <= oo; pp0++) {
447  int pp1 = oo - pp0;
448  if (pp1 >= 0) {
449  if (faceN != NULL)
450  if (faceN[ff] != NULL) {
451  faceN[ff][shift + jj] = v * L0[pp0] * L1[pp1];
452  }
453  if (diff_faceN != NULL)
454  if (diff_faceN[ff] != NULL) {
455  dd = 0;
456  double L0L1 = L0[pp0] * L1[pp1];
457  for (; dd < 3; dd++) {
458  diff_faceN[ff][3 * shift + 3 * jj + dd] =
459  (L0[pp0] * diffL1[dd * (p[ff] + 1) + pp1] +
460  diffL0[dd * (p[ff] + 1) + pp0] * L1[pp1]) *
461  v;
462  diff_faceN[ff][3 * shift + 3 * jj + dd] += L0L1 * v2[dd];
463  }
464  }
465  jj++;
466  }
467  }
468  }
469  if (jj != P[ff])
470  SETERRQ2(PETSC_COMM_SELF, 1, "wrong order %d != %d", jj, P[ff]);
471  }
472  }
473  MoFEMFunctionReturnHot(0);
474 }

◆ H1_FaceShapeFunctions_MBTRI()

PetscErrorCode H1_FaceShapeFunctions_MBTRI ( const int *  face_nodes,
int  p,
double *  N,
double *  diffN,
double *  faceN,
double *  diff_faceN,
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)  base_polynomials 
)

Definition at line 191 of file h1.c.

196  {
197  MoFEMFunctionBeginHot;
198 
199  int P = NBFACETRI_H1(p);
200  if (P == 0)
201  MoFEMFunctionReturnHot(0);
202  double diff_ksiL0[2], diff_ksiL1[2];
203  double *diff_ksi_faces[] = {diff_ksiL0, diff_ksiL1};
204  int dd = 0;
205  if (diff_faceN != NULL) {
206  for (; dd < 2; dd++) {
207  diff_ksi_faces[0][dd] =
208  diffN[face_nodes[1] * 2 + dd] - diffN[face_nodes[0] * 2 + dd];
209  diff_ksi_faces[1][dd] =
210  diffN[face_nodes[2] * 2 + dd] - diffN[face_nodes[0] * 2 + dd];
211  }
212  }
213  int ii = 0;
214  for (; ii < GDIM; ii++) {
215  int node_shift = ii * 3;
216  double ksi_faces[2];
217  ksi_faces[0] =
218  N[node_shift + face_nodes[1]] - N[node_shift + face_nodes[0]];
219  ksi_faces[1] =
220  N[node_shift + face_nodes[2]] - N[node_shift + face_nodes[0]];
221  double L0[p + 1], L1[p + 1];
222  double diffL0[2 * (p + 1)], diffL1[2 * (p + 1)];
223  ierr = base_polynomials(p, ksi_faces[0], diff_ksi_faces[0], L0, diffL0, 2);
224  CHKERRQ(ierr);
225  ierr = base_polynomials(p, ksi_faces[1], diff_ksi_faces[1], L1, diffL1, 2);
226  CHKERRQ(ierr);
227  double v = N[node_shift + face_nodes[0]] * N[node_shift + face_nodes[1]] *
228  N[node_shift + face_nodes[2]];
229  double v2[2] = {0, 0};
230  if (diff_faceN != NULL) {
231  double n1n2 =
232  N[node_shift + face_nodes[1]] * N[node_shift + face_nodes[2]];
233  double n0n2 =
234  N[node_shift + face_nodes[0]] * N[node_shift + face_nodes[2]];
235  double n0n1 =
236  N[node_shift + face_nodes[0]] * N[node_shift + face_nodes[1]];
237  dd = 0;
238  for (; dd < 2; dd++) {
239  v2[dd] = diffN[face_nodes[0] * 2 + dd] * n1n2 +
240  diffN[face_nodes[1] * 2 + dd] * n0n2 +
241  diffN[face_nodes[2] * 2 + dd] * n0n1;
242  }
243  }
244  int shift = ii * P;
245  int jj = 0;
246  int oo = 0;
247  for (; oo <= (p - 3); oo++) {
248  int pp0 = 0;
249  for (; pp0 <= oo; pp0++) {
250  int pp1 = oo - pp0;
251  if (pp1 >= 0) {
252  if (faceN != NULL) {
253  faceN[shift + jj] = v * L0[pp0] * L1[pp1];
254  }
255  if (diff_faceN != NULL) {
256  dd = 0;
257  for (; dd < 2; dd++) {
258  double *val = &diff_faceN[2 * shift + 2 * jj + dd];
259  *val = (L0[pp0] * diffL1[dd * (p + 1) + pp1] +
260  diffL0[dd * (p + 1) + pp0] * L1[pp1]) *
261  v;
262  *val += L0[pp0] * L1[pp1] * v2[dd];
263  }
264  }
265  jj++;
266  }
267  }
268  }
269  if (jj != P)
270  SETERRQ2(PETSC_COMM_SELF, 1, "wrong order %d != %d", jj, P);
271  }
272  MoFEMFunctionReturnHot(0);
273 }

◆ H1_QuadShapeFunctions_MBPRISM()

PetscErrorCode H1_QuadShapeFunctions_MBPRISM ( int *  faces_nodes,
int *  p,
double *  N,
double *  diffN,
double *  faceN[],
double *  diff_faceN[],
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)  base_polynomials 
)

Definition at line 642 of file h1.c.

647  {
648  MoFEMFunctionBeginHot;
649  // TODO: return separately components of the tensor product between two edges
650 
651  int P[3];
652  int ff = 0;
653  for (; ff < 3; ff++)
654  P[ff] = NBFACEQUAD_H1(p[ff]);
655 
656  double ksi_faces[6];
657  double diff_ksiL0F0[3], diff_ksiL3F0[3];
658  double diff_ksiL0F1[3], diff_ksiL3F1[3];
659  double diff_ksiL0F2[3], diff_ksiL3F2[3];
660  double *diff_ksi_faces[] = {diff_ksiL0F0, diff_ksiL3F0, diff_ksiL0F1,
661  diff_ksiL3F1, diff_ksiL0F2, diff_ksiL3F2};
662  int ii = 0;
663  for (; ii < GDIM; ii++) {
664  int node_shift = ii * 6;
665  int node_diff_shift = 3 * node_shift;
666  int ff = 0;
667  for (; ff < 3; ff++) {
668  if (P[ff] == 0)
669  continue;
670  int n0 = faces_nodes[4 * ff + 0];
671  int n1 = faces_nodes[4 * ff + 1];
672  int n2 = faces_nodes[4 * ff + 2];
673  int n3 = faces_nodes[4 * ff + 3];
674  int e0 = 2 * ff + 0;
675  int e1 = 2 * ff + 1;
676  double ksi0 = N[node_shift + n0] + N[node_shift + n3];
677  double ksi1 = N[node_shift + n1] + N[node_shift + n2];
678  double eta0 = N[node_shift + n0] + N[node_shift + n1];
679  double eta1 = N[node_shift + n2] + N[node_shift + n3];
680  ksi_faces[e0] = ksi1 - ksi0;
681  ksi_faces[e1] = eta1 - eta0;
682 
683  int dd = 0;
684  for (; dd < 3; dd++) {
685  double diff_ksi0 = diffN[node_diff_shift + 3 * n0 + dd] +
686  diffN[node_diff_shift + 3 * n3 + dd];
687  double diff_ksi1 = diffN[node_diff_shift + 3 * n1 + dd] +
688  diffN[node_diff_shift + 3 * n2 + dd];
689  double diff_eta0 = diffN[node_diff_shift + 3 * n0 + dd] +
690  diffN[node_diff_shift + 3 * n1 + dd];
691  double diff_eta1 = diffN[node_diff_shift + 3 * n2 + dd] +
692  diffN[node_diff_shift + 3 * n3 + dd];
693  diff_ksi_faces[e0][dd] = diff_ksi1 - diff_ksi0;
694  diff_ksi_faces[e1][dd] = diff_eta1 - diff_eta0;
695  }
696  double L0[p[ff] + 1], L1[p[ff] + 1];
697  double diffL0[3 * (p[ff] + 1)], diffL1[3 * (p[ff] + 1)];
698  ierr = base_polynomials(p[ff], ksi_faces[e0], diff_ksi_faces[e0], L0,
699  diffL0, 3);
700  CHKERRQ(ierr);
701  ierr = base_polynomials(p[ff], ksi_faces[e1], diff_ksi_faces[e1], L1,
702  diffL1, 3);
703  CHKERRQ(ierr);
704 
705  double v = N[node_shift + n0] * N[node_shift + n2] +
706  N[node_shift + n1] * N[node_shift + n3];
707 
708  double diff_v[3];
709  dd = 0;
710  for (; dd < 3; ++dd)
711  diff_v[dd] =
712 
713  diffN[node_diff_shift + 3 * n0 + dd] * N[node_shift + n2] +
714  N[node_shift + n0] * diffN[node_diff_shift + 3 * n2 + dd] +
715 
716  diffN[node_diff_shift + 3 * n1 + dd] * N[node_shift + n3] +
717  N[node_shift + n1] * diffN[node_diff_shift + 3 * n3 + dd];
718 
719  int shift;
720  shift = ii * P[ff];
721 
722  if (faceN != NULL) {
723  if (faceN[ff] != NULL) {
724  int jj = 0;
725  int oo = 0;
726  for (; oo <= p[ff] - 2; ++oo) {
727  int pp0 = 0;
728  for (; pp0 < oo; ++pp0) {
729  int pp1 = oo;
730  faceN[ff][shift + jj] = L0[pp0] * L1[pp1] * v;
731  ++jj;
732  faceN[ff][shift + jj] = L0[pp1] * L1[pp0] * v;
733  ++jj;
734  }
735  faceN[ff][shift + jj] = L0[oo] * L1[oo] * v;
736  ++jj;
737  }
738  if (jj != P[ff])
739  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
740  "Inconsistent implementation (bug in the code) %d != %d",
741  jj, P);
742  }
743  }
744 
745  if (diff_faceN != NULL) {
746  if (diff_faceN[ff] != NULL) {
747  int jj = 0;
748  int oo = 0;
749  for (; oo <= p[ff] - 2; ++oo) {
750  int pp0 = 0;
751  for (; pp0 < oo; ++pp0) {
752  int pp1 = oo;
753  int dd;
754  for (dd = 0; dd < 3; dd++) {
755  diff_faceN[ff][3 * shift + 3 * jj + dd] =
756  (L0[pp0] * diffL1[dd * (p[ff] + 1) + pp1] +
757  diffL0[dd * (p[ff] + 1) + pp0] * L1[pp1]) *
758  v +
759  L0[pp0] * L1[pp1] * diff_v[dd];
760  }
761  ++jj;
762  for (dd = 0; dd < 3; dd++) {
763  diff_faceN[ff][3 * shift + 3 * jj + dd] =
764  (L0[pp1] * diffL1[dd * (p[ff] + 1) + pp0] +
765  diffL0[dd * (p[ff] + 1) + pp1] * L1[pp0]) *
766  v +
767  L0[pp1] * L1[pp0] * diff_v[dd];
768  }
769  ++jj;
770  }
771  for (dd = 0; dd < 3; dd++) {
772  diff_faceN[ff][3 * shift + 3 * jj + dd] =
773  (L0[oo] * diffL1[dd * (p[ff] + 1) + oo] +
774  diffL0[dd * (p[ff] + 1) + oo] * L1[oo]) *
775  v +
776  L0[oo] * L1[oo] * diff_v[dd];
777  }
778  ++jj;
779  }
780  if (jj != P[ff])
781  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
782  "Inconsistent implementation (bug in the code) %d != %d",
783  jj, P);
784  }
785  }
786  }
787  }
788  MoFEMFunctionReturnHot(0);
789 }

◆ H1_QuadShapeFunctions_MBQUAD()

PetscErrorCode H1_QuadShapeFunctions_MBQUAD ( int *  faces_nodes,
int  p,
double *  N,
double *  diffN,
double *  faceN,
double *  diff_faceN,
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)  base_polynomials 
)

Definition at line 959 of file h1.c.

964  {
965  MoFEMFunctionBeginHot;
966  // TODO: return separately components of the tensor product between two edges
967 
968  int P = NBFACEQUAD_H1(p);
969  if (P == 0)
970  MoFEMFunctionReturnHot(0);
971  double diff_ksiL0F0[2], diff_ksiL1F0[2];
972  double *diff_ksi_faces[] = {diff_ksiL0F0, diff_ksiL1F0};
973  int ii = 0;
974  for (; ii < GDIM; ii++) {
975  int node_shift = ii * 4;
976  int node_diff_shift = 2 * node_shift;
977 
978  int n0 = faces_nodes[0];
979  int n1 = faces_nodes[1];
980  int n2 = faces_nodes[2];
981  int n3 = faces_nodes[3];
982  double ksi0 = N[node_shift + n0] + N[node_shift + n3];
983  double ksi1 = N[node_shift + n1] + N[node_shift + n2];
984  double eta0 = N[node_shift + n0] + N[node_shift + n1];
985  double eta1 = N[node_shift + n2] + N[node_shift + n3];
986  double ksi_faces = ksi1 - ksi0;
987  double eta_faces = eta1 - eta0;
988  double diff_ksi_faces[2];
989  double diff_eta_faces[2];
990  int dd = 0;
991  for (; dd < 2; dd++) {
992  double diff_ksi0 = diffN[node_diff_shift + 2 * n0 + dd] +
993  diffN[node_diff_shift + 2 * n3 + dd];
994  double diff_ksi1 = diffN[node_diff_shift + 2 * n1 + dd] +
995  diffN[node_diff_shift + 2 * n2 + dd];
996  double diff_eta0 = diffN[node_diff_shift + 2 * n0 + dd] +
997  diffN[node_diff_shift + 2 * n1 + dd];
998  double diff_eta1 = diffN[node_diff_shift + 2 * n2 + dd] +
999  diffN[node_diff_shift + 2 * n3 + dd];
1000 
1001  diff_ksi_faces[dd] = diff_ksi1 - diff_ksi0;
1002  diff_eta_faces[dd] = diff_eta1 - diff_eta0;
1003  }
1004  double L0[p + 1], L1[p + 1];
1005  double diffL0[2 * (p + 1)], diffL1[2 * (p + 1)];
1006  ierr = base_polynomials(p, ksi_faces, diff_ksi_faces, L0, diffL0, 2);
1007  CHKERRQ(ierr);
1008  ierr = base_polynomials(p, eta_faces, diff_eta_faces, L1, diffL1, 2);
1009  CHKERRQ(ierr);
1010 
1011  double v = N[node_shift + n0] * N[node_shift + n2] +
1012  N[node_shift + n1] * N[node_shift + n3];
1013 
1014  double diff_v[2];
1015  dd = 0;
1016  for (; dd < 2; ++dd)
1017  diff_v[dd] =
1018 
1019  diffN[node_diff_shift + 2 * n0 + dd] * N[node_shift + n2] +
1020  N[node_shift + n0] * diffN[node_diff_shift + 2 * n2 + dd] +
1021 
1022  diffN[node_diff_shift + 2 * n1 + dd] * N[node_shift + n3] +
1023  N[node_shift + n1] * diffN[node_diff_shift + 2 * n3 + dd];
1024 
1025  const int shift = ii * P;
1026 
1027  if (faceN != NULL) {
1028  int jj = 0;
1029  int oo = 0;
1030  for (; oo <= p - 2; ++oo) {
1031  int pp0 = 0;
1032  for (; pp0 < oo; ++pp0) {
1033  int pp1 = oo;
1034  faceN[shift + jj] = L0[pp0] * L1[pp1] * v;
1035  ++jj;
1036  faceN[shift + jj] = L0[pp1] * L1[pp0] * v;
1037  ++jj;
1038  }
1039  faceN[shift + jj] = L0[oo] * L1[oo] * v;
1040  ++jj;
1041  }
1042  if (jj != P)
1043  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1044  "Inconsistent implementation (bug in the code) %d != %d", jj,
1045  P);
1046  }
1047 
1048  if (diff_faceN != NULL) {
1049  int jj = 0;
1050  int oo = 0;
1051  for (; oo <= p - 2; ++oo) {
1052  int pp0 = 0;
1053  for (; pp0 < oo; ++pp0) {
1054  int pp1 = oo;
1055  int dd;
1056  for (dd = 0; dd < 2; dd++) {
1057  diff_faceN[2 * shift + 2 * jj + dd] =
1058  (L0[pp0] * diffL1[dd * (p + 1) + pp1] +
1059  diffL0[dd * (p + 1) + pp0] * L1[pp1]) *
1060  v +
1061  L0[pp0] * L1[pp1] * diff_v[dd];
1062  }
1063  ++jj;
1064  for (dd = 0; dd < 2; dd++) {
1065  diff_faceN[2 * shift + 2 * jj + dd] =
1066  (L0[pp1] * diffL1[dd * (p + 1) + pp0] +
1067  diffL0[dd * (p + 1) + pp1] * L1[pp0]) *
1068  v +
1069  L0[pp1] * L1[pp0] * diff_v[dd];
1070  }
1071  ++jj;
1072  }
1073  for (dd = 0; dd < 2; dd++) {
1074  diff_faceN[2 * shift + 2 * jj + dd] =
1075  (L0[oo] * diffL1[dd * (p + 1) + oo] +
1076  diffL0[dd * (p + 1) + oo] * L1[oo]) *
1077  v +
1078  L0[oo] * L1[oo] * diff_v[dd];
1079  }
1080  ++jj;
1081  }
1082  if (jj != P)
1083  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1084  "Inconsistent implementation (bug in the code) %d != %d", jj,
1085  P);
1086  }
1087  }
1088  MoFEMFunctionReturnHot(0);
1089 }

◆ H1_VolumeGradientOfDeformation_hierarchical()

PetscErrorCode H1_VolumeGradientOfDeformation_hierarchical ( int  p,
double *  diffN,
double *  dofs,
double *  F 
)

Definition at line 631 of file h1.c.

633  {
634  MoFEMFunctionBeginHot;
635  int col, row = 0;
636  for (; row < 3; row++)
637  for (col = 0; col < 3; col++)
638  F[3 * row + col] =
639  cblas_ddot(NBVOLUMETET_H1(p), &diffN[col], 3, &dofs[row], 3);
640  MoFEMFunctionReturnHot(0);
641 }

◆ H1_VolumeShapeDiffMBTETinvJ()

PetscErrorCode H1_VolumeShapeDiffMBTETinvJ ( int  base_p,
int  p,
double *  volume_diffN,
double *  invJac,
double *  volume_diffNinvJac,
int  GDIM 
)

Definition at line 592 of file h1.c.

595  {
596  MoFEMFunctionBeginHot;
597  int ii, gg;
598  for (ii = 0; ii < NBVOLUMETET_H1(p); ii++) {
599  for (gg = 0; gg < GDIM; gg++) {
600  int shift1 = NBVOLUMETET_H1(base_p) * gg;
601  int shift2 = NBVOLUMETET_H1(p) * gg;
602  cblas_dgemv(CblasRowMajor, CblasTrans, 3, 3, 1., invJac, 3,
603  &(volume_diffN)[3 * shift1 + 3 * ii], 1, 0.,
604  &(volume_diffNinvJac)[3 * shift2 + 3 * ii], 1);
605  }
606  }
607  MoFEMFunctionReturnHot(0);
608 }

◆ H1_VolumeShapeFunctions_MBPRISM()

PetscErrorCode H1_VolumeShapeFunctions_MBPRISM ( int  p,
double *  N,
double *  diffN,
double *  volumeN,
double *  diff_volumeN,
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)  base_polynomials 
)

Definition at line 790 of file h1.c.

795  {
796  MoFEMFunctionBeginHot;
797 
798  int P = NBVOLUMEPRISM_H1(p);
799  if (P == 0)
800  MoFEMFunctionReturnHot(0);
801  double diff_ksiL0[3], diff_ksiL1[3], diff_ksiL2[3];
802  int ii = 0;
803  for (; ii < GDIM; ii++) {
804  int node_shift = ii * 6;
805  int node_diff_shift = ii * 18;
806 
807  double n0 = N[node_shift + 0];
808  double n1 = N[node_shift + 1];
809  double n2 = N[node_shift + 2];
810  double n3 = N[node_shift + 3];
811  double n4 = N[node_shift + 4];
812  double n5 = N[node_shift + 5];
813 
814  double ksiL0 = n1 + n4 - n0 - n3;
815  double ksiL1 = n2 + n5 - n0 - n3;
816  double ksiL2 = (n3 + n4 + n5) - (n0 + n1 + n2);
817 
818  int dd = 0;
819  for (; dd < 3; dd++) {
820  double diff_n0 = diffN[node_diff_shift + 3 * 0 + dd];
821  double diff_n1 = diffN[node_diff_shift + 3 * 1 + dd];
822  double diff_n2 = diffN[node_diff_shift + 3 * 2 + dd];
823  double diff_n3 = diffN[node_diff_shift + 3 * 3 + dd];
824  double diff_n4 = diffN[node_diff_shift + 3 * 4 + dd];
825  double diff_n5 = diffN[node_diff_shift + 3 * 5 + dd];
826  diff_ksiL0[dd] = (diff_n1 + diff_n4) - (diff_n0 + diff_n3);
827  diff_ksiL1[dd] = (diff_n2 + diff_n5) - (diff_n0 + diff_n3);
828  diff_ksiL2[dd] =
829  (diff_n3 + diff_n4 + diff_n5) - (diff_n0 + diff_n1 + diff_n2);
830  }
831 
832  double L0[p + 1], L1[p + 1], L2[p + 1];
833  double diffL0[3 * (p + 1)], diffL1[3 * (p + 1)], diffL2[3 * (p + 1)];
834  ierr = base_polynomials(p, ksiL0, diff_ksiL0, L0, diffL0, 3);
835  CHKERRQ(ierr);
836  ierr = base_polynomials(p, ksiL1, diff_ksiL1, L1, diffL1, 3);
837  CHKERRQ(ierr);
838  ierr = base_polynomials(p, ksiL2, diff_ksiL2, L2, diffL2, 3);
839  CHKERRQ(ierr);
840 
841  double v_tri = (n0 + n3) * (n1 + n4) * (n2 + n5);
842  double v_edge = (n0 + n1 + n2) * (n3 + n4 + n5);
843  double v = v_tri * v_edge;
844 
845  double diff_v_tri[3];
846  double diff_v_edge[3];
847  double diff_v[3];
848  dd = 0;
849  for (; dd < 3; dd++) {
850  double diff_n0 = diffN[node_diff_shift + 3 * 0 + dd];
851  double diff_n1 = diffN[node_diff_shift + 3 * 1 + dd];
852  double diff_n2 = diffN[node_diff_shift + 3 * 2 + dd];
853  double diff_n3 = diffN[node_diff_shift + 3 * 3 + dd];
854  double diff_n4 = diffN[node_diff_shift + 3 * 4 + dd];
855  double diff_n5 = diffN[node_diff_shift + 3 * 5 + dd];
856  diff_v_tri[dd] = (diff_n0 + diff_n3) * (n1 + n4) * (n2 + n5) +
857  (diff_n1 + diff_n4) * (n0 + n3) * (n2 + n5) +
858  (diff_n2 + diff_n5) * (n0 + n3) * (n1 + n4);
859  diff_v_edge[dd] = (diff_n0 + diff_n1 + diff_n2) * (n3 + n4 + n5) +
860  (diff_n3 + diff_n4 + diff_n5) * (n0 + n1 + n2);
861  diff_v[dd] = diff_v_tri[dd] * v_edge + v_tri * diff_v_edge[dd];
862  }
863 
864  int shift = ii * P;
865 
866  if (volumeN != NULL) {
867  int jj = 0;
868  int oo, pp0, pp1, zz;
869  for (oo = 0; oo <= p - 3; ++oo) {
870  for (pp0 = 0; pp0 < oo; ++pp0) {
871  for (pp1 = 0; pp1 < oo; ++pp1) {
872  const int perm[3][3] = {
873  {oo, pp0, pp1}, {pp0, oo, pp1}, {pp0, pp1, oo}};
874  for (zz = 0; zz != 3; ++zz) {
875  volumeN[shift + jj] =
876  L0[perm[zz][0]] * L1[perm[zz][1]] * L2[perm[zz][2]] * v;
877  ++jj;
878  }
879  }
880  const int perm[3][3] = {{pp0, oo, oo}, {oo, pp0, oo}, {oo, oo, pp0}};
881  for (zz = 0; zz != 3; ++zz) {
882  volumeN[shift + jj] =
883  L0[perm[zz][0]] * L1[perm[zz][1]] * L2[perm[zz][2]] * v;
884  ++jj;
885  }
886  }
887  volumeN[shift + jj] = L0[oo] * L1[oo] * L2[oo] * v;
888  ++jj;
889  }
890  if (jj != P)
891  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
892  "Inconsistent implementation (bug in the code) %d != %d", jj,
893  P);
894  }
895 
896  if (diff_volumeN != NULL) {
897  int jj = 0;
898  int oo, pp0, pp1, zz, dd;
899  for (oo = 0; oo <= p - 3; ++oo) {
900  for (pp0 = 0; pp0 < oo; ++pp0) {
901  for (pp1 = 0; pp1 < oo; ++pp1) {
902  const int perm[3][3] = {
903  {oo, pp0, pp1}, {pp0, oo, pp1}, {pp0, pp1, oo}};
904  for (zz = 0; zz != 3; ++zz) {
905  const int i = perm[zz][0];
906  const int j = perm[zz][1];
907  const int k = perm[zz][2];
908  for (dd = 0; dd != 3; ++dd) {
909  diff_volumeN[3 * shift + 3 * jj + dd] =
910  (diffL0[dd * (p + 1) + i] * L1[j] * L2[k] +
911  L0[i] * diffL1[dd * (p + 1) + j] * L2[k] +
912  L0[i] * L1[j] * diffL2[dd * (p + 1) + k]) *
913  v +
914  L0[i] * L1[j] * L2[k] * diff_v[dd];
915  }
916  ++jj;
917  }
918  }
919  const int perm[3][3] = {{pp0, oo, oo}, {oo, pp0, oo}, {oo, oo, pp0}};
920  for (zz = 0; zz != 3; ++zz) {
921  const int i = perm[zz][0];
922  const int j = perm[zz][1];
923  const int k = perm[zz][2];
924  for (dd = 0; dd != 3; ++dd) {
925  diff_volumeN[3 * shift + 3 * jj + dd] =
926  (diffL0[dd * (p + 1) + i] * L1[j] * L2[k] +
927  L0[i] * diffL1[dd * (p + 1) + j] * L2[k] +
928  L0[i] * L1[j] * diffL2[dd * (p + 1) + k]) *
929  v +
930  L0[i] * L1[j] * L2[k] * diff_v[dd];
931  }
932  ++jj;
933  }
934  }
935 
936  const int i = oo;
937  const int j = oo;
938  const int k = oo;
939  for (dd = 0; dd != 3; ++dd) {
940  diff_volumeN[3 * shift + 3 * jj + dd] =
941  (diffL0[dd * (p + 1) + i] * L1[j] * L2[k] +
942  L0[i] * diffL1[dd * (p + 1) + j] * L2[k] +
943  L0[i] * L1[j] * diffL2[dd * (p + 1) + k]) *
944  v +
945  L0[i] * L1[j] * L2[k] * diff_v[dd];
946  }
947 
948  ++jj;
949  }
950  if (jj != P)
951  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
952  "Inconsistent implementation (bug in the code) %d != %d", jj,
953  P);
954  }
955  }
956  MoFEMFunctionReturnHot(0);
957 }

◆ H1_VolumeShapeFunctions_MBTET()

PetscErrorCode H1_VolumeShapeFunctions_MBTET ( int  p,
double *  N,
double *  diffN,
double *  volumeN,
double *  diff_volumeN,
int  GDIM,
PetscErrorCode(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)  base_polynomials 
)

Definition at line 475 of file h1.c.

480  {
481  MoFEMFunctionBeginHot;
482 
483  int P = NBVOLUMETET_H1(p);
484  if (P == 0)
485  MoFEMFunctionReturnHot(0);
486  double diff_ksiL0[3], diff_ksiL1[3], diff_ksiL2[3];
487  int dd = 0;
488  for (; dd < 3; dd++) {
489  diff_ksiL0[dd] = (diffN[1 * 3 + dd] - diffN[0 * 3 + dd]);
490  diff_ksiL1[dd] = (diffN[2 * 3 + dd] - diffN[0 * 3 + dd]);
491  diff_ksiL2[dd] = (diffN[3 * 3 + dd] - diffN[0 * 3 + dd]);
492  }
493  int ii = 0;
494  for (; ii < GDIM; ii++) {
495  int node_shift = ii * 4;
496  double ksiL0 = N[node_shift + 1] - N[node_shift + 0];
497  double ksiL1 = N[node_shift + 2] - N[node_shift + 0];
498  double ksiL2 = N[node_shift + 3] - N[node_shift + 0];
499  double L0[p + 1], L1[p + 1], L2[p + 1];
500  double diffL0[3 * (p + 1)], diffL1[3 * (p + 1)], diffL2[3 * (p + 1)];
501  ierr = base_polynomials(p, ksiL0, diff_ksiL0, L0, diffL0, 3);
502  CHKERRQ(ierr);
503  ierr = base_polynomials(p, ksiL1, diff_ksiL1, L1, diffL1, 3);
504  CHKERRQ(ierr);
505  ierr = base_polynomials(p, ksiL2, diff_ksiL2, L2, diffL2, 3);
506  CHKERRQ(ierr);
507  double v = N[node_shift + 0] * N[node_shift + 1] * N[node_shift + 2] *
508  N[node_shift + 3];
509  double v2[3] = {0, 0, 0};
510  dd = 0;
511  for (; dd < 3; dd++) {
512  v2[dd] = diffN[3 * 0 + dd] * N[node_shift + 1] * N[node_shift + 2] *
513  N[node_shift + 3] +
514  N[node_shift + 0] * diffN[3 * 1 + dd] * N[node_shift + 2] *
515  N[node_shift + 3] +
516  N[node_shift + 0] * N[node_shift + 1] * diffN[3 * 2 + dd] *
517  N[node_shift + 3] +
518  N[node_shift + 0] * N[node_shift + 1] * N[node_shift + 2] *
519  diffN[3 * 3 + dd];
520  }
521  int shift = ii * P;
522  int jj = 0;
523  int oo = 0;
524  for (; oo <= (p - 4); oo++) {
525  int pp0 = 0;
526  for (; pp0 <= oo; pp0++) {
527  int pp1 = 0;
528  for (; (pp0 + pp1) <= oo; pp1++) {
529  int pp2 = oo - pp0 - pp1;
530  if (pp2 >= 0) {
531  if (volumeN != NULL) {
532  volumeN[shift + jj] = L0[pp0] * L1[pp1] * L2[pp2] * v;
533  }
534  if (diff_volumeN != NULL) {
535  dd = 0;
536  for (; dd < 3; dd++) {
537  diff_volumeN[3 * shift + 3 * jj + dd] =
538  (diffL0[dd * (p + 1) + pp0] * L1[pp1] * L2[pp2] +
539  L0[pp0] * diffL1[dd * (p + 1) + pp1] * L2[pp2] +
540  L0[pp0] * L1[pp1] * diffL2[dd * (p + 1) + pp2]) *
541  v;
542  diff_volumeN[3 * shift + 3 * jj + dd] +=
543  L0[pp0] * L1[pp1] * L2[pp2] * v2[dd];
544  }
545  }
546  jj++;
547  }
548  }
549  }
550  }
551  if (jj != P)
552  SETERRQ1(PETSC_COMM_SELF, 1, "wrong order %d", jj);
553  }
554  MoFEMFunctionReturnHot(0);
555 }

Variable Documentation

◆ ierr

PetscErrorCode ierr
static

Definition at line 15 of file h1.c.

MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
NBEDGE_H1
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
Definition: h1_hdiv_hcurl_l2.h:55
sdf_hertz.d
float d
Definition: sdf_hertz.py:5
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
NBVOLUMETET_H1
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
Definition: h1_hdiv_hcurl_l2.h:75
EshelbianPlasticity::P
@ P
Definition: EshelbianContact.cpp:193
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
NBFACEQUAD_H1
#define NBFACEQUAD_H1(P)
Number of base functions on quad for H1 space.
Definition: h1_hdiv_hcurl_l2.h:65
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
N
const int N
Definition: speed_test.cpp:3
ierr
static PetscErrorCode ierr
Definition: h1.c:15
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
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
NBFACETRI_H1
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
Definition: h1_hdiv_hcurl_l2.h:60
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
NBVOLUMEPRISM_H1
#define NBVOLUMEPRISM_H1(P)
Number of base functions on prism for H1 space.
Definition: h1_hdiv_hcurl_l2.h:80
invJac
MatrixDouble invJac
Definition: HookeElement.hpp:683
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
F
@ F
Definition: free_surface.cpp:394
Generated by Doxygen 1.8.17 and hosted at University of Glasgow