v0.14.0
Loading...
Searching...
No Matches
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.
 
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}
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:460
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:453
F
@ F
Definition free_surface.cpp:394
NBEDGE_H1
#define NBEDGE_H1(P)
Number of base function on edge for H1 space.
Definition h1_hdiv_hcurl_l2.h:55

◆ 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}
invJac
MatrixDouble invJac
Definition HookeElement.hpp:683

◆ 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(* base_polynomials )(int p, double s, double *diff_s, double *L, double *diffL, const int dim) )

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}
ierr
static PetscErrorCode ierr
Definition h1.c:15
sdf_hertz_2d_axisymm_plane.d
float d
Definition sdf_hertz_2d_axisymm_plane.py:4
N
const int N
Definition speed_test.cpp:3

◆ H1_EdgeShapeFunctions_MBTET()

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) )

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}
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

◆ 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(* base_polynomials )(int p, double s, double *diff_s, double *L, double *diffL, const int dim) )

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}
NBFACETRI_H1
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
Definition h1_hdiv_hcurl_l2.h:60

◆ 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(* base_polynomials )(int p, double s, double *diff_s, double *L, double *diffL, const int dim) )

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(* base_polynomials )(int p, double s, double *diff_s, double *L, double *diffL, const int dim) )

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(* base_polynomials )(int p, double s, double *diff_s, double *L, double *diffL, const int dim) )

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}
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
NBFACEQUAD_H1
#define NBFACEQUAD_H1(P)
Number of base functions on quad for H1 space.
Definition h1_hdiv_hcurl_l2.h:65

◆ H1_QuadShapeFunctions_MBQUAD()

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) )

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}
NBVOLUMETET_H1
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
Definition h1_hdiv_hcurl_l2.h:75

◆ 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(* base_polynomials )(int p, double s, double *diff_s, double *L, double *diffL, const int dim) )

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}
L2
@ L2
field with C-1 continuity
Definition definitions.h:88
NBVOLUMEPRISM_H1
#define NBVOLUMEPRISM_H1(P)
Number of base functions on prism for H1 space.
Definition h1_hdiv_hcurl_l2.h:80
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:19
k
FTensor::Index< 'k', 3 > k
Definition matrix_function.cpp:20

◆ H1_VolumeShapeFunctions_MBTET()

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) )

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.

Generated by Doxygen 1.12.0 and hosted at University of Glasgow