v0.14.0
Loading...
Searching...
No Matches
Macros | Functions | Variables
fem_tools.c File Reference

Loose implementation of some useful functions. More...

#include <definitions.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <complex.h>
#include <assert.h>
#include <string.h>
#include <fem_tools.h>
#include <gm_rule.h>
#include <h1_hdiv_hcurl_l2.h>

Go to the source code of this file.

Macros

#define N_MBTETQ0(x, y, z)   ((2. * (1. - x - y - z) - 1.) * (1. - x - y - z))
 
#define N_MBTETQ1(x, y, z)   ((2. * x - 1.) * x)
 
#define N_MBTETQ2(x, y, z)   ((2. * y - 1.) * y)
 
#define N_MBTETQ3(x, y, z)   ((2. * z - 1.) * z)
 
#define N_MBTETQ4(x, y, z)   (4. * (1. - x - y - z) * x)
 
#define N_MBTETQ5(x, y, z)   (4. * x * y)
 
#define N_MBTETQ6(x, y, z)   (4. * (1. - x - y - z) * y)
 
#define N_MBTETQ7(x, y, z)   (4. * (1. - x - y - z) * z)
 
#define N_MBTETQ8(x, y, z)   (4. * x * z)
 
#define N_MBTETQ9(x, y, z)   (4. * y * z)
 
#define diffN_MBTETQ0x(x, y, z)   (-3. + 4. * x + 4. * y + 4. * z)
 
#define diffN_MBTETQ0y(x, y, z)   (-3. + 4. * x + 4. * y + 4. * z)
 
#define diffN_MBTETQ0z(x, y, z)   (-3. + 4. * x + 4. * y + 4. * z)
 
#define diffN_MBTETQ1x(x, y, z)   (4. * x - 1.)
 
#define diffN_MBTETQ1y(x, y, z)   (0.)
 
#define diffN_MBTETQ1z(x, y, z)   (0.)
 
#define diffN_MBTETQ2x(x, y, z)   (0.)
 
#define diffN_MBTETQ2y(x, y, z)   (4. * y - 1.)
 
#define diffN_MBTETQ2z(x, y, z)   (0.)
 
#define diffN_MBTETQ3x(x, y, z)   (0.)
 
#define diffN_MBTETQ3y(x, y, z)   (0.)
 
#define diffN_MBTETQ3z(x, y, z)   (4. * z - 1.)
 
#define diffN_MBTETQ4x(x, y, z)   (-8. * x + 4. - 4. * y - 4. * z)
 
#define diffN_MBTETQ4y(x, y, z)   (-4. * x)
 
#define diffN_MBTETQ4z(x, y, z)   (-4. * x)
 
#define diffN_MBTETQ5x(x, y, z)   (4. * y)
 
#define diffN_MBTETQ5y(x, y, z)   (4. * x)
 
#define diffN_MBTETQ5z(x, y, z)   (0.)
 
#define diffN_MBTETQ6x(x, y, z)   (-4. * y)
 
#define diffN_MBTETQ6y(x, y, z)   (-8. * y + 4. - 4. * x - 4. * z)
 
#define diffN_MBTETQ6z(x, y, z)   (-4. * y)
 
#define diffN_MBTETQ7x(x, y, z)   (-4. * z)
 
#define diffN_MBTETQ7y(x, y, z)   (-4. * z)
 
#define diffN_MBTETQ7z(x, y, z)   (-8. * z + 4. - 4. * x - 4. * y)
 
#define diffN_MBTETQ8x(x, y, z)   (4. * z)
 
#define diffN_MBTETQ8y(x, y, z)   (0.)
 
#define diffN_MBTETQ8z(x, y, z)   (4. * x)
 
#define diffN_MBTETQ9x(x, y, z)   (0.)
 
#define diffN_MBTETQ9y(x, y, z)   (4. * z)
 
#define diffN_MBTETQ9z(x, y, z)   (4. * y)
 

Functions

double ShapeDetJacVolume (double *jac)
 determined of jacobian
 
PetscErrorCode ShapeInvJacVolume (double *jac)
 
PetscErrorCode Grundmann_Moeller_integration_points_1D_EDGE (int rule, double *G_TRI_X, double *G_TRI_W)
 Compute weights and integration points for edge using Grundmann_Moeller rule.
 
PetscErrorCode Grundmann_Moeller_integration_points_2D_TRI (int rule, double *G_TRI_X, double *G_TRI_Y, double *G_TRI_W)
 
PetscErrorCode Grundmann_Moeller_integration_points_3D_TET (int rule, double *G_TET_X, double *G_TET_Y, double *G_TET_Z, double *G_TET_W)
 
PetscErrorCode ShapeMBTRI (double *N, const double *X, const double *Y, const int G_DIM)
 calculate shape functions on triangle
 
PetscErrorCode ShapeDiffMBTRI (double *diffN)
 calculate derivatives of shape functions
 
PetscErrorCode ShapeFaceBaseMBTRI (double *diffN, const double *coords, double *normal, double *s1, double *s2)
 
PetscErrorCode ShapeFaceNormalMBTRI (double *diffN, const double *coords, double *normal)
 
PetscErrorCode ShapeFaceDiffNormalMBTRI (double *diffN, const double *coords, double *diff_normal)
 calculate derivative of normal in respect to nodal positions
 
PetscErrorCode ShapeJacMBTET (double *diffN, const double *coords, double *jac)
 calculate jacobian
 
double ShapeVolumeMBTET (double *diffN, const double *coords)
 calculate TET volume
 
PetscErrorCode ShapeMBTET (double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
 calculate shape functions
 
PetscErrorCode ShapeDiffMBTET (double *diffN)
 calculate derivatives of shape functions
 
PetscErrorCode ShapeMBTET_inverse (double *N, double *diffN, const double *elem_coords, const double *glob_coords, double *loc_coords)
 calculate local coordinates for given global coordinates
 
PetscErrorCode ShapeMBTRI_inverse (double *N, double *diffN, const double *elem_coords, const double *glob_coords, double *loc_coords)
 calculate local coordinates of triangle element for given global coordinates in 2D (Assume e.g. z=0)
 
PetscErrorCode ShapeDiffMBTETinvJ (double *diffN, double *invJac, double *diffNinvJac)
 calculate shape functions derivatives in space
 
PetscErrorCode GradientOfDeformation (double *diffN, double *dofs, double *F)
 calculate gradient of deformation
 
void ShapeDiffMBTETinvJ_complex (double *diffN, __CLPK_doublecomplex *invJac, __CLPK_doublecomplex *diffNinvJac, enum CBLAS_TRANSPOSE Trans)
 
PetscErrorCode ShapeFaceNormalMBTRI_complex (double *diffN, __CLPK_doublecomplex *xcoords, __CLPK_doublecomplex *xnormal)
 
PetscErrorCode MakeComplexTensor (double *reA, double *imA, __CLPK_doublecomplex *xA)
 
PetscErrorCode InvertComplexGradient (__CLPK_doublecomplex *xF)
 
PetscErrorCode InvertComplexSymmMatrix3by3 (__CLPK_doublecomplex *xC)
 
PetscErrorCode DeterminantComplexGradient (__CLPK_doublecomplex *xF, __CLPK_doublecomplex *det_xF)
 
PetscErrorCode Spin (double *spinOmega, double *vecOmega)
 calculate spin matrix from vector
 
PetscErrorCode make_complex_matrix (double *reA, double *imA, __CLPK_doublecomplex *xA)
 Compose complex matrix (3x3) from two real matrices.
 
PetscErrorCode Normal_hierarchical (int order_approx, int *order_edge_approx, int order, int *order_edge, double *diffN, double *diffN_face, double *diffN_edge[], double *dofs, double *dofs_edge[], double *dofs_face, double *idofs, double *idofs_edge[], double *idofs_face, __CLPK_doublecomplex *xnormal, __CLPK_doublecomplex *xs1, __CLPK_doublecomplex *xs2, int gg)
 Complex normal.
 
PetscErrorCode Base_scale (__CLPK_doublecomplex *xnormal, __CLPK_doublecomplex *xs1, __CLPK_doublecomplex *xs2)
 
PetscErrorCode ShapeMBEDGE (double *N, const double *G_X, int DIM)
 
PetscErrorCode ShapeDiffMBEDGE (double *diffN)
 
PetscErrorCode ShapeMBTRIQ (double *N, const double *X, const double *Y, const int G_DIM)
 
PetscErrorCode ShapeDiffMBTRIQ (double *diffN, const double *X, const double *Y, const int G_DIM)
 
PetscErrorCode ShapeMBTETQ (double *N, const double x, const double y, const double z)
 
PetscErrorCode ShapeDiffMBTETQ (double *diffN, const double x, const double y, const double z)
 
PetscErrorCode ShapeMBTETQ_GAUSS (double *N, const double *X, const double *Y, const double *Z, const int G_DIM)
 
PetscErrorCode ShapeDiffMBTETQ_GAUSS (double *diffN, const double *X, const double *Y, const double *Z, const int G_DIM)
 
PetscErrorCode ShapeJacMBTETQ (const double *diffN, const double *coords, double *Jac)
 
PetscErrorCode ShapeMBTETQ_detJac_at_Gauss_Points (double *detJac_at_Gauss_Points, const double *diffN, const double *coords, int G_DIM)
 
double ShapeVolumeMBTETQ (const double *diffN, const double *coords, int G_DIM, double *G_TET_W)
 
PetscErrorCode ShapeMBTETQ_inverse (double *N, double *diffN, const double *elem_coords, const double *glob_coords, double *loc_coords, const double eps)
 

Variables

static PetscErrorCode ierr
 

Detailed Description

Loose implementation of some useful functions.

Definition in file fem_tools.c.

Macro Definition Documentation

◆ diffN_MBTETQ0x

#define diffN_MBTETQ0x ( x,
y,
z )   (-3. + 4. * x + 4. * y + 4. * z)

Definition at line 828 of file fem_tools.c.

◆ diffN_MBTETQ0y

#define diffN_MBTETQ0y ( x,
y,
z )   (-3. + 4. * x + 4. * y + 4. * z)

Definition at line 829 of file fem_tools.c.

◆ diffN_MBTETQ0z

#define diffN_MBTETQ0z ( x,
y,
z )   (-3. + 4. * x + 4. * y + 4. * z)

Definition at line 830 of file fem_tools.c.

◆ diffN_MBTETQ1x

#define diffN_MBTETQ1x ( x,
y,
z )   (4. * x - 1.)

Definition at line 831 of file fem_tools.c.

◆ diffN_MBTETQ1y

#define diffN_MBTETQ1y ( x,
y,
z )   (0.)

Definition at line 832 of file fem_tools.c.

◆ diffN_MBTETQ1z

#define diffN_MBTETQ1z ( x,
y,
z )   (0.)

Definition at line 833 of file fem_tools.c.

◆ diffN_MBTETQ2x

#define diffN_MBTETQ2x ( x,
y,
z )   (0.)

Definition at line 834 of file fem_tools.c.

◆ diffN_MBTETQ2y

#define diffN_MBTETQ2y ( x,
y,
z )   (4. * y - 1.)

Definition at line 835 of file fem_tools.c.

◆ diffN_MBTETQ2z

#define diffN_MBTETQ2z ( x,
y,
z )   (0.)

Definition at line 836 of file fem_tools.c.

◆ diffN_MBTETQ3x

#define diffN_MBTETQ3x ( x,
y,
z )   (0.)

Definition at line 837 of file fem_tools.c.

◆ diffN_MBTETQ3y

#define diffN_MBTETQ3y ( x,
y,
z )   (0.)

Definition at line 838 of file fem_tools.c.

◆ diffN_MBTETQ3z

#define diffN_MBTETQ3z ( x,
y,
z )   (4. * z - 1.)

Definition at line 839 of file fem_tools.c.

◆ diffN_MBTETQ4x

#define diffN_MBTETQ4x ( x,
y,
z )   (-8. * x + 4. - 4. * y - 4. * z)

Definition at line 840 of file fem_tools.c.

◆ diffN_MBTETQ4y

#define diffN_MBTETQ4y ( x,
y,
z )   (-4. * x)

Definition at line 841 of file fem_tools.c.

◆ diffN_MBTETQ4z

#define diffN_MBTETQ4z ( x,
y,
z )   (-4. * x)

Definition at line 842 of file fem_tools.c.

◆ diffN_MBTETQ5x

#define diffN_MBTETQ5x ( x,
y,
z )   (4. * y)

Definition at line 843 of file fem_tools.c.

◆ diffN_MBTETQ5y

#define diffN_MBTETQ5y ( x,
y,
z )   (4. * x)

Definition at line 844 of file fem_tools.c.

◆ diffN_MBTETQ5z

#define diffN_MBTETQ5z ( x,
y,
z )   (0.)

Definition at line 845 of file fem_tools.c.

◆ diffN_MBTETQ6x

#define diffN_MBTETQ6x ( x,
y,
z )   (-4. * y)

Definition at line 846 of file fem_tools.c.

◆ diffN_MBTETQ6y

#define diffN_MBTETQ6y ( x,
y,
z )   (-8. * y + 4. - 4. * x - 4. * z)

Definition at line 847 of file fem_tools.c.

◆ diffN_MBTETQ6z

#define diffN_MBTETQ6z ( x,
y,
z )   (-4. * y)

Definition at line 848 of file fem_tools.c.

◆ diffN_MBTETQ7x

#define diffN_MBTETQ7x ( x,
y,
z )   (-4. * z)

Definition at line 849 of file fem_tools.c.

◆ diffN_MBTETQ7y

#define diffN_MBTETQ7y ( x,
y,
z )   (-4. * z)

Definition at line 850 of file fem_tools.c.

◆ diffN_MBTETQ7z

#define diffN_MBTETQ7z ( x,
y,
z )   (-8. * z + 4. - 4. * x - 4. * y)

Definition at line 851 of file fem_tools.c.

◆ diffN_MBTETQ8x

#define diffN_MBTETQ8x ( x,
y,
z )   (4. * z)

Definition at line 852 of file fem_tools.c.

◆ diffN_MBTETQ8y

#define diffN_MBTETQ8y ( x,
y,
z )   (0.)

Definition at line 853 of file fem_tools.c.

◆ diffN_MBTETQ8z

#define diffN_MBTETQ8z ( x,
y,
z )   (4. * x)

Definition at line 854 of file fem_tools.c.

◆ diffN_MBTETQ9x

#define diffN_MBTETQ9x ( x,
y,
z )   (0.)

Definition at line 855 of file fem_tools.c.

◆ diffN_MBTETQ9y

#define diffN_MBTETQ9y ( x,
y,
z )   (4. * z)

Definition at line 856 of file fem_tools.c.

◆ diffN_MBTETQ9z

#define diffN_MBTETQ9z ( x,
y,
z )   (4. * y)

Definition at line 857 of file fem_tools.c.

◆ N_MBTETQ0

#define N_MBTETQ0 ( x,
y,
z )   ((2. * (1. - x - y - z) - 1.) * (1. - x - y - z))

Definition at line 818 of file fem_tools.c.

◆ N_MBTETQ1

#define N_MBTETQ1 ( x,
y,
z )   ((2. * x - 1.) * x)

Definition at line 819 of file fem_tools.c.

◆ N_MBTETQ2

#define N_MBTETQ2 ( x,
y,
z )   ((2. * y - 1.) * y)

Definition at line 820 of file fem_tools.c.

◆ N_MBTETQ3

#define N_MBTETQ3 ( x,
y,
z )   ((2. * z - 1.) * z)

Definition at line 821 of file fem_tools.c.

◆ N_MBTETQ4

#define N_MBTETQ4 ( x,
y,
z )   (4. * (1. - x - y - z) * x)

Definition at line 822 of file fem_tools.c.

◆ N_MBTETQ5

#define N_MBTETQ5 ( x,
y,
z )   (4. * x * y)

Definition at line 823 of file fem_tools.c.

◆ N_MBTETQ6

#define N_MBTETQ6 ( x,
y,
z )   (4. * (1. - x - y - z) * y)

Definition at line 824 of file fem_tools.c.

◆ N_MBTETQ7

#define N_MBTETQ7 ( x,
y,
z )   (4. * (1. - x - y - z) * z)

Definition at line 825 of file fem_tools.c.

◆ N_MBTETQ8

#define N_MBTETQ8 ( x,
y,
z )   (4. * x * z)

Definition at line 826 of file fem_tools.c.

◆ N_MBTETQ9

#define N_MBTETQ9 ( x,
y,
z )   (4. * y * z)

Definition at line 827 of file fem_tools.c.

Function Documentation

◆ Base_scale()

PetscErrorCode Base_scale ( __CLPK_doublecomplex * xnormal,
__CLPK_doublecomplex * xs1,
__CLPK_doublecomplex * xs2 )

Definition at line 741 of file fem_tools.c.

743 {
745 complex double xnrm2_normal = csqrt(cpow(xnormal[0].r + I * xnormal[0].i, 2) +
746 cpow(xnormal[1].r + I * xnormal[1].i, 2) +
747 cpow(xnormal[2].r + I * xnormal[2].i, 2));
748 int dd = 0;
749 for (; dd < 3; dd++) {
750 complex double s1 = (xs1[dd].r + I * xs1[dd].i) * xnrm2_normal;
751 complex double s2 = (xs2[dd].r + I * xs2[dd].i) * xnrm2_normal;
752 xs1[dd].r = creal(s1);
753 xs1[dd].i = cimag(s1);
754 xs2[dd].r = creal(s2);
755 xs2[dd].i = cimag(s2);
756 }
758}
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
FTensor::Index< 'i', SPACE_DIM > i
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
constexpr IntegrationType I
__CLPK_doublereal i
Definition lapack_wrap.h:35
__CLPK_doublereal r
Definition lapack_wrap.h:35

◆ DeterminantComplexGradient()

PetscErrorCode DeterminantComplexGradient ( __CLPK_doublecomplex * xF,
__CLPK_doublecomplex * det_xF )

Definition at line 526 of file fem_tools.c.

527 {
529 __CLPK_integer IPIV[4];
530 if (lapack_zgetrf(3, 3, xF, 3, IPIV) != 0) {
531 SETERRQ(PETSC_COMM_SELF, 1, "lapack_zgetrf(3,3,xF,3,IPIV) != 0");
532 }
533 double complex det = 1;
534 int i = 0, j = 0;
535 for (; i < 3; i++) {
536 det *= xF[3 * i + i].r + I * xF[3 * i + i].i;
537 if (IPIV[i] != i + 1)
538 j++;
539 }
540 if ((j - (j / 2) * 2) != 0)
541 det = -det;
542 (*det_xF).r = creal(det);
543 (*det_xF).i = cimag(det);
545}
static __CLPK_integer lapack_zgetrf(__CLPK_integer m, __CLPK_integer n, __CLPK_doublecomplex *a, __CLPK_integer lda, __CLPK_integer *ipiv)
long int __CLPK_integer
Definition lapack_wrap.h:23
FTensor::Index< 'j', 3 > j

◆ GradientOfDeformation()

PetscErrorCode GradientOfDeformation ( double * diffN,
double * dofs,
double * F )

calculate gradient of deformation

Definition at line 425 of file fem_tools.c.

425 {
427 int col, row = 0;
428 for (; row < 3; row++)
429 for (col = 0; col < 3; col++) {
430 F[3 * row + col] = cblas_ddot(4, &diffN[col], 3, &dofs[row], 3);
431 }
433}
@ F

◆ Grundmann_Moeller_integration_points_1D_EDGE()

PetscErrorCode Grundmann_Moeller_integration_points_1D_EDGE ( int rule,
double * G_TRI_X,
double * G_TRI_W )

Compute weights and integration points for edge using Grundmann_Moeller rule.

Definition at line 55 of file fem_tools.c.

57 {
59
60 int dim_num = 1;
61 int point;
62 int point_num;
63 double *w;
64 double *x;
65
66 // GM_RULE_SET determines the weights and abscissas
67 // pof a Grundmann-Moeller quadrature rule for
68 // the DIM_NUM dimensional simplex,
69 // using a rule of in index RULE,
70 // which will have degree of exactness 2*RULE+1.
71
72 // printf ( " Here we use DIM_NUM = %d\n", dim_num );
73 // printf ( " RULE = %d\n", rule );
74 // printf ( " DEGREE = %d\n", 2 * rule + 1 );
75
76 point_num = gm_rule_size(rule, dim_num);
77
78 ierr = PetscMalloc(point_num * sizeof(double), &w);
79 CHKERRQ(ierr);
80 ierr = PetscMalloc(dim_num * point_num * sizeof(double), &x);
81 CHKERRQ(ierr);
82
83 gm_rule_set(rule, dim_num, point_num, w, x);
84
85 for (point = 0; point < point_num; point++) {
86 G_TRI_X[point] = x[0 + point * dim_num];
87 G_TRI_W[point] = w[point];
88 }
89
90 ierr = PetscFree(w);
91 CHKERRQ(ierr);
92 ierr = PetscFree(x);
93 CHKERRQ(ierr);
94
96}
static PetscErrorCode ierr
Definition fem_tools.c:20
void gm_rule_set(int rule, int dim_num, int point_num, double w[], double x[])
Definition gm_rule.c:152
int gm_rule_size(int rule, int dim_num)
Definition gm_rule.c:294
double w(const double g, const double t)

◆ Grundmann_Moeller_integration_points_2D_TRI()

PetscErrorCode Grundmann_Moeller_integration_points_2D_TRI ( int rule,
double * G_TRI_X,
double * G_TRI_Y,
double * G_TRI_W )

Compute weights and integration points for 2D Triangle using Grundmann_Moeller rule

Definition at line 98 of file fem_tools.c.

101 {
103
104 int dim_num = 2;
105 int point;
106 int point_num;
107 double *w;
108 double *x;
109
110 // GM_RULE_SET determines the weights and abscissas
111 // pof a Grundmann-Moeller quadrature rule for
112 // the DIM_NUM dimensional simplex,
113 // using a rule of in index RULE,
114 // which will have degree of exactness 2*RULE+1.
115
116 // printf ( " Here we use DIM_NUM = %d\n", dim_num );
117 // printf ( " RULE = %d\n", rule );
118 // printf ( " DEGREE = %d\n", 2 * rule + 1 );
119
120 point_num = gm_rule_size(rule, dim_num);
121
122 ierr = PetscMalloc(point_num * sizeof(double), &w);
123 CHKERRQ(ierr);
124 ierr = PetscMalloc(dim_num * point_num * sizeof(double), &x);
125 CHKERRQ(ierr);
126
127 gm_rule_set(rule, dim_num, point_num, w, x);
128
129 for (point = 0; point < point_num; point++) {
130 G_TRI_X[point] = x[0 + point * dim_num];
131 G_TRI_Y[point] = x[1 + point * dim_num];
132 G_TRI_W[point] = w[point];
133 }
134
135 ierr = PetscFree(w);
136 CHKERRQ(ierr);
137 ierr = PetscFree(x);
138 CHKERRQ(ierr);
139
141}

◆ Grundmann_Moeller_integration_points_3D_TET()

PetscErrorCode Grundmann_Moeller_integration_points_3D_TET ( int rule,
double * G_TET_X,
double * G_TET_Y,
double * G_TET_Z,
double * G_TET_W )

Compute weights and integration points for 3D Tet using Grundmann_Moeller rule

Definition at line 143 of file fem_tools.c.

147 {
149
150 int dim_num = 3;
151 int point;
152 int point_num;
153 double *w;
154 double *x;
155
156 // printf ( " Here we use DIM_NUM = %d\n", dim_num );
157 // printf ( " RULE = %d\n", rule );
158 // printf ( " DEGREE = %d\n", 2 * rule + 1 );
159
160 point_num = gm_rule_size(rule, dim_num);
161
162 ierr = PetscMalloc(point_num * sizeof(double), &w);
163 CHKERRQ(ierr);
164 ierr = PetscMalloc(dim_num * point_num * sizeof(double), &x);
165 CHKERRQ(ierr);
166
167 gm_rule_set(rule, dim_num, point_num, w, x);
168 for (point = 0; point < point_num; point++) {
169 G_TET_X[point] = x[0 + point * dim_num];
170 G_TET_Y[point] = x[1 + point * dim_num];
171 G_TET_Z[point] = x[2 + point * dim_num];
172 G_TET_W[point] = w[point];
173 }
174
175 ierr = PetscFree(w);
176 CHKERRQ(ierr);
177 ierr = PetscFree(x);
178 CHKERRQ(ierr);
179
181}

◆ InvertComplexGradient()

PetscErrorCode InvertComplexGradient ( __CLPK_doublecomplex * xF)

Definition at line 499 of file fem_tools.c.

499 {
501 __CLPK_integer IPIV[4];
502 __CLPK_doublecomplex WORK[4];
503 __CLPK_integer LWORK = 4;
504 __CLPK_integer info;
505 info = lapack_zgetrf(3, 3, xF, 3, IPIV);
506 if (info != 0)
507 SETERRQ(PETSC_COMM_SELF, 1, "info == 0");
508 info = lapack_zgetri(3, xF, 3, IPIV, WORK, LWORK);
509 if (info != 0)
510 SETERRQ(PETSC_COMM_SELF, 1, "info == 0");
512}
static __CLPK_integer lapack_zgetri(__CLPK_integer n, __CLPK_doublecomplex *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublecomplex *work, __CLPK_integer lwork)

◆ InvertComplexSymmMatrix3by3()

PetscErrorCode InvertComplexSymmMatrix3by3 ( __CLPK_doublecomplex * xC)

Definition at line 513 of file fem_tools.c.

513 {
515 __CLPK_integer info;
516 info = lapack_zpotrf('L', 3, xC, 3);
517 if (info == 0)
518 SETERRQ(PETSC_COMM_SELF, 1, "info == 0");
519 // assert(info == 0);
520 info = lapack_zpotri('L', 3, xC, 3);
521 if (info == 0)
522 SETERRQ(PETSC_COMM_SELF, 1, "info == 0");
523 // assert(info == 0);
525}
static __CLPK_integer lapack_zpotri(char uplo, __CLPK_integer n, __CLPK_doublecomplex *a, __CLPK_integer lda)
static __CLPK_integer lapack_zpotrf(char uplo, __CLPK_integer n, __CLPK_doublecomplex *a, __CLPK_integer lda)

◆ make_complex_matrix()

PetscErrorCode make_complex_matrix ( double * reA,
double * imA,
__CLPK_doublecomplex * xA )

Compose complex matrix (3x3) from two real matrices.

Definition at line 557 of file fem_tools.c.

558 {
560 int ii = 0, jj;
561 for (; ii < 3; ii++) {
562 for (jj = 0; jj < 3; jj++) {
563 xA[3 * ii + jj].r = reA[3 * ii + jj];
564 xA[3 * ii + jj].i = imA[3 * ii + jj];
565 }
566 }
568}

◆ MakeComplexTensor()

PetscErrorCode MakeComplexTensor ( double * reA,
double * imA,
__CLPK_doublecomplex * xA )

Definition at line 487 of file fem_tools.c.

488 {
490 int ii = 0, jj;
491 for (; ii < 3; ii++) {
492 for (jj = 0; jj < 3; jj++) {
493 xA[3 * ii + jj].r = reA[3 * ii + jj];
494 xA[3 * ii + jj].i = imA[3 * ii + jj];
495 }
496 }
498}

◆ Normal_hierarchical()

PetscErrorCode Normal_hierarchical ( int order_approx,
int * order_edge_approx,
int order,
int * order_edge,
double * diffN,
double * diffN_face,
double * diffN_edge[],
double * dofs,
double * dofs_edge[],
double * dofs_face,
double * idofs,
double * idofs_edge[],
double * idofs_face,
__CLPK_doublecomplex * xnormal,
__CLPK_doublecomplex * xs1,
__CLPK_doublecomplex * xs2,
int gg )

Complex normal.

Definition at line 569 of file fem_tools.c.

574 {
576 int nn, ee, dd;
577 // node
578 double complex diffX_x_node, diffX_y_node, diffX_z_node;
579 double complex diffY_x_node, diffY_y_node, diffY_z_node;
580 diffX_x_node = 0.;
581 diffX_y_node = 0.;
582 diffX_z_node = 0.;
583 diffY_x_node = 0.;
584 diffY_y_node = 0.;
585 diffY_z_node = 0.;
586 if (dofs != NULL || idofs != NULL) {
587 nn = 0;
588 for (; nn < 3; nn++) {
589 if (dofs != NULL) {
590 diffX_x_node += dofs[3 * nn + 0] * diffN[2 * nn + 0];
591 diffX_y_node += dofs[3 * nn + 1] * diffN[2 * nn + 0];
592 diffX_z_node += dofs[3 * nn + 2] * diffN[2 * nn + 0];
593 diffY_x_node += dofs[3 * nn + 0] * diffN[2 * nn + 1];
594 diffY_y_node += dofs[3 * nn + 1] * diffN[2 * nn + 1];
595 diffY_z_node += dofs[3 * nn + 2] * diffN[2 * nn + 1];
596 }
597 if (idofs != NULL) {
598 diffX_x_node += I * idofs[3 * nn + 0] * diffN[2 * nn + 0];
599 diffX_y_node += I * idofs[3 * nn + 1] * diffN[2 * nn + 0];
600 diffX_z_node += I * idofs[3 * nn + 2] * diffN[2 * nn + 0];
601 diffY_x_node += I * idofs[3 * nn + 0] * diffN[2 * nn + 1];
602 diffY_y_node += I * idofs[3 * nn + 1] * diffN[2 * nn + 1];
603 diffY_z_node += I * idofs[3 * nn + 2] * diffN[2 * nn + 1];
604 }
605 }
606 }
607 double complex diffX_x, diffX_y, diffX_z;
608 double complex diffY_x, diffY_y, diffY_z;
609 diffX_x = diffX_x_node;
610 diffX_y = diffX_y_node;
611 diffX_z = diffX_z_node;
612 diffY_x = diffY_x_node;
613 diffY_y = diffY_y_node;
614 diffY_z = diffY_z_node;
615 if (dofs_face != NULL || idofs_face != NULL) {
616 int nb_dofs_face = NBFACETRI_H1(order);
617 int nb_dofs_approx_face = NBFACETRI_H1(order_approx);
618 if (nb_dofs_face > 0) {
619 if (dofs_face != NULL) {
620 diffX_x += cblas_ddot(nb_dofs_face, &dofs_face[0], 3,
621 &diffN_face[gg * 2 * nb_dofs_approx_face + 0], 2);
622 diffX_y += cblas_ddot(nb_dofs_face, &dofs_face[1], 3,
623 &diffN_face[gg * 2 * nb_dofs_approx_face + 0], 2);
624 diffX_z += cblas_ddot(nb_dofs_face, &dofs_face[2], 3,
625 &diffN_face[gg * 2 * nb_dofs_approx_face + 0], 2);
626 diffY_x += cblas_ddot(nb_dofs_face, &dofs_face[0], 3,
627 &diffN_face[gg * 2 * nb_dofs_approx_face + 1], 2);
628 diffY_y += cblas_ddot(nb_dofs_face, &dofs_face[1], 3,
629 &diffN_face[gg * 2 * nb_dofs_approx_face + 1], 2);
630 diffY_z += cblas_ddot(nb_dofs_face, &dofs_face[2], 3,
631 &diffN_face[gg * 2 * nb_dofs_approx_face + 1], 2);
632 }
633 if (idofs_face != NULL) {
634 diffX_x +=
635 I * cblas_ddot(nb_dofs_face, &idofs_face[0], 3,
636 &diffN_face[gg * 2 * nb_dofs_approx_face + 0], 2);
637 diffX_y +=
638 I * cblas_ddot(nb_dofs_face, &idofs_face[1], 3,
639 &diffN_face[gg * 2 * nb_dofs_approx_face + 0], 2);
640 diffX_z +=
641 I * cblas_ddot(nb_dofs_face, &idofs_face[2], 3,
642 &diffN_face[gg * 2 * nb_dofs_approx_face + 0], 2);
643 diffY_x +=
644 I * cblas_ddot(nb_dofs_face, &idofs_face[0], 3,
645 &diffN_face[gg * 2 * nb_dofs_approx_face + 1], 2);
646 diffY_y +=
647 I * cblas_ddot(nb_dofs_face, &idofs_face[1], 3,
648 &diffN_face[gg * 2 * nb_dofs_approx_face + 1], 2);
649 diffY_z +=
650 I * cblas_ddot(nb_dofs_face, &idofs_face[2], 3,
651 &diffN_face[gg * 2 * nb_dofs_approx_face + 1], 2);
652 }
653 }
654 }
655 ee = 0;
656 if (dofs_edge != NULL || idofs_edge != NULL) {
657 for (; ee < 3; ee++) {
658 int nb_dofs_edge = NBEDGE_H1(order_edge[ee]);
659 int nb_dofs_approx_edge = NBEDGE_H1(order_edge_approx[ee]);
660 if (nb_dofs_edge > 0) {
661 if (dofs_edge != NULL) {
662 if (dofs_edge[ee] != NULL) {
663 diffX_x += cblas_ddot(
664 nb_dofs_edge, &(dofs_edge[ee])[0], 3,
665 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 0], 2);
666 diffX_y += cblas_ddot(
667 nb_dofs_edge, &(dofs_edge[ee])[1], 3,
668 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 0], 2);
669 diffX_z += cblas_ddot(
670 nb_dofs_edge, &(dofs_edge[ee])[2], 3,
671 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 0], 2);
672 diffY_x += cblas_ddot(
673 nb_dofs_edge, &(dofs_edge[ee])[0], 3,
674 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 1], 2);
675 diffY_y += cblas_ddot(
676 nb_dofs_edge, &(dofs_edge[ee])[1], 3,
677 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 1], 2);
678 diffY_z += cblas_ddot(
679 nb_dofs_edge, &(dofs_edge[ee])[2], 3,
680 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 1], 2);
681 }
682 }
683 if (idofs_edge != NULL) {
684 if (idofs_edge[ee] == NULL)
685 continue;
686 diffX_x +=
687 I * cblas_ddot(
688 nb_dofs_edge, &(idofs_edge[ee])[0], 3,
689 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 0], 2);
690 diffX_y +=
691 I * cblas_ddot(
692 nb_dofs_edge, &(idofs_edge[ee])[1], 3,
693 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 0], 2);
694 diffX_z +=
695 I * cblas_ddot(
696 nb_dofs_edge, &(idofs_edge[ee])[2], 3,
697 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 0], 2);
698 diffY_x +=
699 I * cblas_ddot(
700 nb_dofs_edge, &(idofs_edge[ee])[0], 3,
701 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 1], 2);
702 diffY_y +=
703 I * cblas_ddot(
704 nb_dofs_edge, &(idofs_edge[ee])[1], 3,
705 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 1], 2);
706 diffY_z +=
707 I * cblas_ddot(
708 nb_dofs_edge, &(idofs_edge[ee])[2], 3,
709 &(diffN_edge[ee])[gg * 2 * nb_dofs_approx_edge + 1], 2);
710 }
711 }
712 }
713 }
714 double complex normal[3];
715 normal[0] = diffX_y * diffY_z - diffX_z * diffY_y;
716 normal[1] = diffX_z * diffY_x - diffX_x * diffY_z;
717 normal[2] = diffX_x * diffY_y - diffX_y * diffY_x;
718 dd = 0;
719 for (; dd < 3; dd++) {
720 xnormal[dd].r = creal(normal[dd]);
721 xnormal[dd].i = cimag(normal[dd]);
722 }
723 if (xs1 != NULL) {
724 xs1[0].r = creal(diffX_x);
725 xs1[0].i = cimag(diffX_x);
726 xs1[1].r = creal(diffX_y);
727 xs1[1].i = cimag(diffX_y);
728 xs1[2].r = creal(diffX_z);
729 xs1[2].i = cimag(diffX_z);
730 }
731 if (xs2 != NULL) {
732 xs2[0].r = creal(diffY_x);
733 xs2[0].i = cimag(diffY_x);
734 xs2[1].r = creal(diffY_y);
735 xs2[1].i = cimag(diffY_y);
736 xs2[2].r = creal(diffY_z);
737 xs2[2].i = cimag(diffY_z);
738 }
740}
constexpr int order
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.

◆ ShapeDetJacVolume()

double ShapeDetJacVolume ( double * jac)

determined of jacobian

Definition at line 22 of file fem_tools.c.

22 {
23 double det_jac;
24 __CLPK_integer ipiv[4];
25 __CLPK_integer info = lapack_dgetrf(3, 3, jac, 3, ipiv);
26 if (info != 0)
27 return -1;
28 int i = 0, j = 0;
29 det_jac = 1.;
30 for (; i < 3; i++) {
31 det_jac *= jac[3 * i + i];
32 if (ipiv[i] != i + 1)
33 j++;
34 }
35 if ((j - (j / 2) * 2) != 0)
36 det_jac = -det_jac;
37 return det_jac;
38}
static __CLPK_integer lapack_dgetrf(__CLPK_integer m, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv)

◆ ShapeDiffMBEDGE()

PetscErrorCode ShapeDiffMBEDGE ( double * diffN)

Definition at line 771 of file fem_tools.c.

771 {
773 diffN[0] = diffN_MBEDGE0;
774 diffN[1] = diffN_MBEDGE1;
776}
#define diffN_MBEDGE0
derivative of edge shape function
Definition fem_tools.h:107
#define diffN_MBEDGE1
derivative of edge shape function
Definition fem_tools.h:108

◆ ShapeDiffMBTET()

PetscErrorCode ShapeDiffMBTET ( double * diffN)

calculate derivatives of shape functions

Definition at line 319 of file fem_tools.c.

319 {
321 diffN[0] = diffN_MBTET0x;
322 diffN[1] = diffN_MBTET0y;
323 diffN[2] = diffN_MBTET0z;
324 diffN[3] = diffN_MBTET1x;
325 diffN[4] = diffN_MBTET1y;
326 diffN[5] = diffN_MBTET1z;
327 diffN[6] = diffN_MBTET2x;
328 diffN[7] = diffN_MBTET2y;
329 diffN[8] = diffN_MBTET2z;
330 diffN[9] = diffN_MBTET3x;
331 diffN[10] = diffN_MBTET3y;
332 diffN[11] = diffN_MBTET3z;
334}
#define diffN_MBTET0z
derivative of tetrahedral shape function
Definition fem_tools.h:34
#define diffN_MBTET2z
derivative of tetrahedral shape function
Definition fem_tools.h:40
#define diffN_MBTET3y
derivative of tetrahedral shape function
Definition fem_tools.h:42
#define diffN_MBTET0y
derivative of tetrahedral shape function
Definition fem_tools.h:33
#define diffN_MBTET1z
derivative of tetrahedral shape function
Definition fem_tools.h:37
#define diffN_MBTET3z
derivative of tetrahedral shape function
Definition fem_tools.h:43
#define diffN_MBTET1y
derivative of tetrahedral shape function
Definition fem_tools.h:36
#define diffN_MBTET2y
derivative of tetrahedral shape function
Definition fem_tools.h:39
#define diffN_MBTET2x
derivative of tetrahedral shape function
Definition fem_tools.h:38
#define diffN_MBTET3x
derivative of tetrahedral shape function
Definition fem_tools.h:41
#define diffN_MBTET1x
derivative of tetrahedral shape function
Definition fem_tools.h:35
#define diffN_MBTET0x
derivative of tetrahedral shape function
Definition fem_tools.h:32

◆ ShapeDiffMBTETinvJ()

PetscErrorCode ShapeDiffMBTETinvJ ( double * diffN,
double * invJac,
double * diffNinvJac )

calculate shape functions derivatives in space

Definition at line 415 of file fem_tools.c.

416 {
418 int ii = 0;
419 for (; ii < 4; ii++) {
420 cblas_dgemv(CblasRowMajor, CblasTrans, 3, 3, 1., invJac, 3, &diffN[ii * 3],
421 1, 0., &diffNinvJac[ii * 3], 1);
422 }
424}
MatrixDouble invJac

◆ ShapeDiffMBTETinvJ_complex()

void ShapeDiffMBTETinvJ_complex ( double * diffN,
__CLPK_doublecomplex * invJac,
__CLPK_doublecomplex * diffNinvJac,
enum CBLAS_TRANSPOSE Trans )

Definition at line 437 of file fem_tools.c.

439 {
440 __CLPK_doublecomplex tmp1 = {1., 0.}, tmp2 = {0., 0.};
441 int ii = 0, jj;
442 for (; ii < 4; ii++) {
443 __CLPK_doublecomplex tmp3[3];
444 for (jj = 0; jj < 3; jj++) {
445 tmp3[jj].r = diffN[ii * 3 + jj];
446 tmp3[jj].i = 0;
447 }
448 cblas_zgemv(CblasRowMajor, Trans, 3, 3, &tmp1, invJac, 3, tmp3, 1, &tmp2,
449 &diffNinvJac[ii * 3], 1);
450 }
451}

◆ ShapeDiffMBTETQ()

PetscErrorCode ShapeDiffMBTETQ ( double * diffN,
const double x,
const double y,
const double z )

Definition at line 873 of file fem_tools.c.

874 {
876 diffN[0] = diffN_MBTETQ0x(x, y, z);
877 diffN[1] = diffN_MBTETQ0y(x, y, z);
878 diffN[2] = diffN_MBTETQ0z(x, y, z);
879 diffN[3] = diffN_MBTETQ1x(x, y, z);
880 diffN[4] = diffN_MBTETQ1y(x, y, z);
881 diffN[5] = diffN_MBTETQ1z(x, y, z);
882 diffN[6] = diffN_MBTETQ2x(x, y, z);
883 diffN[7] = diffN_MBTETQ2y(x, y, z);
884 diffN[8] = diffN_MBTETQ2z(x, y, z);
885 diffN[9] = diffN_MBTETQ3x(x, y, z);
886 diffN[10] = diffN_MBTETQ3y(x, y, z);
887 diffN[11] = diffN_MBTETQ3z(x, y, z);
888 diffN[12] = diffN_MBTETQ4x(x, y, z);
889 diffN[13] = diffN_MBTETQ4y(x, y, z);
890 diffN[14] = diffN_MBTETQ4z(x, y, z);
891 diffN[15] = diffN_MBTETQ5x(x, y, z);
892 diffN[16] = diffN_MBTETQ5y(x, y, z);
893 diffN[17] = diffN_MBTETQ5z(x, y, z);
894 diffN[18] = diffN_MBTETQ6x(x, y, z);
895 diffN[19] = diffN_MBTETQ6y(x, y, z);
896 diffN[20] = diffN_MBTETQ6z(x, y, z);
897 diffN[21] = diffN_MBTETQ7x(x, y, z);
898 diffN[22] = diffN_MBTETQ7y(x, y, z);
899 diffN[23] = diffN_MBTETQ7z(x, y, z);
900 diffN[24] = diffN_MBTETQ8x(x, y, z);
901 diffN[25] = diffN_MBTETQ8y(x, y, z);
902 diffN[26] = diffN_MBTETQ8z(x, y, z);
903 diffN[27] = diffN_MBTETQ9x(x, y, z);
904 diffN[28] = diffN_MBTETQ9y(x, y, z);
905 diffN[29] = diffN_MBTETQ9z(x, y, z);
907}
#define diffN_MBTETQ5x(x, y, z)
Definition fem_tools.c:843
#define diffN_MBTETQ1y(x, y, z)
Definition fem_tools.c:832
#define diffN_MBTETQ1z(x, y, z)
Definition fem_tools.c:833
#define diffN_MBTETQ8z(x, y, z)
Definition fem_tools.c:854
#define diffN_MBTETQ7y(x, y, z)
Definition fem_tools.c:850
#define diffN_MBTETQ3y(x, y, z)
Definition fem_tools.c:838
#define diffN_MBTETQ0y(x, y, z)
Definition fem_tools.c:829
#define diffN_MBTETQ9y(x, y, z)
Definition fem_tools.c:856
#define diffN_MBTETQ9x(x, y, z)
Definition fem_tools.c:855
#define diffN_MBTETQ0x(x, y, z)
Definition fem_tools.c:828
#define diffN_MBTETQ6y(x, y, z)
Definition fem_tools.c:847
#define diffN_MBTETQ6x(x, y, z)
Definition fem_tools.c:846
#define diffN_MBTETQ2z(x, y, z)
Definition fem_tools.c:836
#define diffN_MBTETQ6z(x, y, z)
Definition fem_tools.c:848
#define diffN_MBTETQ4z(x, y, z)
Definition fem_tools.c:842
#define diffN_MBTETQ0z(x, y, z)
Definition fem_tools.c:830
#define diffN_MBTETQ3z(x, y, z)
Definition fem_tools.c:839
#define diffN_MBTETQ4x(x, y, z)
Definition fem_tools.c:840
#define diffN_MBTETQ2y(x, y, z)
Definition fem_tools.c:835
#define diffN_MBTETQ9z(x, y, z)
Definition fem_tools.c:857
#define diffN_MBTETQ5y(x, y, z)
Definition fem_tools.c:844
#define diffN_MBTETQ4y(x, y, z)
Definition fem_tools.c:841
#define diffN_MBTETQ3x(x, y, z)
Definition fem_tools.c:837
#define diffN_MBTETQ2x(x, y, z)
Definition fem_tools.c:834
#define diffN_MBTETQ7x(x, y, z)
Definition fem_tools.c:849
#define diffN_MBTETQ1x(x, y, z)
Definition fem_tools.c:831
#define diffN_MBTETQ8y(x, y, z)
Definition fem_tools.c:853
#define diffN_MBTETQ7z(x, y, z)
Definition fem_tools.c:851
#define diffN_MBTETQ5z(x, y, z)
Definition fem_tools.c:845
#define diffN_MBTETQ8x(x, y, z)
Definition fem_tools.c:852

◆ ShapeDiffMBTETQ_GAUSS()

PetscErrorCode ShapeDiffMBTETQ_GAUSS ( double * diffN,
const double * X,
const double * Y,
const double * Z,
const int G_DIM )

Definition at line 927 of file fem_tools.c.

929 {
931 int ii = 0;
932 for (; ii < G_DIM; ii++) {
933 double x = X[ii], y = Y[ii], z = Z[ii];
934 diffN[30 * ii + 0] = diffN_MBTETQ0x(x, y, z);
935 diffN[30 * ii + 1] = diffN_MBTETQ0y(x, y, z);
936 diffN[30 * ii + 2] = diffN_MBTETQ0z(x, y, z);
937 diffN[30 * ii + 3] = diffN_MBTETQ1x(x, y, z);
938 diffN[30 * ii + 4] = diffN_MBTETQ1y(x, y, z);
939 diffN[30 * ii + 5] = diffN_MBTETQ1z(x, y, z);
940 diffN[30 * ii + 6] = diffN_MBTETQ2x(x, y, z);
941 diffN[30 * ii + 7] = diffN_MBTETQ2y(x, y, z);
942 diffN[30 * ii + 8] = diffN_MBTETQ2z(x, y, z);
943 diffN[30 * ii + 9] = diffN_MBTETQ3x(x, y, z);
944 diffN[30 * ii + 10] = diffN_MBTETQ3y(x, y, z);
945 diffN[30 * ii + 11] = diffN_MBTETQ3z(x, y, z);
946 diffN[30 * ii + 12] = diffN_MBTETQ4x(x, y, z);
947 diffN[30 * ii + 13] = diffN_MBTETQ4y(x, y, z);
948 diffN[30 * ii + 14] = diffN_MBTETQ4z(x, y, z);
949 diffN[30 * ii + 15] = diffN_MBTETQ5x(x, y, z);
950 diffN[30 * ii + 16] = diffN_MBTETQ5y(x, y, z);
951 diffN[30 * ii + 17] = diffN_MBTETQ5z(x, y, z);
952 diffN[30 * ii + 18] = diffN_MBTETQ6x(x, y, z);
953 diffN[30 * ii + 19] = diffN_MBTETQ6y(x, y, z);
954 diffN[30 * ii + 20] = diffN_MBTETQ6z(x, y, z);
955 diffN[30 * ii + 21] = diffN_MBTETQ7x(x, y, z);
956 diffN[30 * ii + 22] = diffN_MBTETQ7y(x, y, z);
957 diffN[30 * ii + 23] = diffN_MBTETQ7z(x, y, z);
958 diffN[30 * ii + 24] = diffN_MBTETQ8x(x, y, z);
959 diffN[30 * ii + 25] = diffN_MBTETQ8y(x, y, z);
960 diffN[30 * ii + 26] = diffN_MBTETQ8z(x, y, z);
961 diffN[30 * ii + 27] = diffN_MBTETQ9x(x, y, z);
962 diffN[30 * ii + 28] = diffN_MBTETQ9y(x, y, z);
963 diffN[30 * ii + 29] = diffN_MBTETQ9z(x, y, z);
964 }
966}

◆ ShapeDiffMBTRI()

PetscErrorCode ShapeDiffMBTRI ( double * diffN)

calculate derivatives of shape functions

Definition at line 194 of file fem_tools.c.

194 {
196 diffN[0] = diffN_MBTRI0x;
197 diffN[1] = diffN_MBTRI0y;
198 diffN[2] = diffN_MBTRI1x;
199 diffN[3] = diffN_MBTRI1y;
200 diffN[4] = diffN_MBTRI2x;
201 diffN[5] = diffN_MBTRI2y;
203}
#define diffN_MBTRI2x
derivative of triangle shape function
Definition fem_tools.h:53
#define diffN_MBTRI0y
derivative of triangle shape function
Definition fem_tools.h:50
#define diffN_MBTRI1y
derivative of triangle shape function
Definition fem_tools.h:52
#define diffN_MBTRI0x
derivative of triangle shape function
Definition fem_tools.h:49
#define diffN_MBTRI2y
derivative of triangle shape function
Definition fem_tools.h:54
#define diffN_MBTRI1x
derivative of triangle shape function
Definition fem_tools.h:51

◆ ShapeDiffMBTRIQ()

PetscErrorCode ShapeDiffMBTRIQ ( double * diffN,
const double * X,
const double * Y,
const int G_DIM )

Definition at line 795 of file fem_tools.c.

796 {
798 int ii = 0;
799 for (; ii < G_DIM; ii++) {
800 double x = X[ii], y = Y[ii];
801 diffN[12 * ii + 0] = diffN_MBTRIQ0x(x, y);
802 diffN[12 * ii + 1] = diffN_MBTRIQ0y(x, y);
803 diffN[12 * ii + 2] = diffN_MBTRIQ1x(x, y);
804 diffN[12 * ii + 3] = diffN_MBTRIQ1y(x, y);
805 diffN[12 * ii + 4] = diffN_MBTRIQ2x(x, y);
806 diffN[12 * ii + 5] = diffN_MBTRIQ2y(x, y);
807 diffN[12 * ii + 6] = diffN_MBTRIQ3x(x, y);
808 diffN[12 * ii + 7] = diffN_MBTRIQ3y(x, y);
809 diffN[12 * ii + 8] = diffN_MBTRIQ4x(x, y);
810 diffN[12 * ii + 9] = diffN_MBTRIQ4y(x, y);
811 diffN[12 * ii + 10] = diffN_MBTRIQ5x(x, y);
812 diffN[12 * ii + 11] = diffN_MBTRIQ5y(x, y);
813 }
815}
#define diffN_MBTRIQ0x(x, y)
Definition fem_tools.h:117
#define diffN_MBTRIQ1x(x, y)
Definition fem_tools.h:119
#define diffN_MBTRIQ5x(x, y)
Definition fem_tools.h:127
#define diffN_MBTRIQ2x(x, y)
Definition fem_tools.h:121
#define diffN_MBTRIQ2y(x, y)
Definition fem_tools.h:122
#define diffN_MBTRIQ1y(x, y)
Definition fem_tools.h:120
#define diffN_MBTRIQ4y(x, y)
Definition fem_tools.h:126
#define diffN_MBTRIQ5y(x, y)
Definition fem_tools.h:128
#define diffN_MBTRIQ3x(x, y)
Definition fem_tools.h:123
#define diffN_MBTRIQ3y(x, y)
Definition fem_tools.h:124
#define diffN_MBTRIQ4x(x, y)
Definition fem_tools.h:125
#define diffN_MBTRIQ0y(x, y)
Definition fem_tools.h:118

◆ ShapeFaceBaseMBTRI()

PetscErrorCode ShapeFaceBaseMBTRI ( double * diffN,
const double * coords,
double * normal,
double * s1,
double * s2 )

Definition at line 204 of file fem_tools.c.

205 {
207
208 double diffX_ksi[3];
209 double diffX_eta[3];
210 int ii = 0;
211 for (; ii < 3; ii++) {
212 diffX_ksi[ii] = cblas_ddot(3, &coords[ii], 3, &diffN[0], 2);
213 diffX_eta[ii] = cblas_ddot(3, &coords[ii], 3, &diffN[1], 2);
214 }
215 if (s1 != NULL) {
216 cblas_dcopy(3, diffX_ksi, 1, s1, 1);
217 }
218 if (s2 != NULL) {
219 cblas_dcopy(3, diffX_eta, 1, s2, 1);
220 }
221 double Spin_diffX_ksi[9];
222 ierr = Spin(Spin_diffX_ksi, diffX_ksi);
223 CHKERRQ(ierr);
224 cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1., Spin_diffX_ksi, 3,
225 diffX_eta, 1, 0., normal, 1);
227}
PetscErrorCode Spin(double *spinOmega, double *vecOmega)
calculate spin matrix from vector
Definition fem_tools.c:546

◆ ShapeFaceDiffNormalMBTRI()

PetscErrorCode ShapeFaceDiffNormalMBTRI ( double * diffN,
const double * coords,
double * diff_normal )

calculate derivative of normal in respect to nodal positions

Definition at line 237 of file fem_tools.c.

238 {
240 // N = Spin(dX/dksi)*dX/deta = -Spin(dX/deta)*dX/dksi
241
242 double diffX_ksi[3];
243 double diffX_eta[3];
244 int ii = 0;
245 for (; ii < 3; ii++) {
246 diffX_ksi[ii] = cblas_ddot(3, &coords[ii], 3, &diffN[0], 2);
247 diffX_eta[ii] = cblas_ddot(3, &coords[ii], 3, &diffN[1], 2);
248 }
249 double Spin_diffX_ksi[9];
250 ierr = Spin(Spin_diffX_ksi, diffX_ksi);
251 CHKERRQ(ierr);
252 double Spin_diffX_eta[9];
253 ierr = Spin(Spin_diffX_eta, diffX_eta);
254 CHKERRQ(ierr);
255 double B_ksi[3 * 9];
256 bzero(B_ksi, 3 * 9 * sizeof(double));
257 double B_eta[3 * 9];
258 bzero(B_eta, 3 * 9 * sizeof(double));
259 // B_ksi[] = [
260 // diffN[2*0+0], 0, 0, diffN[2*1+0], 0, 0, diffN[2*2+0], 0,
261 // 0
262 // 0, diffN[2*0+0], 0, 0, diffN[2*1+0], 0, 0, diffN[2*2+0],
263 // 0
264 // 0, 0, diffM[2*0+0], 0, 0, diffN[2*1+0], 0, 0,
265 // diffN[2*2+0]
266 //]
267 // B_eta[] = [
268 // diffN[2*0+1], 0, 0, diffN[2*1+1], 0, 0, diffN[2*2+1], 0,
269 // 0
270 // 0, diffN[2*0+1], 0, 0, diffN[2*1+1], 0, 0, diffN[2*2+1],
271 // 0
272 // 0, 0, diffM[2*0+1], 0, 0, diffN[2*1+1], 0, 0,
273 // diffN[2*2+1]
274 //]
275 ii = 0;
276 for (; ii < 3; ii++) {
277 cblas_dcopy(3, &diffN[0], 2, &B_ksi[ii * 9 + ii], 3);
278 cblas_dcopy(3, &diffN[1], 2, &B_eta[ii * 9 + ii], 3);
279 }
280 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 3, 9, 3, +1.,
281 Spin_diffX_ksi, 3, B_eta, 9, 0., diff_normal, 9);
282 cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 3, 9, 3, -1.,
283 Spin_diffX_eta, 3, B_ksi, 9, 1., diff_normal, 9);
285}

◆ ShapeFaceNormalMBTRI()

PetscErrorCode ShapeFaceNormalMBTRI ( double * diffN,
const double * coords,
double * normal )

calculate face normal

Parameters
diffNderivatives of shape functions
coordsis position of the nodes
normalvector

Definition at line 229 of file fem_tools.c.

230 {
232 ierr = ShapeFaceBaseMBTRI(diffN, coords, normal, NULL, NULL);
233 CHKERRQ(ierr);
235}
PetscErrorCode ShapeFaceBaseMBTRI(double *diffN, const double *coords, double *normal, double *s1, double *s2)
Definition fem_tools.c:204

◆ ShapeFaceNormalMBTRI_complex()

PetscErrorCode ShapeFaceNormalMBTRI_complex ( double * diffN,
__CLPK_doublecomplex * xcoords,
__CLPK_doublecomplex * xnormal )

Definition at line 452 of file fem_tools.c.

454 {
456 double complex diffX_x, diffX_y, diffX_z;
457 double complex diffY_x, diffY_y, diffY_z;
458 diffX_x = diffX_y = diffX_z = 0.;
459 diffY_x = diffY_y = diffY_z = 0.;
460 int ii;
461 for (ii = 0; ii < 3; ii++) {
462 diffX_x +=
463 (xcoords[3 * ii + 0].r + I * xcoords[3 * ii + 0].i) * diffN[2 * ii + 0];
464 diffX_y +=
465 (xcoords[3 * ii + 1].r + I * xcoords[3 * ii + 1].i) * diffN[2 * ii + 0];
466 diffX_z +=
467 (xcoords[3 * ii + 2].r + I * xcoords[3 * ii + 2].i) * diffN[2 * ii + 0];
468 diffY_x +=
469 (xcoords[3 * ii + 0].r + I * xcoords[3 * ii + 0].i) * diffN[2 * ii + 1];
470 diffY_y +=
471 (xcoords[3 * ii + 1].r + I * xcoords[3 * ii + 1].i) * diffN[2 * ii + 1];
472 diffY_z +=
473 (xcoords[3 * ii + 2].r + I * xcoords[3 * ii + 2].i) * diffN[2 * ii + 1];
474 }
475 double complex tmp;
476 tmp = diffX_y * diffY_z - diffX_z * diffY_y;
477 xnormal[0].r = creal(tmp);
478 xnormal[0].i = cimag(tmp);
479 tmp = diffX_z * diffY_x - diffX_x * diffY_z;
480 xnormal[1].r = creal(tmp);
481 xnormal[1].i = cimag(tmp);
482 tmp = diffX_x * diffY_y - diffX_y * diffY_x;
483 xnormal[2].r = creal(tmp);
484 xnormal[2].i = cimag(tmp);
486}

◆ ShapeInvJacVolume()

PetscErrorCode ShapeInvJacVolume ( double * jac)

Definition at line 39 of file fem_tools.c.

39 {
41 __CLPK_integer ipiv[4];
42 __CLPK_doublereal work[3];
43 __CLPK_integer lwork = 3;
44 __CLPK_integer info;
45 info = lapack_dgetrf(3, 3, jac, 3, ipiv);
46 if (info != 0)
47 SETERRQ1(PETSC_COMM_SELF, 1, "info = %d", info);
48 info = lapack_dgetri(3, jac, 3, ipiv, work, lwork);
49 if (info != 0)
50 SETERRQ1(PETSC_COMM_SELF, 1, "info = %d", info);
52}
static __CLPK_integer lapack_dgetri(__CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *work, __CLPK_integer lwork)

◆ ShapeJacMBTET()

PetscErrorCode ShapeJacMBTET ( double * diffN,
const double * coords,
double * jac )

calculate jacobian

Definition at line 288 of file fem_tools.c.

288 {
290 int ii, jj, kk;
291 bzero(jac, sizeof(double) * 9);
292 for (ii = 0; ii < 4; ii++) // shape func.
293 for (jj = 0; jj < 3; jj++) // space
294 for (kk = 0; kk < 3; kk++) // derivative of shape func.
295 jac[jj * 3 + kk] += diffN[ii * 3 + kk] * coords[ii * 3 + jj];
297}

◆ ShapeJacMBTETQ()

PetscErrorCode ShapeJacMBTETQ ( const double * diffN,
const double * coords,
double * Jac )

Definition at line 967 of file fem_tools.c.

968 {
970 int ii, jj, kk;
971 bzero(Jac, sizeof(double) * 9);
972 for (ii = 0; ii < 10; ii++) // shape func.
973 for (jj = 0; jj < 3; jj++) // space
974 for (kk = 0; kk < 3; kk++) // derivative of shape func.
975 Jac[jj * 3 + kk] += diffN[ii * 3 + kk] * coords[ii * 3 + jj];
977}

◆ ShapeMBEDGE()

PetscErrorCode ShapeMBEDGE ( double * N,
const double * G_X,
int DIM )

Definition at line 761 of file fem_tools.c.

761 {
763 int ii = 0;
764 for (; ii < DIM; ii++) {
765 double x = G_X[ii];
766 N[2 * ii + 0] = N_MBEDGE0(x);
767 N[2 * ii + 1] = N_MBEDGE1(x);
768 }
770}
#define N_MBEDGE0(x)
edge shape function
Definition fem_tools.h:105
#define N_MBEDGE1(x)
edge shape function
Definition fem_tools.h:106
const int N
Definition speed_test.cpp:3

◆ ShapeMBTET()

PetscErrorCode ShapeMBTET ( double * N,
const double * G_X,
const double * G_Y,
const double * G_Z,
int DIM )

calculate shape functions

Definition at line 306 of file fem_tools.c.

307 {
309 int ii = 0;
310 for (; ii < DIM; ii++) {
311 double x = G_X[ii], y = G_Y[ii], z = G_Z[ii];
312 N[4 * ii + 0] = N_MBTET0(x, y, z);
313 N[4 * ii + 1] = N_MBTET1(x, y, z);
314 N[4 * ii + 2] = N_MBTET2(x, y, z);
315 N[4 * ii + 3] = N_MBTET3(x, y, z);
316 }
318}
#define N_MBTET0(x, y, z)
tetrahedral shape function
Definition fem_tools.h:28
#define N_MBTET1(x, y, z)
tetrahedral shape function
Definition fem_tools.h:29
#define N_MBTET2(x, y, z)
tetrahedral shape function
Definition fem_tools.h:30
#define N_MBTET3(x, y, z)
tetrahedral shape function
Definition fem_tools.h:31

◆ ShapeMBTET_inverse()

PetscErrorCode ShapeMBTET_inverse ( double * N,
double * diffN,
const double * elem_coords,
const double * glob_coords,
double * loc_coords )

calculate local coordinates for given global coordinates

new version for multiple points need to be implemented

Definition at line 335 of file fem_tools.c.

338 {
340 double A[3 * 3];
341 int IPIV[3];
342 // COL MAJOR
343 // X
344 A[0 + 3 * 0] =
345 cblas_ddot(4, &diffN[0 * 3 + 0], 3, &elem_coords[0 * 3 + 0], 3);
346 A[0 + 3 * 1] =
347 cblas_ddot(4, &diffN[0 * 3 + 1], 3, &elem_coords[0 * 3 + 0], 3);
348 A[0 + 3 * 2] =
349 cblas_ddot(4, &diffN[0 * 3 + 2], 3, &elem_coords[0 * 3 + 0], 3);
350 loc_coords[0] =
351 glob_coords[0] - cblas_ddot(4, &N[0], 1, &elem_coords[0 * 3 + 0], 3);
352 // printf("A\n[ %3.2f %3.2f %3.2f ] %3.2f \n",A[0*3],A[1*3],A[2*3],R[0]);
353 // Y
354 A[1 + 3 * 0] =
355 cblas_ddot(4, &diffN[0 * 3 + 0], 3, &elem_coords[0 * 3 + 1], 3);
356 A[1 + 3 * 1] =
357 cblas_ddot(4, &diffN[0 * 3 + 1], 3, &elem_coords[0 * 3 + 1], 3);
358 A[1 + 3 * 2] =
359 cblas_ddot(4, &diffN[0 * 3 + 2], 3, &elem_coords[0 * 3 + 1], 3);
360 loc_coords[1] =
361 glob_coords[1] - cblas_ddot(4, &N[0], 1, &elem_coords[0 * 3 + 1], 3);
362 // printf("[ %3.2f %3.2f %3.2f ] %3.2f \n",A[1+3*0],A[1+3*1],A[1+3*2],R[1]);
363 // Z
364 A[2 + 3 * 0] =
365 cblas_ddot(4, &diffN[0 * 3 + 0], 3, &elem_coords[0 * 3 + 2], 3);
366 A[2 + 3 * 1] =
367 cblas_ddot(4, &diffN[0 * 3 + 1], 3, &elem_coords[0 * 3 + 2], 3);
368 A[2 + 3 * 2] =
369 cblas_ddot(4, &diffN[0 * 3 + 2], 3, &elem_coords[0 * 3 + 2], 3);
370 loc_coords[2] =
371 glob_coords[2] - cblas_ddot(4, &N[0], 1, &elem_coords[0 * 3 + 2], 3);
372 // printf("[ %3.2f %3.2f %3.2f ] %3.2f \n",A[2+3*0],A[2+3*1],A[2+3*2],R[1]);
373 int info =
374 lapack_dgesv(3, 1, &A[0], 3, (__CLPK_integer *)IPIV, loc_coords, 3);
375 if (info != 0)
376 SETERRQ1(PETSC_COMM_SELF, 1, "info == %d", info);
378}
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)

◆ ShapeMBTETQ()

PetscErrorCode ShapeMBTETQ ( double * N,
const double x,
const double y,
const double z )

Definition at line 858 of file fem_tools.c.

859 {
861 N[0] = N_MBTETQ0(x, y, z);
862 N[1] = N_MBTETQ1(x, y, z);
863 N[2] = N_MBTETQ2(x, y, z);
864 N[3] = N_MBTETQ3(x, y, z);
865 N[4] = N_MBTETQ4(x, y, z);
866 N[5] = N_MBTETQ5(x, y, z);
867 N[6] = N_MBTETQ6(x, y, z);
868 N[7] = N_MBTETQ7(x, y, z);
869 N[8] = N_MBTETQ8(x, y, z);
870 N[9] = N_MBTETQ9(x, y, z);
872}
#define N_MBTETQ2(x, y, z)
Definition fem_tools.c:820
#define N_MBTETQ9(x, y, z)
Definition fem_tools.c:827
#define N_MBTETQ4(x, y, z)
Definition fem_tools.c:822
#define N_MBTETQ8(x, y, z)
Definition fem_tools.c:826
#define N_MBTETQ6(x, y, z)
Definition fem_tools.c:824
#define N_MBTETQ7(x, y, z)
Definition fem_tools.c:825
#define N_MBTETQ5(x, y, z)
Definition fem_tools.c:823
#define N_MBTETQ1(x, y, z)
Definition fem_tools.c:819
#define N_MBTETQ3(x, y, z)
Definition fem_tools.c:821
#define N_MBTETQ0(x, y, z)
Definition fem_tools.c:818

◆ ShapeMBTETQ_detJac_at_Gauss_Points()

PetscErrorCode ShapeMBTETQ_detJac_at_Gauss_Points ( double * detJac_at_Gauss_Points,
const double * diffN,
const double * coords,
int G_DIM )

Definition at line 979 of file fem_tools.c.

981 {
983
984 double Jac[9];
985 int ii = 0;
986 for (; ii < G_DIM; ii++) {
987 ierr = ShapeJacMBTETQ(&diffN[30 * ii], coords, Jac);
988 CHKERRQ(ierr);
989 detJac_at_Gauss_Points[ii] = ShapeDetJacVolume(Jac);
990 }
992}
double ShapeDetJacVolume(double *jac)
determined of jacobian
Definition fem_tools.c:22
PetscErrorCode ShapeJacMBTETQ(const double *diffN, const double *coords, double *Jac)
Definition fem_tools.c:967

◆ ShapeMBTETQ_GAUSS()

PetscErrorCode ShapeMBTETQ_GAUSS ( double * N,
const double * X,
const double * Y,
const double * Z,
const int G_DIM )

Definition at line 908 of file fem_tools.c.

909 {
911 int ii = 0;
912 for (; ii < G_DIM; ii++) {
913 double x = X[ii], y = Y[ii], z = Z[ii];
914 N[10 * ii + 0] = N_MBTETQ0(x, y, z);
915 N[10 * ii + 1] = N_MBTETQ1(x, y, z);
916 N[10 * ii + 2] = N_MBTETQ2(x, y, z);
917 N[10 * ii + 3] = N_MBTETQ3(x, y, z);
918 N[10 * ii + 4] = N_MBTETQ4(x, y, z);
919 N[10 * ii + 5] = N_MBTETQ5(x, y, z);
920 N[10 * ii + 6] = N_MBTETQ6(x, y, z);
921 N[10 * ii + 7] = N_MBTETQ7(x, y, z);
922 N[10 * ii + 8] = N_MBTETQ8(x, y, z);
923 N[10 * ii + 9] = N_MBTETQ9(x, y, z);
924 }
926}

◆ ShapeMBTETQ_inverse()

PetscErrorCode ShapeMBTETQ_inverse ( double * N,
double * diffN,
const double * elem_coords,
const double * glob_coords,
double * loc_coords,
const double eps )

Definition at line 1007 of file fem_tools.c.

1010 {
1012 double A[3 * 3];
1013 double R[3];
1014 int IPIV[3];
1015 float NORM_dR = 1000.;
1016 float NORM_R0;
1017 ShapeMBTETQ(N, 0.1, 0.1, 0.1);
1018 ShapeDiffMBTETQ(diffN, 0.1, 0.1, 0.1);
1019 R[0] = glob_coords[0] - cblas_ddot(10, &N[0], 1, &elem_coords[0], 3);
1020 R[1] = glob_coords[1] - cblas_ddot(10, &N[0], 1, &elem_coords[1], 3);
1021 R[2] = glob_coords[2] - cblas_ddot(10, &N[0], 1, &elem_coords[2], 3);
1022 NORM_R0 = cblas_dnrm2(3, &R[0], 1);
1023 while ((NORM_dR / NORM_R0) > eps) {
1024 // COL MAJOR
1025 // X
1026 A[0 + 3 * 0] = cblas_ddot(10, &diffN[0 * 3 + 0], 3, &elem_coords[0], 3);
1027 A[0 + 3 * 1] = cblas_ddot(10, &diffN[0 * 3 + 1], 3, &elem_coords[0], 3);
1028 A[0 + 3 * 2] = cblas_ddot(10, &diffN[0 * 3 + 2], 3, &elem_coords[0], 3);
1029 R[0] = glob_coords[0] - cblas_ddot(10, &N[0], 1, &elem_coords[0], 3);
1030 // Y
1031 A[1 + 3 * 0] = cblas_ddot(10, &diffN[0 * 3 + 0], 3, &elem_coords[1], 3);
1032 A[1 + 3 * 1] = cblas_ddot(10, &diffN[0 * 3 + 1], 3, &elem_coords[1], 3);
1033 A[1 + 3 * 2] = cblas_ddot(10, &diffN[0 * 3 + 2], 3, &elem_coords[1], 3);
1034 R[1] = glob_coords[1] - cblas_ddot(10, &N[0], 1, &elem_coords[1], 3);
1035 // Z
1036 A[2 + 3 * 0] =
1037 cblas_ddot(10, &diffN[0 * 3 + 0], 3, &elem_coords[0 * 3 + 2], 3);
1038 A[2 + 3 * 1] =
1039 cblas_ddot(10, &diffN[0 * 3 + 1], 3, &elem_coords[0 * 3 + 2], 3);
1040 A[2 + 3 * 2] =
1041 cblas_ddot(10, &diffN[0 * 3 + 2], 3, &elem_coords[0 * 3 + 2], 3);
1042 R[2] = glob_coords[2] - cblas_ddot(10, &N[0], 1, &elem_coords[2], 3);
1043 int info = lapack_dgesv(3, 1, &A[0], 3, (__CLPK_integer *)IPIV, R, 3);
1044 assert(info == 0);
1045 NOT_USED(info);
1046 cblas_daxpy(3, 1., R, 1, loc_coords, 1);
1047 NORM_dR = cblas_dnrm2(3, &R[0], 1);
1048 ShapeMBTETQ(N, loc_coords[0], loc_coords[1], loc_coords[2]);
1049 ShapeDiffMBTETQ(diffN, loc_coords[0], loc_coords[1], loc_coords[2]);
1050 }
1052}
static const double eps
#define NOT_USED(x)
PetscErrorCode ShapeDiffMBTETQ(double *diffN, const double x, const double y, const double z)
Definition fem_tools.c:873
PetscErrorCode ShapeMBTETQ(double *N, const double x, const double y, const double z)
Definition fem_tools.c:858
@ R

◆ ShapeMBTRI()

PetscErrorCode ShapeMBTRI ( double * N,
const double * X,
const double * Y,
const int G_DIM )

calculate shape functions on triangle

Parameters
Nshape function array
Xarray of Gauss X coordinates
Yarray of Gauss Y coordinates
G_DIMnumber of Gauss points

Definition at line 182 of file fem_tools.c.

183 {
185 int ii = 0;
186 for (; ii < G_DIM; ii++) {
187 double x = X[ii], y = Y[ii];
188 N[3 * ii + 0] = N_MBTRI0(x, y);
189 N[3 * ii + 1] = N_MBTRI1(x, y);
190 N[3 * ii + 2] = N_MBTRI2(x, y);
191 }
193}
#define N_MBTRI1(x, y)
triangle shape function
Definition fem_tools.h:47
#define N_MBTRI0(x, y)
triangle shape function
Definition fem_tools.h:46
#define N_MBTRI2(x, y)
triangle shape function
Definition fem_tools.h:48

◆ ShapeMBTRI_inverse()

PetscErrorCode ShapeMBTRI_inverse ( double * N,
double * diffN,
const double * elem_coords,
const double * glob_coords,
double * loc_coords )

calculate local coordinates of triangle element for given global coordinates in 2D (Assume e.g. z=0)

\[ \left[\begin{array}{cc} \frac{\partial N_{1}}{\partial\xi}x_{N_{1}}+\frac{\partial N_{2}}{\partial\xi}x_{N_{2}}+\frac{\partial N_{3}}{\partial\xi}x_{N_{3}} & \frac{\partial N_{1}}{\partial\eta}x_{N_{1}}+\frac{\partial N_{2}}{\partial\eta}x_{N_{2}}+\frac{\partial N_{3}}{\partial\eta}x_{N_{3}}\\ \frac{\partial N_{1}}{\partial\xi}y_{N_{1}}+\frac{\partial N_{2}}{\partial\xi}y_{N_{2}}+\frac{\partial N_{3}}{\partial\xi}y_{N_{3}} & \frac{\partial N_{1}}{\partial\eta}y_{N_{1}}+\frac{\partial N_{2}}{\partial\eta}y_{N_{2}}+\frac{\partial N_{3}}{\partial\eta}y_{N_{3}} \end{array}\right]\left\{ \begin{array}{c} \xi\\ \eta \end{array}\right\} =\left\{ \begin{array}{c} x_{gp}-\left(N_{1}x_{N_{1}}+N_{2}x_{N_{2}}+N_{3}x_{N_{3}}\right)\\ y_{gp}-\left(N_{1}y_{N_{1}}+N_{2}y_{N_{2}}+N_{3}y_{N_{3}}\right) \end{array}\right\} \]

/

Parameters
Nshape function array /
diffNarray of shape function derivative w.r.t to local coordinates /
elem_coordsglobal coordinates of element /
glob_coordsglobal coordinates of required point /
loc_coordslocal coordinates of required point

Definition at line 380 of file fem_tools.c.

383 {
385 double A[2 * 2];
386
387 // 1st and 2nd element of matrix A
388 A[0] = cblas_ddot(3, &diffN[0], 2, &elem_coords[0], 2); // dot product
389 A[1] = cblas_ddot(3, &diffN[1], 2, &elem_coords[0], 2);
390 loc_coords[0] = glob_coords[0] - cblas_ddot(3, &N[0], 1, &elem_coords[0], 2);
391
392 // 3rd and 4th element of matrix A
393 A[2] = cblas_ddot(3, &diffN[0], 2, &elem_coords[1], 2);
394 A[3] = cblas_ddot(3, &diffN[1], 2, &elem_coords[1], 2);
395 loc_coords[1] = glob_coords[1] - cblas_ddot(3, &N[0], 1, &elem_coords[1], 2);
396
397 // calculate directly the solution (as the size of matrix is only 2x2)
398 double invA[2 * 2], detA;
399 detA = A[0] * A[3] - A[1] * A[2];
400 detA = 1.0 / detA;
401 invA[0] = A[3] * detA;
402 invA[1] = -1.0 * A[1] * detA;
403 invA[2] = -1.0 * A[2] * detA;
404 invA[3] = A[0] * detA;
405
406 double loc_coords_new[2];
407 loc_coords_new[0] = invA[0] * loc_coords[0] + invA[1] * loc_coords[1];
408 loc_coords_new[1] = invA[2] * loc_coords[0] + invA[3] * loc_coords[1];
409
410 loc_coords[0] = loc_coords_new[0];
411 loc_coords[1] = loc_coords_new[1];
413}

◆ ShapeMBTRIQ()

PetscErrorCode ShapeMBTRIQ ( double * N,
const double * X,
const double * Y,
const int G_DIM )

Definition at line 780 of file fem_tools.c.

781 {
783 int ii = 0;
784 for (; ii < G_DIM; ii++) {
785 double x = X[ii], y = Y[ii];
786 N[6 * ii + 0] = N_MBTRIQ0(x, y);
787 N[6 * ii + 1] = N_MBTRIQ1(x, y);
788 N[6 * ii + 2] = N_MBTRIQ2(x, y);
789 N[6 * ii + 3] = N_MBTRIQ3(x, y);
790 N[6 * ii + 4] = N_MBTRIQ4(x, y);
791 N[6 * ii + 5] = N_MBTRIQ5(x, y);
792 }
794}
#define N_MBTRIQ2(x, y)
Definition fem_tools.h:113
#define N_MBTRIQ3(x, y)
Definition fem_tools.h:114
#define N_MBTRIQ0(x, y)
Definition fem_tools.h:111
#define N_MBTRIQ4(x, y)
Definition fem_tools.h:115
#define N_MBTRIQ5(x, y)
Definition fem_tools.h:116
#define N_MBTRIQ1(x, y)
Definition fem_tools.h:112

◆ ShapeVolumeMBTET()

double ShapeVolumeMBTET ( double * diffN,
const double * coords )

calculate TET volume

Definition at line 298 of file fem_tools.c.

298 {
299 double Jac[9];
300 ShapeJacMBTET(diffN, coords, Jac);
301 double detJac = ShapeDetJacVolume(Jac);
302 // printf("detJac = +%6.4e\n",detJac);
303 // print_mat(Jac,3,3);
304 return detJac * G_TET_W1[0] / 6.;
305}
PetscErrorCode ShapeJacMBTET(double *diffN, const double *coords, double *jac)
calculate jacobian
Definition fem_tools.c:288
static const double G_TET_W1[]
Definition fem_tools.h:1115

◆ ShapeVolumeMBTETQ()

double ShapeVolumeMBTETQ ( const double * diffN,
const double * coords,
int G_DIM,
double * G_TET_W )

Definition at line 993 of file fem_tools.c.

994 {
995
996 int ii = 0;
997 double vol = 0;
998 double detJac_at_Gauss_Points[G_DIM];
999 ierr = ShapeMBTETQ_detJac_at_Gauss_Points(detJac_at_Gauss_Points, diffN,
1000 coords, G_DIM);
1001 CHKERRQ(ierr);
1002 for (; ii < G_DIM; ii++) {
1003 vol += G_TET_W[ii] * (detJac_at_Gauss_Points[ii]) / 6;
1004 }
1005 return vol;
1006}
PetscErrorCode ShapeMBTETQ_detJac_at_Gauss_Points(double *detJac_at_Gauss_Points, const double *diffN, const double *coords, int G_DIM)
Definition fem_tools.c:979

◆ Spin()

PetscErrorCode Spin ( double * spinOmega,
double * vecOmega )

calculate spin matrix from vector

Definition at line 546 of file fem_tools.c.

546 {
548 bzero(spinOmega, 9 * sizeof(double));
549 spinOmega[0 * 3 + 1] = -vecOmega[2];
550 spinOmega[0 * 3 + 2] = +vecOmega[1];
551 spinOmega[1 * 3 + 0] = +vecOmega[2];
552 spinOmega[1 * 3 + 2] = -vecOmega[0];
553 spinOmega[2 * 3 + 0] = -vecOmega[1];
554 spinOmega[2 * 3 + 1] = +vecOmega[0];
556}

Variable Documentation

◆ ierr

PetscErrorCode ierr
static

Definition at line 20 of file fem_tools.c.