v0.9.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))
 
PetscErrorCode H1_EdgeGradientOfDeformation_hierachical (int p, double *diffN, double *dofs, double *F)
 
PetscErrorCode H1_FaceGradientOfDeformation_hierachical (int p, double *diffN, double *dofs, double *F)
 
PetscErrorCode H1_VolumeGradientOfDeformation_hierachical (int p, double *diffN, double *dofs, double *F)
 

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

PetscErrorCode H1_EdgeGradientOfDeformation_hierachical ( int  p,
double diffN,
double dofs,
double F 
)
Deprecated:
use H1_EdgeGradientOfDeformation_hierarchical

Definition at line 1122 of file h1.c.

1124  {
1125  return H1_EdgeGradientOfDeformation_hierarchical(p, diffN, dofs, F);
1126 }
PetscErrorCode H1_EdgeGradientOfDeformation_hierarchical(int p, double *diffN, double *dofs, double *F)
Definition: h1.c:553

◆ H1_EdgeGradientOfDeformation_hierarchical()

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

Definition at line 553 of file h1.c.

555  {
557  int col, row = 0;
558  for (; row < 3; row++)
559  for (col = 0; col < 3; col++)
560  F[3 * row + col] =
561  cblas_ddot(NBEDGE_H1(p), &diffN[col], 3, &dofs[row], 3);
563 }
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: cblas_ddot.c:12

◆ H1_EdgeShapeDiffMBTETinvJ()

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

Definition at line 500 of file h1.c.

502  {
504  int ee = 0, ii, gg;
505  for (; ee < 6; ee++) {
506  for (ii = 0; ii < NBEDGE_H1(p[ee]); ii++) {
507  for (gg = 0; gg < GDIM; gg++) {
508  int shift1 = NBEDGE_H1(base_p[ee]) * gg;
509  int shift2 = NBEDGE_H1(p[ee]) * gg;
510  cblas_dgemv(CblasRowMajor, CblasTrans, 3, 3, 1., invJac, 3,
511  &(edge_diffN[ee])[3 * shift1 + 3 * ii], 1, 0.,
512  &(edge_diffNinvJac[ee])[3 * shift2 + 3 * ii], 1);
513  }
514  }
515  }
517 }
void cblas_dgemv(const enum CBLAS_ORDER order, const enum CBLAS_TRANSPOSE TransA, const int M, const int N, const double alpha, const double *A, const int lda, const double *X, const int incX, const double beta, double *Y, const int incY)
Definition: cblas_dgemv.c:11
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508

◆ 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 937 of file h1.c.

942  {
944 
945  double *edgeN01 = NULL, *edgeN12 = NULL, *edgeN23 = NULL, *edgeN30 = NULL;
946  if (edgeN != NULL) {
947  edgeN01 = edgeN[0];
948  edgeN12 = edgeN[1];
949  edgeN23 = edgeN[2];
950  edgeN30 = edgeN[3];
951  }
952  double *diff_edgeN01 = NULL, *diff_edgeN12 = NULL, *diff_edgeN23 = NULL,
953  *diff_edgeN30 = NULL;
954  if (diff_edgeN != NULL) {
955  diff_edgeN01 = diff_edgeN[0];
956  diff_edgeN12 = diff_edgeN[1];
957  diff_edgeN23 = diff_edgeN[2];
958  diff_edgeN30 = diff_edgeN[3];
959  }
960  int P[4];
961  int ee = 0;
962  for (; ee < 4; ee++)
963  P[ee] = NBEDGE_H1(p[ee]);
964 
965  int n0 = 0;
966  int n1 = 1;
967  int n2 = 2;
968  int n3 = 3;
969 
970  int ii = 0;
971  for (; ii < GDIM; ii++) {
972  int node_shift = ii * 4;
973  int node_diff_shift = 2 * node_shift;
974 
975  double shape0 = N[node_shift + n0];
976  double shape1 = N[node_shift + n1];
977  double shape2 = N[node_shift + n2];
978  double shape3 = N[node_shift + n3];
979 
980  double ksi01 = (shape1 + shape2 - shape0 - shape3) * sense[n0];
981  double ksi12 = (shape2 + shape3 - shape1 - shape0) * sense[n1];
982  double ksi23 = (shape3 + shape0 - shape2 - shape1) * sense[n2];
983  double ksi30 = (shape0 + shape1 - shape3 - shape2) * sense[n3];
984 
985  double extrude_zeta01 = shape0 + shape1;
986  double extrude_ksi12 = shape1 + shape2;
987  double extrude_zeta23 = shape2 + shape3;
988  double extrude_ksi30 = shape0 + shape3;
989 
990  double bubble_ksi = extrude_ksi12 * extrude_ksi30;
991  double bubble_zeta = extrude_zeta01 * extrude_zeta23;
992 
993  double diff_ksi01[2], diff_ksi12[2], diff_ksi23[2], diff_ksi30[2];
994  double diff_extrude_zeta01[2];
995  double diff_extrude_ksi12[2];
996  double diff_extrude_zeta23[2];
997  double diff_extrude_ksi30[2];
998  double diff_bubble_ksi[2];
999  double diff_bubble_zeta[2];
1000 
1001  int d = 0;
1002  for (; d < 2; d++) {
1003  double diff_shape0 = diffN[node_diff_shift + 2 * n0 + d];
1004  double diff_shape1 = diffN[node_diff_shift + 2 * n1 + d];
1005  double diff_shape2 = diffN[node_diff_shift + 2 * n2 + d];
1006  double diff_shape3 = diffN[node_diff_shift + 2 * n3 + d];
1007  diff_ksi01[d] =
1008  (diff_shape1 + diff_shape2 - diff_shape0 - diff_shape3) * sense[n0];
1009  diff_ksi12[d] =
1010  (diff_shape2 + diff_shape3 - diff_shape1 - diff_shape0) * sense[n1];
1011  diff_ksi23[d] =
1012  (diff_shape3 + diff_shape0 - diff_shape2 - diff_shape1) * sense[n2];
1013  diff_ksi30[d] =
1014  (diff_shape0 + diff_shape1 - diff_shape3 - diff_shape2) * sense[n3];
1015  diff_extrude_zeta01[d] = diff_shape0 + diff_shape1;
1016  diff_extrude_ksi12[d] = diff_shape1 + diff_shape2;
1017  diff_extrude_zeta23[d] = diff_shape2 + diff_shape3;
1018  diff_extrude_ksi30[d] = diff_shape0 + diff_shape3;
1019  diff_bubble_ksi[d] = diff_extrude_ksi12[d] * extrude_ksi30 +
1020  extrude_ksi12 * diff_extrude_ksi30[d];
1021  diff_bubble_zeta[d] = diff_extrude_zeta01[d] * extrude_zeta23 +
1022  extrude_zeta01 * diff_extrude_zeta23[d];
1023  }
1024 
1025  double L01[p[0] + 1], L12[p[1] + 1], L23[p[2] + 1], L30[p[3] + 1];
1026  double diffL01[2 * (p[0] + 1)], diffL12[2 * (p[1] + 1)],
1027  diffL23[2 * (p[2] + 1)], diffL30[2 * (p[3] + 1)];
1028  ierr = base_polynomials(p[0], ksi01, diff_ksi01, L01, diffL01, 2);
1029  CHKERRQ(ierr);
1030  ierr = base_polynomials(p[1], ksi12, diff_ksi12, L12, diffL12, 2);
1031  CHKERRQ(ierr);
1032  ierr = base_polynomials(p[2], ksi23, diff_ksi23, L23, diffL23, 2);
1033  CHKERRQ(ierr);
1034  ierr = base_polynomials(p[3], ksi30, diff_ksi30, L30, diffL30, 2);
1035  CHKERRQ(ierr);
1036 
1037  int shift;
1038  if (edgeN != NULL) {
1039  // edge01
1040  shift = ii * (P[0]);
1041  cblas_dcopy(P[0], L01, 1, &edgeN01[shift], 1);
1042  cblas_dscal(P[0], bubble_ksi * extrude_zeta01, &edgeN01[shift], 1);
1043  // edge12
1044  shift = ii * (P[1]);
1045  cblas_dcopy(P[1], L12, 1, &edgeN12[shift], 1);
1046  cblas_dscal(P[1], bubble_zeta * extrude_ksi12, &edgeN12[shift], 1);
1047  // edge23
1048  shift = ii * (P[2]);
1049  cblas_dcopy(P[2], L23, 1, &edgeN23[shift], 1);
1050  cblas_dscal(P[2], bubble_ksi * extrude_zeta23, &edgeN23[shift], 1);
1051  // edge30
1052  shift = ii * (P[3]);
1053  cblas_dcopy(P[3], L30, 1, &edgeN30[shift], 1);
1054  cblas_dscal(P[3], bubble_zeta * extrude_ksi30, &edgeN30[shift], 1);
1055  }
1056  if (diff_edgeN != NULL) {
1057  if (P[0] > 0) {
1058  // edge01
1059  shift = ii * (P[0]);
1060  bzero(&diff_edgeN01[2 * shift], sizeof(double) * 2 * (P[0]));
1061  int d = 0;
1062  for (; d != 2; ++d) {
1063  cblas_daxpy(P[0], bubble_ksi * extrude_zeta01,
1064  &diffL01[d * (p[0] + 1)], 1, &diff_edgeN01[2 * shift + d],
1065  2);
1066  cblas_daxpy(P[0],
1067  diff_bubble_ksi[d] * extrude_zeta01 +
1068  bubble_ksi * diff_extrude_zeta01[d],
1069  L01, 1, &diff_edgeN01[2 * shift + d], 2);
1070  }
1071  }
1072  if (P[1] > 0) {
1073  // edge12
1074  shift = ii * (P[1]);
1075  bzero(&diff_edgeN12[2 * shift], sizeof(double) * 2 * (P[1]));
1076  int d = 0;
1077  for (; d != 2; ++d) {
1078  cblas_daxpy(P[1], bubble_zeta * extrude_ksi12,
1079  &diffL12[d * (p[1] + 1)], 1, &diff_edgeN12[2 * shift + d],
1080  2);
1081  cblas_daxpy(P[1],
1082  diff_bubble_zeta[d] * extrude_ksi12 +
1083  bubble_zeta * diff_extrude_ksi12[d],
1084  L12, 1, &diff_edgeN12[2 * shift + d], 2);
1085  }
1086  }
1087  if (P[2] > 0) {
1088  // edge23
1089  shift = ii * (P[2]);
1090  bzero(&diff_edgeN23[2 * shift], sizeof(double) * 2 * (P[2]));
1091  int d = 0;
1092  for (; d != 2; ++d) {
1093  cblas_daxpy(P[2], bubble_ksi * extrude_zeta23,
1094  &diffL23[d * (p[2] + 1)], 1, &diff_edgeN23[2 * shift + d],
1095  2);
1096  cblas_daxpy(P[2],
1097  diff_bubble_ksi[d] * extrude_zeta23 +
1098  bubble_ksi * diff_extrude_zeta23[d],
1099  L23, 1, &diff_edgeN23[2 * shift + d], 2);
1100  }
1101  }
1102  if (P[3] > 0) {
1103  // edge30
1104  shift = ii * (P[3]);
1105  bzero(&diff_edgeN30[2 * shift], sizeof(double) * 2 * (P[3]));
1106  int d = 0;
1107  for (; d != 2; ++d) {
1108  cblas_daxpy(P[3], bubble_zeta * extrude_ksi30,
1109  &diffL30[d * (p[3] + 1)], 1, &diff_edgeN30[2 * shift + d],
1110  2);
1111  cblas_daxpy(P[3],
1112  diff_bubble_zeta[d] * extrude_ksi30 +
1113  bubble_zeta * diff_extrude_ksi30[d],
1114  L30, 1, &diff_edgeN30[2 * shift + d], 2);
1115  }
1116  }
1117  }
1118  }
1120 }
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
void cblas_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_daxpy.c:11
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
void cblas_dcopy(const int N, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_dcopy.c:11
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
void cblas_dscal(const int N, const double alpha, double *X, const int incX)
Definition: cblas_dscal.c:11
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
CHKERRQ(ierr)
static PetscErrorCode ierr
Definition: h1.c:25
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(*)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)  base_polynomials 
)

Definition at line 245 of file h1.c.

250  {
252 
253  int P[6];
254  int ee = 0;
255  for (; ee < 6; ee++)
256  P[ee] = NBEDGE_H1(p[ee]);
257  int edges_nodes[2 * 6] = {0, 1, 1, 2, 2, 0, 0, 3, 1, 3, 2, 3};
258  double diff_ksi01[3], diff_ksi12[3], diff_ksi20[3];
259  double diff_ksi03[3], diff_ksi13[3], diff_ksi23[3];
260  double *edges_diff_ksi[6] = {diff_ksi01, diff_ksi12, diff_ksi20,
261  diff_ksi03, diff_ksi13, diff_ksi23};
262  for (ee = 0; ee < 6; ee++) {
263  int dd = 0;
264  for (; dd < 3; dd++) {
265  edges_diff_ksi[ee][dd] = (diffN[edges_nodes[2 * ee + 1] * 3 + dd] -
266  diffN[edges_nodes[2 * ee + 0] * 3 + dd]) *
267  sense[ee];
268  }
269  }
270  int ii = 0;
271  for (; ii < GDIM; ii++) {
272  int node_shift = ii * 4;
273  double edges_ksi[6];
274  for (ee = 0; ee < 6; ee++) {
275  edges_ksi[ee] = (N[node_shift + edges_nodes[2 * ee + 1]] -
276  N[node_shift + edges_nodes[2 * ee + 0]]) *
277  sense[ee];
278  }
279  for (ee = 0; ee < 6; ee++) {
280  if (P[ee] == 0)
281  continue;
282  double L[p[ee] + 1], diffL[3 * (p[ee] + 1)];
283  ierr = base_polynomials(p[ee], edges_ksi[ee], edges_diff_ksi[ee], L,
284  diffL, 3);
285  CHKERRQ(ierr);
286  double v = N[node_shift + edges_nodes[2 * ee + 0]] *
287  N[node_shift + edges_nodes[2 * ee + 1]];
288  if (edgeN != NULL)
289  if (edgeN[ee] != NULL) {
290  int shift = ii * P[ee];
291  cblas_dcopy(P[ee], L, 1, &edgeN[ee][shift], 1);
292  cblas_dscal(P[ee],
293  N[node_shift + edges_nodes[2 * ee + 0]] *
294  N[node_shift + edges_nodes[2 * ee + 1]],
295  &edgeN[ee][shift], 1);
296  }
297  if (diff_edgeN != NULL)
298  if (diff_edgeN[ee] != NULL) {
299  int shift = ii * P[ee];
300  bzero(&diff_edgeN[ee][3 * shift], sizeof(double) * 3 * P[ee]);
301  int dd = 0;
302  for (; dd < 3; dd++) {
303  cblas_daxpy(P[ee], v, &diffL[dd * (p[ee] + 1)], 1,
304  &diff_edgeN[ee][3 * shift + dd], 3);
305  cblas_daxpy(P[ee],
306  diffN[3 * edges_nodes[2 * ee + 0] + dd] *
307  N[node_shift + edges_nodes[2 * ee + 1]] +
308  N[node_shift + edges_nodes[2 * ee + 0]] *
309  diffN[3 * edges_nodes[2 * ee + 1] + dd],
310  L, 1, &diff_edgeN[ee][3 * shift + dd], 3);
311  }
312  }
313  }
314  }
316 }
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
void cblas_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_daxpy.c:11
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
void cblas_dcopy(const int N, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_dcopy.c:11
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
void cblas_dscal(const int N, const double alpha, double *X, const int incX)
Definition: cblas_dscal.c:11
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
CHKERRQ(ierr)
static PetscErrorCode ierr
Definition: h1.c:25
const int N
Definition: speed_test.cpp:3

◆ 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 27 of file h1.c.

32  {
34 
35  double *edgeN01 = NULL, *edgeN12 = NULL, *edgeN20 = NULL;
36  if (edgeN != NULL) {
37  edgeN01 = edgeN[0];
38  edgeN12 = edgeN[1];
39  edgeN20 = edgeN[2];
40  }
41  double *diff_edgeN01 = NULL, *diff_edgeN12 = NULL, *diff_edgeN20 = NULL;
42  if (diff_edgeN != NULL) {
43  diff_edgeN01 = diff_edgeN[0];
44  diff_edgeN12 = diff_edgeN[1];
45  diff_edgeN20 = diff_edgeN[2];
46  }
47  int P[3];
48  int ee = 0;
49  for (; ee < 3; ee++) {
50  P[ee] = NBEDGE_H1(p[ee]);
51  }
52  int dd = 0;
53  double diff_ksi01[2], diff_ksi12[2], diff_ksi20[2];
54  if (diff_edgeN != NULL) {
55  for (; dd < 2; dd++) {
56  diff_ksi01[dd] = (diffN[1 * 2 + dd] - diffN[0 * 2 + dd]) * sense[0];
57  diff_ksi12[dd] = (diffN[2 * 2 + dd] - diffN[1 * 2 + dd]) * sense[1];
58  diff_ksi20[dd] = (diffN[0 * 2 + dd] - diffN[2 * 2 + dd]) * sense[2];
59  }
60  }
61  int ii = 0;
62  for (; ii < GDIM; ii++) {
63  int node_shift = ii * 3;
64  double ksi01 = (N[node_shift + 1] - N[node_shift + 0]) * sense[0];
65  double ksi12 = (N[node_shift + 2] - N[node_shift + 1]) * sense[1];
66  double ksi20 = (N[node_shift + 0] - N[node_shift + 2]) * sense[2];
67  double L01[p[0] + 1], L12[p[1] + 1], L20[p[2] + 1];
68  double diffL01[2 * (p[0] + 1)], diffL12[2 * (p[1] + 1)],
69  diffL20[2 * (p[2] + 1)];
70  ierr = base_polynomials(p[0], ksi01, diff_ksi01, L01, diffL01, 2);
71  CHKERRQ(ierr);
72  ierr = base_polynomials(p[1], ksi12, diff_ksi12, L12, diffL12, 2);
73  CHKERRQ(ierr);
74  ierr = base_polynomials(p[2], ksi20, diff_ksi20, L20, diffL20, 2);
75  CHKERRQ(ierr);
76  int shift;
77  if (edgeN != NULL) {
78  // edge01
79  shift = ii * (P[0]);
80  cblas_dcopy(P[0], L01, 1, &edgeN01[shift], 1);
81  cblas_dscal(P[0], N[node_shift + 0] * N[node_shift + 1], &edgeN01[shift],
82  1);
83  // edge12
84  shift = ii * (P[1]);
85  cblas_dcopy(P[1], L12, 1, &edgeN12[shift], 1);
86  cblas_dscal(P[1], N[node_shift + 1] * N[node_shift + 2], &edgeN12[shift],
87  1);
88  // edge20
89  shift = ii * (P[2]);
90  cblas_dcopy(P[2], L20, 1, &edgeN20[shift], 1);
91  cblas_dscal(P[2], N[node_shift + 0] * N[node_shift + 2], &edgeN20[shift],
92  1);
93  }
94  if (diff_edgeN != NULL) {
95  if (P[0] > 0) {
96  // edge01
97  shift = ii * (P[0]);
98  bzero(&diff_edgeN01[2 * shift], sizeof(double) * 2 * (P[0]));
99  // diffX
100  cblas_daxpy(P[0], N[node_shift + 0] * N[node_shift + 1],
101  &diffL01[0 * (p[0] + 1)], 1, &diff_edgeN01[2 * shift + 0],
102  2);
103  cblas_daxpy(P[0],
104  diffN[2 * 0 + 0] * N[node_shift + 1] +
105  N[node_shift + 0] * diffN[2 * 1 + 0],
106  L01, 1, &diff_edgeN01[2 * shift + 0], 2);
107  // diff Y
108  cblas_daxpy(P[0], N[node_shift + 0] * N[node_shift + 1],
109  &diffL01[1 * (p[0] + 1)], 1, &diff_edgeN01[2 * shift + 1],
110  2);
111  cblas_daxpy(P[0],
112  diffN[2 * 0 + 1] * N[node_shift + 1] +
113  N[node_shift + 0] * diffN[2 * 1 + 1],
114  L01, 1, &diff_edgeN01[2 * shift + 1], 2);
115  }
116  if (P[1] > 0) {
117  // edge12
118  shift = ii * (P[1]);
119  bzero(&diff_edgeN12[2 * shift], sizeof(double) * 2 * (P[1]));
120  // diffX
121  cblas_daxpy(P[1], N[node_shift + 1] * N[node_shift + 2],
122  &diffL12[0 * (p[1] + 1)], 1, &diff_edgeN12[2 * shift + 0],
123  2);
124  cblas_daxpy(P[1],
125  diffN[2 * 1 + 0] * N[node_shift + 2] +
126  N[node_shift + 1] * diffN[2 * 2 + 0],
127  L12, 1, &diff_edgeN12[2 * shift + 0], 2);
128  // diffY
129  cblas_daxpy(P[1], N[node_shift + 1] * N[node_shift + 2],
130  &diffL12[1 * (p[1] + 1)], 1, &diff_edgeN12[2 * shift + 1],
131  2);
132  cblas_daxpy(P[1],
133  diffN[2 * 1 + 1] * N[node_shift + 2] +
134  N[node_shift + 1] * diffN[2 * 2 + 1],
135  L12, 1, &diff_edgeN12[2 * shift + 1], 2);
136  }
137  if (P[2] > 0) {
138  // edge20
139  shift = ii * (P[2]);
140  bzero(&diff_edgeN20[2 * shift], sizeof(double) * 2 * (P[2]));
141  // diffX
142  cblas_daxpy(P[2], N[node_shift + 0] * N[node_shift + 2],
143  &diffL20[0 * (p[2] + 1)], 1, &diff_edgeN20[2 * shift + 0],
144  2);
145  cblas_daxpy(P[2],
146  diffN[2 * 0 + 0] * N[node_shift + 2] +
147  N[node_shift + 0] * diffN[2 * 2 + 0],
148  L20, 1, &diff_edgeN20[2 * shift + 0], 2);
149  // diffY
150  cblas_daxpy(P[2], N[node_shift + 0] * N[node_shift + 2],
151  &diffL20[1 * (p[2] + 1)], 1, &diff_edgeN20[2 * shift + 1],
152  2);
153  cblas_daxpy(P[2],
154  diffN[2 * 0 + 1] * N[node_shift + 2] +
155  N[node_shift + 0] * diffN[2 * 2 + 1],
156  L20, 1, &diff_edgeN20[2 * shift + 1], 2);
157  }
158  }
159  }
161 }
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
void cblas_daxpy(const int N, const double alpha, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_daxpy.c:11
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
void cblas_dcopy(const int N, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_dcopy.c:11
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
void cblas_dscal(const int N, const double alpha, double *X, const int incX)
Definition: cblas_dscal.c:11
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
CHKERRQ(ierr)
static PetscErrorCode ierr
Definition: h1.c:25
const int N
Definition: speed_test.cpp:3

◆ H1_FaceGradientOfDeformation_hierachical()

PetscErrorCode H1_FaceGradientOfDeformation_hierachical ( int  p,
double diffN,
double dofs,
double F 
)
Deprecated:
use H1_FaceGradientOfDeformation_hierarchical

Definition at line 1127 of file h1.c.

1129  {
1130  return H1_FaceGradientOfDeformation_hierarchical(p, diffN, dofs, F);
1131 }
PetscErrorCode H1_FaceGradientOfDeformation_hierarchical(int p, double *diffN, double *dofs, double *F)
Definition: h1.c:564

◆ H1_FaceGradientOfDeformation_hierarchical()

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

Definition at line 564 of file h1.c.

566  {
568  int col, row = 0;
569  for (; row < 3; row++)
570  for (col = 0; col < 3; col++)
571  F[3 * row + col] =
572  cblas_ddot(NBFACETRI_H1(p), &diffN[col], 3, &dofs[row], 3);
574 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: cblas_ddot.c:12
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.

◆ H1_FaceShapeDiffMBTETinvJ()

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

Definition at line 518 of file h1.c.

520  {
522  int ff = 0, ii, gg;
523  for (; ff < 4; ff++) {
524  for (ii = 0; ii < NBFACETRI_H1(p[ff]); ii++) {
525  for (gg = 0; gg < GDIM; gg++) {
526  int shift1 = NBFACETRI_H1(base_p[ff]) * gg;
527  int shift2 = NBFACETRI_H1(p[ff]) * gg;
528  cblas_dgemv(CblasRowMajor, CblasTrans, 3, 3, 1., invJac, 3,
529  &(face_diffN[ff])[3 * shift1 + 3 * ii], 1, 0.,
530  &(face_diffNinvJac[ff])[3 * shift2 + 3 * ii], 1);
531  }
532  }
533  }
535 }
void cblas_dgemv(const enum CBLAS_ORDER order, const enum CBLAS_TRANSPOSE TransA, const int M, const int N, const double alpha, const double *A, const int lda, const double *X, const int incX, const double beta, double *Y, const int incY)
Definition: cblas_dgemv.c:11
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.

◆ 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 317 of file h1.c.

322  {
324 
325  int P[4];
326  int ff = 0;
327  for (; ff < 4; ff++)
328  P[ff] = NBFACETRI_H1(p[ff]);
329  double diff_ksiL0F0[3], diff_ksiL1F0[3];
330  double diff_ksiL0F1[3], diff_ksiL1F1[3];
331  double diff_ksiL0F2[3], diff_ksiL1F2[3];
332  double diff_ksiL0F3[3], diff_ksiL1F3[3];
333  double *diff_ksi_faces[] = {diff_ksiL0F0, diff_ksiL1F0, diff_ksiL0F1,
334  diff_ksiL1F1, diff_ksiL0F2, diff_ksiL1F2,
335  diff_ksiL0F3, diff_ksiL1F3};
336  int dd = 0;
337  for (; dd < 3; dd++) {
338  for (ff = 0; ff < 4; ff++) {
339  diff_ksi_faces[ff * 2 + 0][dd] =
340  (diffN[faces_nodes[3 * ff + 1] * 3 + dd] -
341  diffN[faces_nodes[3 * ff + 0] * 3 + dd]);
342  diff_ksi_faces[ff * 2 + 1][dd] =
343  (diffN[faces_nodes[3 * ff + 2] * 3 + dd] -
344  diffN[faces_nodes[3 * ff + 0] * 3 + dd]);
345  }
346  }
347  int ii = 0;
348  for (; ii < GDIM; ii++) {
349  int node_shift = ii * 4;
350  double ksi_faces[8];
351  for (ff = 0; ff < 4; ff++) {
352  ksi_faces[2 * ff + 0] = N[node_shift + faces_nodes[3 * ff + 1]] -
353  N[node_shift + faces_nodes[3 * ff + 0]];
354  ksi_faces[2 * ff + 1] = N[node_shift + faces_nodes[3 * ff + 2]] -
355  N[node_shift + faces_nodes[3 * ff + 0]];
356  }
357  int shift;
358  for (ff = 0; ff < 4; ff++) {
359  if (P[ff] == 0)
360  continue;
361  double L0[p[ff] + 1], L1[p[ff] + 1];
362  double diffL0[3 * (p[ff] + 1)], diffL1[3 * (p[ff] + 1)];
363  ierr = base_polynomials(p[ff], ksi_faces[ff * 2 + 0],
364  diff_ksi_faces[ff * 2 + 0], L0, diffL0, 3);
365  CHKERRQ(ierr);
366  ierr = base_polynomials(p[ff], ksi_faces[ff * 2 + 1],
367  diff_ksi_faces[ff * 2 + 1], L1, diffL1, 3);
368  CHKERRQ(ierr);
369  double v = N[node_shift + faces_nodes[3 * ff + 0]] *
370  N[node_shift + faces_nodes[3 * ff + 1]] *
371  N[node_shift + faces_nodes[3 * ff + 2]];
372  double v2[3] = {0, 0, 0};
373  dd = 0;
374  double n1n2 = N[node_shift + faces_nodes[3 * ff + 1]] *
375  N[node_shift + faces_nodes[3 * ff + 2]];
376  double n0n2 = N[node_shift + faces_nodes[3 * ff + 0]] *
377  N[node_shift + faces_nodes[3 * ff + 2]];
378  double n0n1 = N[node_shift + faces_nodes[3 * ff + 0]] *
379  N[node_shift + faces_nodes[3 * ff + 1]];
380  for (; dd < 3; dd++) {
381  v2[dd] = diffN[3 * faces_nodes[3 * ff + 0] + dd] * n1n2 +
382  diffN[3 * faces_nodes[3 * ff + 1] + dd] * n0n2 +
383  diffN[3 * faces_nodes[3 * ff + 2] + dd] * n0n1;
384  }
385  shift = ii * P[ff];
386  int jj = 0;
387  int oo = 0;
388  for (; oo <= (p[ff] - 3); oo++) {
389  int pp0 = 0;
390  for (; pp0 <= oo; pp0++) {
391  int pp1 = oo - pp0;
392  if (pp1 >= 0) {
393  if (faceN != NULL)
394  if (faceN[ff] != NULL) {
395  faceN[ff][shift + jj] = v * L0[pp0] * L1[pp1];
396  }
397  if (diff_faceN != NULL)
398  if (diff_faceN[ff] != NULL) {
399  dd = 0;
400  double L0L1 = L0[pp0] * L1[pp1];
401  for (; dd < 3; dd++) {
402  diff_faceN[ff][3 * shift + 3 * jj + dd] =
403  (L0[pp0] * diffL1[dd * (p[ff] + 1) + pp1] +
404  diffL0[dd * (p[ff] + 1) + pp0] * L1[pp1]) *
405  v;
406  diff_faceN[ff][3 * shift + 3 * jj + dd] += L0L1 * v2[dd];
407  }
408  }
409  jj++;
410  }
411  }
412  }
413  if (jj != P[ff])
414  SETERRQ2(PETSC_COMM_SELF, 1, "wrong order %d != %d", jj, P[ff]);
415  }
416  }
418 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
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
CHKERRQ(ierr)
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
static PetscErrorCode ierr
Definition: h1.c:25
const int N
Definition: speed_test.cpp:3

◆ 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 162 of file h1.c.

167  {
169 
170  int P = NBFACETRI_H1(p);
171  if (P == 0)
173  double diff_ksiL0[2], diff_ksiL1[2];
174  double *diff_ksi_faces[] = {diff_ksiL0, diff_ksiL1};
175  int dd = 0;
176  if (diff_faceN != NULL) {
177  for (; dd < 2; dd++) {
178  diff_ksi_faces[0][dd] =
179  diffN[face_nodes[1] * 2 + dd] - diffN[face_nodes[0] * 2 + dd];
180  diff_ksi_faces[1][dd] =
181  diffN[face_nodes[2] * 2 + dd] - diffN[face_nodes[0] * 2 + dd];
182  }
183  }
184  int ii = 0;
185  for (; ii < GDIM; ii++) {
186  int node_shift = ii * 3;
187  double ksi_faces[2];
188  ksi_faces[0] =
189  N[node_shift + face_nodes[1]] - N[node_shift + face_nodes[0]];
190  ksi_faces[1] =
191  N[node_shift + face_nodes[2]] - N[node_shift + face_nodes[0]];
192  double L0[p + 1], L1[p + 1];
193  double diffL0[2 * (p + 1)], diffL1[2 * (p + 1)];
194  ierr = base_polynomials(p, ksi_faces[0], diff_ksi_faces[0], L0, diffL0, 2);
195  CHKERRQ(ierr);
196  ierr = base_polynomials(p, ksi_faces[1], diff_ksi_faces[1], L1, diffL1, 2);
197  CHKERRQ(ierr);
198  double v = N[node_shift + face_nodes[0]] * N[node_shift + face_nodes[1]] *
199  N[node_shift + face_nodes[2]];
200  double v2[2] = {0, 0};
201  if (diff_faceN != NULL) {
202  double n1n2 =
203  N[node_shift + face_nodes[1]] * N[node_shift + face_nodes[2]];
204  double n0n2 =
205  N[node_shift + face_nodes[0]] * N[node_shift + face_nodes[2]];
206  double n0n1 =
207  N[node_shift + face_nodes[0]] * N[node_shift + face_nodes[1]];
208  dd = 0;
209  for (; dd < 2; dd++) {
210  v2[dd] = diffN[face_nodes[0] * 2 + dd] * n1n2 +
211  diffN[face_nodes[1] * 2 + dd] * n0n2 +
212  diffN[face_nodes[2] * 2 + dd] * n0n1;
213  }
214  }
215  int shift = ii * P;
216  int jj = 0;
217  int oo = 0;
218  for (; oo <= (p - 3); oo++) {
219  int pp0 = 0;
220  for (; pp0 <= oo; pp0++) {
221  int pp1 = oo - pp0;
222  if (pp1 >= 0) {
223  if (faceN != NULL) {
224  faceN[shift + jj] = v * L0[pp0] * L1[pp1];
225  }
226  if (diff_faceN != NULL) {
227  dd = 0;
228  for (; dd < 2; dd++) {
229  double *val = &diff_faceN[2 * shift + 2 * jj + dd];
230  *val = (L0[pp0] * diffL1[dd * (p + 1) + pp1] +
231  diffL0[dd * (p + 1) + pp0] * L1[pp1]) *
232  v;
233  *val += L0[pp0] * L1[pp1] * v2[dd];
234  }
235  }
236  jj++;
237  }
238  }
239  }
240  if (jj != P)
241  SETERRQ2(PETSC_COMM_SELF, 1, "wrong order %d != %d", jj, P);
242  }
244 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
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
CHKERRQ(ierr)
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
static PetscErrorCode ierr
Definition: h1.c:25
const int N
Definition: speed_test.cpp:3

◆ 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 586 of file h1.c.

591  {
593  // TODO: return separately components of the tensor product between two edges
594 
595  int P[3];
596  int ff = 0;
597  for (; ff < 3; ff++)
598  P[ff] = NBFACEQUAD_H1(p[ff]);
599 
600  double ksi_faces[6];
601  double diff_ksiL0F0[3], diff_ksiL3F0[3];
602  double diff_ksiL0F1[3], diff_ksiL3F1[3];
603  double diff_ksiL0F2[3], diff_ksiL3F2[3];
604  double *diff_ksi_faces[] = {diff_ksiL0F0, diff_ksiL3F0, diff_ksiL0F1,
605  diff_ksiL3F1, diff_ksiL0F2, diff_ksiL3F2};
606  int ii = 0;
607  for (; ii < GDIM; ii++) {
608  int node_shift = ii * 6;
609  int node_diff_shift = 3 * node_shift;
610  int ff = 0;
611  for (; ff < 3; ff++) {
612  if (P[ff] == 0)
613  continue;
614  int n0 = faces_nodes[4 * ff + 0];
615  int n1 = faces_nodes[4 * ff + 1];
616  int n2 = faces_nodes[4 * ff + 2];
617  int n3 = faces_nodes[4 * ff + 3];
618  int e0 = 2 * ff + 0;
619  int e1 = 2 * ff + 1;
620  double ksi0 = N[node_shift + n0] + N[node_shift + n3];
621  double ksi1 = N[node_shift + n1] + N[node_shift + n2];
622  double eta0 = N[node_shift + n0] + N[node_shift + n1];
623  double eta1 = N[node_shift + n2] + N[node_shift + n3];
624  ksi_faces[e0] = ksi1 - ksi0;
625  ksi_faces[e1] = eta1 - eta0;
626  // double eps = 1e-14;
627  // if (ksi_faces[e0] < -(1 + eps) || ksi_faces[e0] > (1 + eps)) {
628  // SETERRQ(PETSC_COMM_SELF, 1, "Data inconsistency");
629  // }
630  // if (ksi_faces[e1] < -(1 + eps) || ksi_faces[e1] > (1 + eps)) {
631  // SETERRQ(PETSC_COMM_SELF, 1, "Data inconsistency");
632  // }
633  int dd = 0;
634  for (; dd < 3; dd++) {
635  double diff_ksi0 = diffN[node_diff_shift + 3 * n0 + dd] +
636  diffN[node_diff_shift + 3 * n3 + dd];
637  double diff_ksi1 = diffN[node_diff_shift + 3 * n1 + dd] +
638  diffN[node_diff_shift + 3 * n2 + dd];
639  double diff_eta0 = diffN[node_diff_shift + 3 * n0 + dd] +
640  diffN[node_diff_shift + 3 * n1 + dd];
641  double diff_eta1 = diffN[node_diff_shift + 3 * n2 + dd] +
642  diffN[node_diff_shift + 3 * n3 + dd];
643  diff_ksi_faces[e0][dd] = diff_ksi1 - diff_ksi0;
644  diff_ksi_faces[e1][dd] = diff_eta1 - diff_eta0;
645  }
646  double L0[p[ff] + 1], L1[p[ff] + 1];
647  double diffL0[3 * (p[ff] + 1)], diffL1[3 * (p[ff] + 1)];
648  ierr = base_polynomials(p[ff], ksi_faces[e0], diff_ksi_faces[e0], L0,
649  diffL0, 3);
650  CHKERRQ(ierr);
651  ierr = base_polynomials(p[ff], ksi_faces[e1], diff_ksi_faces[e1], L1,
652  diffL1, 3);
653  CHKERRQ(ierr);
654 
655  double v = N[node_shift + n0] * N[node_shift + n2] +
656  N[node_shift + n1] * N[node_shift + n3];
657 
658  double diff_v[3];
659  dd = 0;
660  for (; dd < 3; ++dd)
661  diff_v[dd] =
662 
663  diffN[node_diff_shift + 3 * n0 + dd] * N[node_shift + n2] +
664  N[node_shift + n0] * diffN[node_diff_shift + 3 * n2 + dd] +
665 
666  diffN[node_diff_shift + 3 * n1 + dd] * N[node_shift + n3] +
667  N[node_shift + n1] * diffN[node_diff_shift + 3 * n3 + dd];
668 
669  int shift;
670  shift = ii * P[ff];
671  int jj = 0;
672  int oo = 0;
673  for (; oo <= (p[ff] - 4); oo++) {
674  int pp0 = 0;
675  for (; pp0 <= oo; pp0++) {
676  int pp1 = oo - pp0;
677  if (pp1 >= 0) {
678  if (faceN != NULL) {
679  if (faceN[ff] != NULL) {
680  faceN[ff][shift + jj] = v * L0[pp0] * L1[pp1];
681  }
682  }
683  if (diff_faceN != NULL) {
684  if (diff_faceN[ff] != NULL) {
685  dd = 0;
686  for (; dd < 3; dd++) {
687  diff_faceN[ff][3 * shift + 3 * jj + dd] =
688 
689  (L0[pp0] * diffL1[dd * (p[ff] + 1) + pp1] +
690  diffL0[dd * (p[ff] + 1) + pp0] * L1[pp1]) *
691  v +
692 
693  L0[pp0] * L1[pp1] * diff_v[dd];
694  }
695  }
696  }
697  jj++;
698  }
699  }
700  }
701  if (jj != P[ff])
702  SETERRQ2(PETSC_COMM_SELF, 1, "wrong order %d != %d", jj, P[ff]);
703  }
704  }
706 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
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
CHKERRQ(ierr)
#define NBFACEQUAD_H1(P)
Number of base functions on quad for H1 space.
static PetscErrorCode ierr
Definition: h1.c:25
const int N
Definition: speed_test.cpp:3

◆ 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 827 of file h1.c.

832  {
834  // TODO: return separately components of the tensor product between two edges
835 
836  int P = NBFACEQUAD_H1(p);
837  if (P == 0)
839  double diff_ksiL0F0[2], diff_ksiL1F0[2];
840  double *diff_ksi_faces[] = {diff_ksiL0F0, diff_ksiL1F0};
841  int ii = 0;
842  for (; ii < GDIM; ii++) {
843  int node_shift = ii * 4;
844  int node_diff_shift = 2 * node_shift;
845 
846  int n0 = faces_nodes[0];
847  int n1 = faces_nodes[1];
848  int n2 = faces_nodes[2];
849  int n3 = faces_nodes[3];
850  double ksi0 = N[node_shift + n0] + N[node_shift + n3];
851  double ksi1 = N[node_shift + n1] + N[node_shift + n2];
852  double eta0 = N[node_shift + n0] + N[node_shift + n1];
853  double eta1 = N[node_shift + n2] + N[node_shift + n3];
854  double ksi_faces = ksi1 - ksi0;
855  double eta_faces = eta1 - eta0;
856  // double eps = 1e-14;
857  // if (ksi_faces < -(1 + eps) || ksi_faces > (1 + eps)) {
858  // SETERRQ(PETSC_COMM_SELF, 1, "Data inconsistency");
859  // }
860  // if (eta_faces < -(1 + eps) || eta_faces > (1 + eps)) {
861  // SETERRQ(PETSC_COMM_SELF, 1, "Data inconsistency");
862  // }
863  double diff_ksi_faces[2];
864  double diff_eta_faces[2];
865  int dd = 0;
866  for (; dd < 2; dd++) {
867  double diff_ksi0 = diffN[node_diff_shift + 2 * n0 + dd] +
868  diffN[node_diff_shift + 2 * n3 + dd];
869  double diff_ksi1 = diffN[node_diff_shift + 2 * n1 + dd] +
870  diffN[node_diff_shift + 2 * n2 + dd];
871  double diff_eta0 = diffN[node_diff_shift + 2 * n0 + dd] +
872  diffN[node_diff_shift + 2 * n1 + dd];
873  double diff_eta1 = diffN[node_diff_shift + 2 * n2 + dd] +
874  diffN[node_diff_shift + 2 * n3 + dd];
875 
876  diff_ksi_faces[dd] = diff_ksi1 - diff_ksi0;
877  diff_eta_faces[dd] = diff_eta1 - diff_eta0;
878  }
879  double L0[p + 1], L1[p + 1];
880  double diffL0[2 * (p + 1)], diffL1[2 * (p + 1)];
881  ierr =
882  base_polynomials(p, ksi_faces, diff_ksi_faces, L0, diffL0, 2);
883  CHKERRQ(ierr);
884  ierr =
885  base_polynomials(p, eta_faces, diff_eta_faces, L1, diffL1, 2);
886  CHKERRQ(ierr);
887 
888  double v = N[node_shift + n0] * N[node_shift + n2] +
889  N[node_shift + n1] * N[node_shift + n3];
890 
891  double diff_v[2];
892  dd = 0;
893  for (; dd < 2; ++dd)
894  diff_v[dd] =
895 
896  diffN[node_diff_shift + 2 * n0 + dd] * N[node_shift + n2] +
897  N[node_shift + n0] * diffN[node_diff_shift + 2 * n2 + dd] +
898 
899  diffN[node_diff_shift + 2 * n1 + dd] * N[node_shift + n3] +
900  N[node_shift + n1] * diffN[node_diff_shift + 2 * n3 + dd];
901 
902 
903  int shift;
904  shift = ii * P;
905  int jj = 0;
906  int oo = 0;
907  for (; oo <= (p - 4); oo++) {
908  int pp0 = 0;
909  for (; pp0 <= oo; pp0++) {
910  int pp1 = oo - pp0;
911  if (pp1 >= 0) {
912  if (faceN != NULL) {
913  faceN[shift + jj] = L0[pp0] * L1[pp1] * v;
914  }
915  if (diff_faceN != NULL) {
916  dd = 0;
917  for (; dd < 2; dd++) {
918  diff_faceN[2 * shift + 2 * jj + dd] =
919  (L0[pp0] * diffL1[dd * (p + 1) + pp1] +
920  diffL0[dd * (p + 1) + pp0] * L1[pp1]) *
921  v;
922  diff_faceN[2 * shift + 2 * jj + dd] +=
923  L0[pp0] * L1[pp1] * diff_v[dd];
924  }
925  }
926  jj++;
927  }
928  }
929  }
930  if (jj != P)
931  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
932  "wrong order %d != %d", jj, P);
933  }
935 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
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
CHKERRQ(ierr)
#define NBFACEQUAD_H1(P)
Number of base functions on quad for H1 space.
static PetscErrorCode ierr
Definition: h1.c:25
const int N
Definition: speed_test.cpp:3

◆ H1_VolumeGradientOfDeformation_hierachical()

PetscErrorCode H1_VolumeGradientOfDeformation_hierachical ( int  p,
double diffN,
double dofs,
double F 
)
Deprecated:
use H1_VolumeGradientOfDeformation_hierarchical

Definition at line 1132 of file h1.c.

1134  {
1135  return H1_VolumeGradientOfDeformation_hierarchical(p, diffN, dofs, F);
1136 }
PetscErrorCode H1_VolumeGradientOfDeformation_hierarchical(int p, double *diffN, double *dofs, double *F)
Definition: h1.c:575

◆ H1_VolumeGradientOfDeformation_hierarchical()

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

Definition at line 575 of file h1.c.

577  {
579  int col, row = 0;
580  for (; row < 3; row++)
581  for (col = 0; col < 3; col++)
582  F[3 * row + col] =
583  cblas_ddot(NBVOLUMETET_H1(p), &diffN[col], 3, &dofs[row], 3);
585 }
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: cblas_ddot.c:12

◆ H1_VolumeShapeDiffMBTETinvJ()

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

Definition at line 536 of file h1.c.

539  {
541  int ii, gg;
542  for (ii = 0; ii < NBVOLUMETET_H1(p); ii++) {
543  for (gg = 0; gg < GDIM; gg++) {
544  int shift1 = NBVOLUMETET_H1(base_p) * gg;
545  int shift2 = NBVOLUMETET_H1(p) * gg;
546  cblas_dgemv(CblasRowMajor, CblasTrans, 3, 3, 1., invJac, 3,
547  &(volume_diffN)[3 * shift1 + 3 * ii], 1, 0.,
548  &(volume_diffNinvJac)[3 * shift2 + 3 * ii], 1);
549  }
550  }
552 }
void cblas_dgemv(const enum CBLAS_ORDER order, const enum CBLAS_TRANSPOSE TransA, const int M, const int N, const double alpha, const double *A, const int lda, const double *X, const int incX, const double beta, double *Y, const int incY)
Definition: cblas_dgemv.c:11
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508

◆ 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 707 of file h1.c.

712  {
714 
715  int P = NBVOLUMEPRISM_H1(p);
716  if (P == 0)
718  double diff_ksiL0[3], diff_ksiL1[3], diff_ksiL2[3];
719  int ii = 0;
720  for (; ii < GDIM; ii++) {
721  int node_shift = ii * 6;
722  int node_diff_shift = ii * 18;
723 
724  double n0 = N[node_shift + 0];
725  double n1 = N[node_shift + 1];
726  double n2 = N[node_shift + 2];
727  double n3 = N[node_shift + 3];
728  double n4 = N[node_shift + 4];
729  double n5 = N[node_shift + 5];
730 
731  double ksiL0 = n1 + n4 - n0 - n3;
732  double ksiL1 = n2 + n5 - n0 - n3;
733  double ksiL2 = (n3 + n4 + n5) - (n0 + n1 + n2);
734 
735  int dd = 0;
736  for (; dd < 3; dd++) {
737  double diff_n0 = diffN[node_diff_shift + 3 * 0 + dd];
738  double diff_n1 = diffN[node_diff_shift + 3 * 1 + dd];
739  double diff_n2 = diffN[node_diff_shift + 3 * 2 + dd];
740  double diff_n3 = diffN[node_diff_shift + 3 * 3 + dd];
741  double diff_n4 = diffN[node_diff_shift + 3 * 4 + dd];
742  double diff_n5 = diffN[node_diff_shift + 3 * 5 + dd];
743  diff_ksiL0[dd] = (diff_n1 + diff_n4) - (diff_n0 + diff_n3);
744  diff_ksiL1[dd] = (diff_n2 + diff_n5) - (diff_n0 + diff_n3);
745  diff_ksiL2[dd] =
746  (diff_n3 + diff_n4 + diff_n5) - (diff_n0 + diff_n1 + diff_n2);
747  }
748 
749  double L0[p + 1], L1[p + 1], L2[p + 1];
750  double diffL0[3 * (p + 1)], diffL1[3 * (p + 1)], diffL2[3 * (p + 1)];
751  ierr = base_polynomials(p, ksiL0, diff_ksiL0, L0, diffL0, 3);
752  CHKERRQ(ierr);
753  ierr = base_polynomials(p, ksiL1, diff_ksiL1, L1, diffL1, 3);
754  CHKERRQ(ierr);
755  ierr = base_polynomials(p, ksiL2, diff_ksiL2, L2, diffL2, 3);
756  CHKERRQ(ierr);
757 
758  double v_tri = (n0 + n3) * (n1 + n4) * (n2 + n5);
759  double v_edge = (n0 + n1 + n2) * (n3 + n4 + n5);
760  double v = v_tri * v_edge;
761 
762  double diff_v_tri[3];
763  double diff_v_edge[3];
764  double diff_v[3];
765  dd = 0;
766  for (; dd < 3; dd++) {
767  double diff_n0 = diffN[node_diff_shift + 3 * 0 + dd];
768  double diff_n1 = diffN[node_diff_shift + 3 * 1 + dd];
769  double diff_n2 = diffN[node_diff_shift + 3 * 2 + dd];
770  double diff_n3 = diffN[node_diff_shift + 3 * 3 + dd];
771  double diff_n4 = diffN[node_diff_shift + 3 * 4 + dd];
772  double diff_n5 = diffN[node_diff_shift + 3 * 5 + dd];
773  diff_v_tri[dd] = (diff_n0 + diff_n3) * (n1 + n4) * (n2 + n5) +
774  (diff_n1 + diff_n4) * (n0 + n3) * (n2 + n5) +
775  (diff_n2 + diff_n5) * (n0 + n3) * (n1 + n4);
776  diff_v_edge[dd] = (diff_n0 + diff_n1 + diff_n2) * (n3 + n4 + n5) +
777  (diff_n3 + diff_n4 + diff_n5) * (n0 + n1 + n2);
778  diff_v[dd] = diff_v_tri[dd] * v_edge + v_tri * diff_v_edge[dd];
779  }
780 
781  int shift = ii * P;
782 
783  int jj = 0;
784  int oo = 5;
785  for (; oo <= p; ++oo) {
786 
787  int oo_ksi_eta = 3;
788  for (; oo_ksi_eta <= oo; ++oo_ksi_eta) {
789 
790  int oo_zeta = oo - oo_ksi_eta;
791  if (oo_zeta >= 2) {
792 
793  int k = oo_zeta - 2;
794 
795  int i = 0;
796  for (; i <= (oo_ksi_eta - 3); ++i) {
797  int j = (oo_ksi_eta - 3) - i;
798 
799  if (volumeN != NULL)
800  volumeN[shift + jj] = L0[i] * L1[j] * L2[k] * v;
801 
802  if (diff_volumeN != NULL) {
803  dd = 0;
804  for (; dd < 3; dd++) {
805  diff_volumeN[3 * shift + 3 * jj + dd] =
806  (diffL0[dd * (p + 1) + i] * L1[j] * L2[k] +
807  L0[i] * diffL1[dd * (p + 1) + j] * L2[k] +
808  L0[i] * L1[j] * diffL2[dd * (p + 1) + k]) *
809  v +
810  L0[i] * L1[j] * L2[k] * diff_v[dd];
811  }
812  }
813 
814  ++jj;
815  }
816 
817  }
818  }
819  }
820 
821  if (jj != P)
822  SETERRQ3(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
823  "wrong order %d != %d (%d order)", jj, P, p);
824  }
826 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
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
CHKERRQ(ierr)
#define NBVOLUMEPRISM_H1(P)
Number of base functions on prism for H1 space.
static PetscErrorCode ierr
Definition: h1.c:25
const int N
Definition: speed_test.cpp:3
field with C-1 continuity
Definition: definitions.h:174

◆ 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 419 of file h1.c.

424  {
426 
427  int P = NBVOLUMETET_H1(p);
428  if (P == 0)
430  double diff_ksiL0[3], diff_ksiL1[3], diff_ksiL2[3];
431  int dd = 0;
432  for (; dd < 3; dd++) {
433  diff_ksiL0[dd] = (diffN[1 * 3 + dd] - diffN[0 * 3 + dd]);
434  diff_ksiL1[dd] = (diffN[2 * 3 + dd] - diffN[0 * 3 + dd]);
435  diff_ksiL2[dd] = (diffN[3 * 3 + dd] - diffN[0 * 3 + dd]);
436  }
437  int ii = 0;
438  for (; ii < GDIM; ii++) {
439  int node_shift = ii * 4;
440  double ksiL0 = N[node_shift + 1] - N[node_shift + 0];
441  double ksiL1 = N[node_shift + 2] - N[node_shift + 0];
442  double ksiL2 = N[node_shift + 3] - N[node_shift + 0];
443  double L0[p + 1], L1[p + 1], L2[p + 1];
444  double diffL0[3 * (p + 1)], diffL1[3 * (p + 1)], diffL2[3 * (p + 1)];
445  ierr = base_polynomials(p, ksiL0, diff_ksiL0, L0, diffL0, 3);
446  CHKERRQ(ierr);
447  ierr = base_polynomials(p, ksiL1, diff_ksiL1, L1, diffL1, 3);
448  CHKERRQ(ierr);
449  ierr = base_polynomials(p, ksiL2, diff_ksiL2, L2, diffL2, 3);
450  CHKERRQ(ierr);
451  double v = N[node_shift + 0] * N[node_shift + 1] * N[node_shift + 2] *
452  N[node_shift + 3];
453  double v2[3] = {0, 0, 0};
454  dd = 0;
455  for (; dd < 3; dd++) {
456  v2[dd] = diffN[3 * 0 + dd] * N[node_shift + 1] * N[node_shift + 2] *
457  N[node_shift + 3] +
458  N[node_shift + 0] * diffN[3 * 1 + dd] * N[node_shift + 2] *
459  N[node_shift + 3] +
460  N[node_shift + 0] * N[node_shift + 1] * diffN[3 * 2 + dd] *
461  N[node_shift + 3] +
462  N[node_shift + 0] * N[node_shift + 1] * N[node_shift + 2] *
463  diffN[3 * 3 + dd];
464  }
465  int shift = ii * P;
466  int jj = 0;
467  int oo = 0;
468  for (; oo <= (p - 4); oo++) {
469  int pp0 = 0;
470  for (; pp0 <= oo; pp0++) {
471  int pp1 = 0;
472  for (; (pp0 + pp1) <= oo; pp1++) {
473  int pp2 = oo - pp0 - pp1;
474  if (pp2 >= 0) {
475  if (volumeN != NULL) {
476  volumeN[shift + jj] = L0[pp0] * L1[pp1] * L2[pp2] * v;
477  }
478  if (diff_volumeN != NULL) {
479  dd = 0;
480  for (; dd < 3; dd++) {
481  diff_volumeN[3 * shift + 3 * jj + dd] =
482  (diffL0[dd * (p + 1) + pp0] * L1[pp1] * L2[pp2] +
483  L0[pp0] * diffL1[dd * (p + 1) + pp1] * L2[pp2] +
484  L0[pp0] * L1[pp1] * diffL2[dd * (p + 1) + pp2]) *
485  v;
486  diff_volumeN[3 * shift + 3 * jj + dd] +=
487  L0[pp0] * L1[pp1] * L2[pp2] * v2[dd];
488  }
489  }
490  jj++;
491  }
492  }
493  }
494  }
495  if (jj != P)
496  SETERRQ1(PETSC_COMM_SELF, 1, "wrong order %d", jj);
497  }
499 }
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
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
CHKERRQ(ierr)
static PetscErrorCode ierr
Definition: h1.c:25
const int N
Definition: speed_test.cpp:3
field with C-1 continuity
Definition: definitions.h:174

Variable Documentation

◆ ierr

PetscErrorCode ierr
static

MoFEM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with MoFEM. If not, see http://www.gnu.org/licenses/

Definition at line 25 of file h1.c.