v0.9.0
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 More...
 
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. More...
 
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. More...
 
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. More...
 
PetscErrorCode ShapeMBTRI (double *N, const double *X, const double *Y, const int G_DIM)
 calculate shape functions on triangle More...
 
PetscErrorCode ShapeDiffMBTRI (double *diffN)
 calculate derivatives of shape functions More...
 
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 More...
 
PetscErrorCode ShapeJacMBTET (double *diffN, const double *coords, double *jac)
 calculate jacobian More...
 
double ShapeVolumeMBTET (double *diffN, const double *coords)
 calculate TET volume More...
 
PetscErrorCode ShapeMBTET (double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
 calculate shape functions More...
 
PetscErrorCode ShapeDiffMBTET (double *diffN)
 calculate derivatives of shape functions More...
 
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 More...
 
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\} \]

More...
 
PetscErrorCode ShapeDiffMBTETinvJ (double *diffN, double *invJac, double *diffNinvJac)
 calculate shape functions derivatives in space More...
 
PetscErrorCode GradientOfDeformation (double *diffN, double *dofs, double *F)
 calculate gradient of deformation More...
 
void ShapeDiffMBTETinvJ_complex (double *diffN, __CLPK_doublecomplex *invJac, __CLPK_doublecomplex *diffNinvJac, const 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 More...
 
PetscErrorCode make_complex_matrix (double *reA, double *imA, __CLPK_doublecomplex *xA)
 Compose complex matrix (3x3) from two real matrices. More...
 
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. More...
 
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,
 
)    (-3. + 4. * x + 4. * y + 4. * z)

Definition at line 840 of file fem_tools.c.

◆ diffN_MBTETQ0y

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

Definition at line 841 of file fem_tools.c.

◆ diffN_MBTETQ0z

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

Definition at line 842 of file fem_tools.c.

◆ diffN_MBTETQ1x

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

Definition at line 843 of file fem_tools.c.

◆ diffN_MBTETQ1y

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

Definition at line 844 of file fem_tools.c.

◆ diffN_MBTETQ1z

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

Definition at line 845 of file fem_tools.c.

◆ diffN_MBTETQ2x

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

Definition at line 846 of file fem_tools.c.

◆ diffN_MBTETQ2y

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

Definition at line 847 of file fem_tools.c.

◆ diffN_MBTETQ2z

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

Definition at line 848 of file fem_tools.c.

◆ diffN_MBTETQ3x

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

Definition at line 849 of file fem_tools.c.

◆ diffN_MBTETQ3y

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

Definition at line 850 of file fem_tools.c.

◆ diffN_MBTETQ3z

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

Definition at line 851 of file fem_tools.c.

◆ diffN_MBTETQ4x

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

Definition at line 852 of file fem_tools.c.

◆ diffN_MBTETQ4y

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

Definition at line 853 of file fem_tools.c.

◆ diffN_MBTETQ4z

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

Definition at line 854 of file fem_tools.c.

◆ diffN_MBTETQ5x

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

Definition at line 855 of file fem_tools.c.

◆ diffN_MBTETQ5y

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

Definition at line 856 of file fem_tools.c.

◆ diffN_MBTETQ5z

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

Definition at line 857 of file fem_tools.c.

◆ diffN_MBTETQ6x

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

Definition at line 858 of file fem_tools.c.

◆ diffN_MBTETQ6y

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

Definition at line 859 of file fem_tools.c.

◆ diffN_MBTETQ6z

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

Definition at line 860 of file fem_tools.c.

◆ diffN_MBTETQ7x

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

Definition at line 861 of file fem_tools.c.

◆ diffN_MBTETQ7y

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

Definition at line 862 of file fem_tools.c.

◆ diffN_MBTETQ7z

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

Definition at line 863 of file fem_tools.c.

◆ diffN_MBTETQ8x

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

Definition at line 864 of file fem_tools.c.

◆ diffN_MBTETQ8y

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

Definition at line 865 of file fem_tools.c.

◆ diffN_MBTETQ8z

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

Definition at line 866 of file fem_tools.c.

◆ diffN_MBTETQ9x

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

Definition at line 867 of file fem_tools.c.

◆ diffN_MBTETQ9y

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

Definition at line 868 of file fem_tools.c.

◆ diffN_MBTETQ9z

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

Definition at line 869 of file fem_tools.c.

◆ N_MBTETQ0

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

Definition at line 830 of file fem_tools.c.

◆ N_MBTETQ1

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

Definition at line 831 of file fem_tools.c.

◆ N_MBTETQ2

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

Definition at line 832 of file fem_tools.c.

◆ N_MBTETQ3

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

Definition at line 833 of file fem_tools.c.

◆ N_MBTETQ4

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

Definition at line 834 of file fem_tools.c.

◆ N_MBTETQ5

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

Definition at line 835 of file fem_tools.c.

◆ N_MBTETQ6

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

Definition at line 836 of file fem_tools.c.

◆ N_MBTETQ7

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

Definition at line 837 of file fem_tools.c.

◆ N_MBTETQ8

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

Definition at line 838 of file fem_tools.c.

◆ N_MBTETQ9

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

Definition at line 839 of file fem_tools.c.

Function Documentation

◆ Base_scale()

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

Definition at line 753 of file fem_tools.c.

755  {
757  complex double xnrm2_normal = csqrt(cpow(xnormal[0].r + I * xnormal[0].i, 2) +
758  cpow(xnormal[1].r + I * xnormal[1].i, 2) +
759  cpow(xnormal[2].r + I * xnormal[2].i, 2));
760  int dd = 0;
761  for (; dd < 3; dd++) {
762  complex double s1 = (xs1[dd].r + I * xs1[dd].i) * xnrm2_normal;
763  complex double s2 = (xs2[dd].r + I * xs2[dd].i) * xnrm2_normal;
764  xs1[dd].r = creal(s1);
765  xs1[dd].i = cimag(s1);
766  xs2[dd].r = creal(s2);
767  xs2[dd].i = cimag(s2);
768  }
770 }
#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
__CLPK_doublereal i
Definition: lapack_wrap.h:47
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
__CLPK_doublereal r
Definition: lapack_wrap.h:47

◆ DeterminantComplexGradient()

PetscErrorCode DeterminantComplexGradient ( __CLPK_doublecomplex xF,
__CLPK_doublecomplex det_xF 
)

Definition at line 538 of file fem_tools.c.

539  {
541  __CLPK_integer IPIV[4];
542  if (lapack_zgetrf(3, 3, xF, 3, IPIV) != 0) {
543  SETERRQ(PETSC_COMM_SELF, 1, "lapack_zgetrf(3,3,xF,3,IPIV) != 0");
544  }
545  double complex det = 1;
546  int i = 0, j = 0;
547  for (; i < 3; i++) {
548  det *= xF[3 * i + i].r + I * xF[3 * i + i].i;
549  if (IPIV[i] != i + 1)
550  j++;
551  }
552  if ((j - (j / 2) * 2) != 0)
553  det = -det;
554  (*det_xF).r = creal(det);
555  (*det_xF).i = cimag(det);
557 }
#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
__CLPK_doublereal i
Definition: lapack_wrap.h:47
__CLPK_doublereal r
Definition: lapack_wrap.h:47
static __CLPK_integer lapack_zgetrf(__CLPK_integer m, __CLPK_integer n, __CLPK_doublecomplex *a, __CLPK_integer lda, __CLPK_integer *ipiv)
Definition: lapack_wrap.h:335
long int __CLPK_integer
Definition: lapack_wrap.h:35

◆ GradientOfDeformation()

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

calculate gradient of deformation

Definition at line 437 of file fem_tools.c.

437  {
439  int col, row = 0;
440  for (; row < 3; row++)
441  for (col = 0; col < 3; col++) {
442  F[3 * row + col] = cblas_ddot(4, &diffN[col], 3, &dofs[row], 3);
443  }
445 }
#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

◆ 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 67 of file fem_tools.c.

69  {
71 
72  int dim_num = 1;
73  int point;
74  int point_num;
75  double *w;
76  double *x;
77 
78  // GM_RULE_SET determines the weights and abscissas
79  // pof a Grundmann-Moeller quadrature rule for
80  // the DIM_NUM dimensional simplex,
81  // using a rule of in index RULE,
82  // which will have degree of exactness 2*RULE+1.
83 
84  // printf ( " Here we use DIM_NUM = %d\n", dim_num );
85  // printf ( " RULE = %d\n", rule );
86  // printf ( " DEGREE = %d\n", 2 * rule + 1 );
87 
88  point_num = gm_rule_size(rule, dim_num);
89 
90  ierr = PetscMalloc(point_num * sizeof(double), &w);
91  CHKERRQ(ierr);
92  ierr = PetscMalloc(dim_num * point_num * sizeof(double), &x);
93  CHKERRQ(ierr);
94 
95  gm_rule_set(rule, dim_num, point_num, w, x);
96 
97  for (point = 0; point < point_num; point++) {
98  G_TRI_X[point] = x[0 + point * dim_num];
99  G_TRI_W[point] = w[point];
100  }
101 
102  ierr = PetscFree(w);
103  CHKERRQ(ierr);
104  ierr = PetscFree(x);
105  CHKERRQ(ierr);
106 
108 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
int gm_rule_size(int rule, int dim_num)
Definition: gm_rule.c:294
void gm_rule_set(int rule, int dim_num, int point_num, double w[], double x[])
Definition: gm_rule.c:152
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
CHKERRQ(ierr)
static PetscErrorCode ierr
Definition: fem_tools.c:32

◆ 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 110 of file fem_tools.c.

113  {
115 
116  int dim_num = 2;
117  int point;
118  int point_num;
119  double *w;
120  double *x;
121 
122  // GM_RULE_SET determines the weights and abscissas
123  // pof a Grundmann-Moeller quadrature rule for
124  // the DIM_NUM dimensional simplex,
125  // using a rule of in index RULE,
126  // which will have degree of exactness 2*RULE+1.
127 
128  // printf ( " Here we use DIM_NUM = %d\n", dim_num );
129  // printf ( " RULE = %d\n", rule );
130  // printf ( " DEGREE = %d\n", 2 * rule + 1 );
131 
132  point_num = gm_rule_size(rule, dim_num);
133 
134  ierr = PetscMalloc(point_num * sizeof(double), &w);
135  CHKERRQ(ierr);
136  ierr = PetscMalloc(dim_num * point_num * sizeof(double), &x);
137  CHKERRQ(ierr);
138 
139  gm_rule_set(rule, dim_num, point_num, w, x);
140 
141  for (point = 0; point < point_num; point++) {
142  G_TRI_X[point] = x[0 + point * dim_num];
143  G_TRI_Y[point] = x[1 + point * dim_num];
144  G_TRI_W[point] = w[point];
145  }
146 
147  ierr = PetscFree(w);
148  CHKERRQ(ierr);
149  ierr = PetscFree(x);
150  CHKERRQ(ierr);
151 
153 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
int gm_rule_size(int rule, int dim_num)
Definition: gm_rule.c:294
void gm_rule_set(int rule, int dim_num, int point_num, double w[], double x[])
Definition: gm_rule.c:152
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
CHKERRQ(ierr)
static PetscErrorCode ierr
Definition: fem_tools.c:32

◆ 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 155 of file fem_tools.c.

159  {
161 
162  int dim_num = 3;
163  int point;
164  int point_num;
165  double *w;
166  double *x;
167 
168  // printf ( " Here we use DIM_NUM = %d\n", dim_num );
169  // printf ( " RULE = %d\n", rule );
170  // printf ( " DEGREE = %d\n", 2 * rule + 1 );
171 
172  point_num = gm_rule_size(rule, dim_num);
173 
174  ierr = PetscMalloc(point_num * sizeof(double), &w);
175  CHKERRQ(ierr);
176  ierr = PetscMalloc(dim_num * point_num * sizeof(double), &x);
177  CHKERRQ(ierr);
178 
179  gm_rule_set(rule, dim_num, point_num, w, x);
180  for (point = 0; point < point_num; point++) {
181  G_TET_X[point] = x[0 + point * dim_num];
182  G_TET_Y[point] = x[1 + point * dim_num];
183  G_TET_Z[point] = x[2 + point * dim_num];
184  G_TET_W[point] = w[point];
185  }
186 
187  ierr = PetscFree(w);
188  CHKERRQ(ierr);
189  ierr = PetscFree(x);
190  CHKERRQ(ierr);
191 
193 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
int gm_rule_size(int rule, int dim_num)
Definition: gm_rule.c:294
void gm_rule_set(int rule, int dim_num, int point_num, double w[], double x[])
Definition: gm_rule.c:152
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
CHKERRQ(ierr)
static PetscErrorCode ierr
Definition: fem_tools.c:32

◆ InvertComplexGradient()

PetscErrorCode InvertComplexGradient ( __CLPK_doublecomplex xF)

Definition at line 511 of file fem_tools.c.

511  {
513  __CLPK_integer IPIV[4];
514  __CLPK_doublecomplex WORK[4];
515  __CLPK_integer LWORK = 4;
516  __CLPK_integer info;
517  info = lapack_zgetrf(3, 3, xF, 3, IPIV);
518  if (info != 0)
519  SETERRQ(PETSC_COMM_SELF, 1, "info == 0");
520  info = lapack_zgetri(3, xF, 3, IPIV, WORK, LWORK);
521  if (info != 0)
522  SETERRQ(PETSC_COMM_SELF, 1, "info == 0");
524 }
#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
static __CLPK_integer lapack_zgetrf(__CLPK_integer m, __CLPK_integer n, __CLPK_doublecomplex *a, __CLPK_integer lda, __CLPK_integer *ipiv)
Definition: lapack_wrap.h:335
static __CLPK_integer lapack_zgetri(__CLPK_integer n, __CLPK_doublecomplex *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublecomplex *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:345
long int __CLPK_integer
Definition: lapack_wrap.h:35

◆ InvertComplexSymmMatrix3by3()

PetscErrorCode InvertComplexSymmMatrix3by3 ( __CLPK_doublecomplex xC)

Definition at line 525 of file fem_tools.c.

525  {
527  __CLPK_integer info;
528  info = lapack_zpotrf('L', 3, xC, 3);
529  if (info == 0)
530  SETERRQ(PETSC_COMM_SELF, 1, "info == 0");
531  // assert(info == 0);
532  info = lapack_zpotri('L', 3, xC, 3);
533  if (info == 0)
534  SETERRQ(PETSC_COMM_SELF, 1, "info == 0");
535  // assert(info == 0);
537 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
static __CLPK_integer lapack_zpotri(char uplo, __CLPK_integer n, __CLPK_doublecomplex *a, __CLPK_integer lda)
Definition: lapack_wrap.h:353
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
static __CLPK_integer lapack_zpotrf(char uplo, __CLPK_integer n, __CLPK_doublecomplex *a, __CLPK_integer lda)
Definition: lapack_wrap.h:360
long int __CLPK_integer
Definition: lapack_wrap.h:35

◆ 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 569 of file fem_tools.c.

570  {
572  int ii = 0, jj;
573  for (; ii < 3; ii++) {
574  for (jj = 0; jj < 3; jj++) {
575  xA[3 * ii + jj].r = reA[3 * ii + jj];
576  xA[3 * ii + jj].i = imA[3 * ii + jj];
577  }
578  }
580 }
#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
__CLPK_doublereal i
Definition: lapack_wrap.h:47
__CLPK_doublereal r
Definition: lapack_wrap.h:47

◆ MakeComplexTensor()

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

Definition at line 499 of file fem_tools.c.

500  {
502  int ii = 0, jj;
503  for (; ii < 3; ii++) {
504  for (jj = 0; jj < 3; jj++) {
505  xA[3 * ii + jj].r = reA[3 * ii + jj];
506  xA[3 * ii + jj].i = imA[3 * ii + jj];
507  }
508  }
510 }
#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
__CLPK_doublereal i
Definition: lapack_wrap.h:47
__CLPK_doublereal r
Definition: lapack_wrap.h:47

◆ 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 581 of file fem_tools.c.

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

◆ ShapeDetJacVolume()

double ShapeDetJacVolume ( double jac)

determined of jacobian

Definition at line 34 of file fem_tools.c.

34  {
35  double det_jac;
36  __CLPK_integer ipiv[4];
37  __CLPK_integer info = lapack_dgetrf(3, 3, jac, 3, ipiv);
38  if (info != 0)
39  return -1;
40  int i = 0, j = 0;
41  det_jac = 1.;
42  for (; i < 3; i++) {
43  det_jac *= jac[3 * i + i];
44  if (ipiv[i] != i + 1)
45  j++;
46  }
47  if ((j - (j / 2) * 2) != 0)
48  det_jac = -det_jac;
49  return det_jac;
50 }
static __CLPK_integer lapack_dgetrf(__CLPK_integer m, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv)
Definition: lapack_wrap.h:169
long int __CLPK_integer
Definition: lapack_wrap.h:35

◆ ShapeDiffMBEDGE()

PetscErrorCode ShapeDiffMBEDGE ( double diffN)

Definition at line 783 of file fem_tools.c.

783  {
785  diffN[0] = diffN_MBEDGE0;
786  diffN[1] = diffN_MBEDGE1;
788 }
#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 diffN_MBEDGE1
derivative of edge shape function
Definition: fem_tools.h:82
#define diffN_MBEDGE0
derivative of edge shape function
Definition: fem_tools.h:81

◆ ShapeDiffMBTET()

PetscErrorCode ShapeDiffMBTET ( double diffN)

calculate derivatives of shape functions

Definition at line 331 of file fem_tools.c.

331  {
333  diffN[0] = diffN_MBTET0x;
334  diffN[1] = diffN_MBTET0y;
335  diffN[2] = diffN_MBTET0z;
336  diffN[3] = diffN_MBTET1x;
337  diffN[4] = diffN_MBTET1y;
338  diffN[5] = diffN_MBTET1z;
339  diffN[6] = diffN_MBTET2x;
340  diffN[7] = diffN_MBTET2y;
341  diffN[8] = diffN_MBTET2z;
342  diffN[9] = diffN_MBTET3x;
343  diffN[10] = diffN_MBTET3y;
344  diffN[11] = diffN_MBTET3z;
346 }
#define diffN_MBTET1z
derivative of tetrahedral shape function
Definition: fem_tools.h:45
#define diffN_MBTET0z
derivative of tetrahedral shape function
Definition: fem_tools.h:42
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define diffN_MBTET2x
derivative of tetrahedral shape function
Definition: fem_tools.h:46
#define diffN_MBTET3y
derivative of tetrahedral shape function
Definition: fem_tools.h:50
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
#define diffN_MBTET1x
derivative of tetrahedral shape function
Definition: fem_tools.h:43
#define diffN_MBTET3z
derivative of tetrahedral shape function
Definition: fem_tools.h:51
#define diffN_MBTET0y
derivative of tetrahedral shape function
Definition: fem_tools.h:41
#define diffN_MBTET2y
derivative of tetrahedral shape function
Definition: fem_tools.h:47
#define diffN_MBTET3x
derivative of tetrahedral shape function
Definition: fem_tools.h:49
#define diffN_MBTET1y
derivative of tetrahedral shape function
Definition: fem_tools.h:44
#define diffN_MBTET2z
derivative of tetrahedral shape function
Definition: fem_tools.h:48
#define diffN_MBTET0x
derivative of tetrahedral shape function
Definition: fem_tools.h:40

◆ ShapeDiffMBTETinvJ()

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

calculate shape functions derivatives in space

Definition at line 427 of file fem_tools.c.

428  {
430  int ii = 0;
431  for (; ii < 4; ii++) {
432  cblas_dgemv(CblasRowMajor, CblasTrans, 3, 3, 1., invJac, 3, &diffN[ii * 3],
433  1, 0., &diffNinvJac[ii * 3], 1);
434  }
436 }
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

◆ ShapeDiffMBTETinvJ_complex()

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

Definition at line 449 of file fem_tools.c.

451  {
452  __CLPK_doublecomplex tmp1 = {1., 0.}, tmp2 = {0., 0.};
453  int ii = 0, jj;
454  for (; ii < 4; ii++) {
455  __CLPK_doublecomplex tmp3[3];
456  for (jj = 0; jj < 3; jj++) {
457  tmp3[jj].r = diffN[ii * 3 + jj];
458  tmp3[jj].i = 0;
459  }
460  cblas_zgemv(CblasRowMajor, Trans, 3, 3, &tmp1, invJac, 3, tmp3, 1, &tmp2,
461  &diffNinvJac[ii * 3], 1);
462  }
463 }
__CLPK_doublereal i
Definition: lapack_wrap.h:47
void cblas_zgemv(const enum CBLAS_ORDER order, const enum CBLAS_TRANSPOSE TransA, const int M, const int N, const void *alpha, const void *A, const int lda, const void *X, const int incX, const void *beta, void *Y, const int incY)
Definition: cblas_zgemv.c:12
__CLPK_doublereal r
Definition: lapack_wrap.h:47

◆ ShapeDiffMBTETQ()

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

Definition at line 885 of file fem_tools.c.

886  {
888  diffN[0] = diffN_MBTETQ0x(x, y, z);
889  diffN[1] = diffN_MBTETQ0y(x, y, z);
890  diffN[2] = diffN_MBTETQ0z(x, y, z);
891  diffN[3] = diffN_MBTETQ1x(x, y, z);
892  diffN[4] = diffN_MBTETQ1y(x, y, z);
893  diffN[5] = diffN_MBTETQ1z(x, y, z);
894  diffN[6] = diffN_MBTETQ2x(x, y, z);
895  diffN[7] = diffN_MBTETQ2y(x, y, z);
896  diffN[8] = diffN_MBTETQ2z(x, y, z);
897  diffN[9] = diffN_MBTETQ3x(x, y, z);
898  diffN[10] = diffN_MBTETQ3y(x, y, z);
899  diffN[11] = diffN_MBTETQ3z(x, y, z);
900  diffN[12] = diffN_MBTETQ4x(x, y, z);
901  diffN[13] = diffN_MBTETQ4y(x, y, z);
902  diffN[14] = diffN_MBTETQ4z(x, y, z);
903  diffN[15] = diffN_MBTETQ5x(x, y, z);
904  diffN[16] = diffN_MBTETQ5y(x, y, z);
905  diffN[17] = diffN_MBTETQ5z(x, y, z);
906  diffN[18] = diffN_MBTETQ6x(x, y, z);
907  diffN[19] = diffN_MBTETQ6y(x, y, z);
908  diffN[20] = diffN_MBTETQ6z(x, y, z);
909  diffN[21] = diffN_MBTETQ7x(x, y, z);
910  diffN[22] = diffN_MBTETQ7y(x, y, z);
911  diffN[23] = diffN_MBTETQ7z(x, y, z);
912  diffN[24] = diffN_MBTETQ8x(x, y, z);
913  diffN[25] = diffN_MBTETQ8y(x, y, z);
914  diffN[26] = diffN_MBTETQ8z(x, y, z);
915  diffN[27] = diffN_MBTETQ9x(x, y, z);
916  diffN[28] = diffN_MBTETQ9y(x, y, z);
917  diffN[29] = diffN_MBTETQ9z(x, y, z);
919 }
#define diffN_MBTETQ4y(x, y, z)
Definition: fem_tools.c:853
#define diffN_MBTETQ2y(x, y, z)
Definition: fem_tools.c:847
#define diffN_MBTETQ6y(x, y, z)
Definition: fem_tools.c:859
#define diffN_MBTETQ2z(x, y, z)
Definition: fem_tools.c:848
#define diffN_MBTETQ6x(x, y, z)
Definition: fem_tools.c:858
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define diffN_MBTETQ7z(x, y, z)
Definition: fem_tools.c:863
#define diffN_MBTETQ3y(x, y, z)
Definition: fem_tools.c:850
#define diffN_MBTETQ8z(x, y, z)
Definition: fem_tools.c:866
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
#define diffN_MBTETQ2x(x, y, z)
Definition: fem_tools.c:846
#define diffN_MBTETQ8x(x, y, z)
Definition: fem_tools.c:864
#define diffN_MBTETQ1x(x, y, z)
Definition: fem_tools.c:843
#define diffN_MBTETQ5y(x, y, z)
Definition: fem_tools.c:856
#define diffN_MBTETQ9x(x, y, z)
Definition: fem_tools.c:867
#define diffN_MBTETQ0x(x, y, z)
Definition: fem_tools.c:840
#define diffN_MBTETQ3z(x, y, z)
Definition: fem_tools.c:851
#define diffN_MBTETQ7y(x, y, z)
Definition: fem_tools.c:862
#define diffN_MBTETQ4z(x, y, z)
Definition: fem_tools.c:854
#define diffN_MBTETQ7x(x, y, z)
Definition: fem_tools.c:861
#define diffN_MBTETQ0y(x, y, z)
Definition: fem_tools.c:841
#define diffN_MBTETQ8y(x, y, z)
Definition: fem_tools.c:865
#define diffN_MBTETQ5z(x, y, z)
Definition: fem_tools.c:857
#define diffN_MBTETQ6z(x, y, z)
Definition: fem_tools.c:860
#define diffN_MBTETQ5x(x, y, z)
Definition: fem_tools.c:855
#define diffN_MBTETQ0z(x, y, z)
Definition: fem_tools.c:842
#define diffN_MBTETQ1z(x, y, z)
Definition: fem_tools.c:845
#define diffN_MBTETQ9y(x, y, z)
Definition: fem_tools.c:868
#define diffN_MBTETQ3x(x, y, z)
Definition: fem_tools.c:849
#define diffN_MBTETQ1y(x, y, z)
Definition: fem_tools.c:844
#define diffN_MBTETQ9z(x, y, z)
Definition: fem_tools.c:869
#define diffN_MBTETQ4x(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 939 of file fem_tools.c.

941  {
943  int ii = 0;
944  for (; ii < G_DIM; ii++) {
945  double x = X[ii], y = Y[ii], z = Z[ii];
946  diffN[30 * ii + 0] = diffN_MBTETQ0x(x, y, z);
947  diffN[30 * ii + 1] = diffN_MBTETQ0y(x, y, z);
948  diffN[30 * ii + 2] = diffN_MBTETQ0z(x, y, z);
949  diffN[30 * ii + 3] = diffN_MBTETQ1x(x, y, z);
950  diffN[30 * ii + 4] = diffN_MBTETQ1y(x, y, z);
951  diffN[30 * ii + 5] = diffN_MBTETQ1z(x, y, z);
952  diffN[30 * ii + 6] = diffN_MBTETQ2x(x, y, z);
953  diffN[30 * ii + 7] = diffN_MBTETQ2y(x, y, z);
954  diffN[30 * ii + 8] = diffN_MBTETQ2z(x, y, z);
955  diffN[30 * ii + 9] = diffN_MBTETQ3x(x, y, z);
956  diffN[30 * ii + 10] = diffN_MBTETQ3y(x, y, z);
957  diffN[30 * ii + 11] = diffN_MBTETQ3z(x, y, z);
958  diffN[30 * ii + 12] = diffN_MBTETQ4x(x, y, z);
959  diffN[30 * ii + 13] = diffN_MBTETQ4y(x, y, z);
960  diffN[30 * ii + 14] = diffN_MBTETQ4z(x, y, z);
961  diffN[30 * ii + 15] = diffN_MBTETQ5x(x, y, z);
962  diffN[30 * ii + 16] = diffN_MBTETQ5y(x, y, z);
963  diffN[30 * ii + 17] = diffN_MBTETQ5z(x, y, z);
964  diffN[30 * ii + 18] = diffN_MBTETQ6x(x, y, z);
965  diffN[30 * ii + 19] = diffN_MBTETQ6y(x, y, z);
966  diffN[30 * ii + 20] = diffN_MBTETQ6z(x, y, z);
967  diffN[30 * ii + 21] = diffN_MBTETQ7x(x, y, z);
968  diffN[30 * ii + 22] = diffN_MBTETQ7y(x, y, z);
969  diffN[30 * ii + 23] = diffN_MBTETQ7z(x, y, z);
970  diffN[30 * ii + 24] = diffN_MBTETQ8x(x, y, z);
971  diffN[30 * ii + 25] = diffN_MBTETQ8y(x, y, z);
972  diffN[30 * ii + 26] = diffN_MBTETQ8z(x, y, z);
973  diffN[30 * ii + 27] = diffN_MBTETQ9x(x, y, z);
974  diffN[30 * ii + 28] = diffN_MBTETQ9y(x, y, z);
975  diffN[30 * ii + 29] = diffN_MBTETQ9z(x, y, z);
976  }
978 }
#define diffN_MBTETQ4y(x, y, z)
Definition: fem_tools.c:853
#define diffN_MBTETQ2y(x, y, z)
Definition: fem_tools.c:847
#define diffN_MBTETQ6y(x, y, z)
Definition: fem_tools.c:859
#define diffN_MBTETQ2z(x, y, z)
Definition: fem_tools.c:848
#define diffN_MBTETQ6x(x, y, z)
Definition: fem_tools.c:858
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define diffN_MBTETQ7z(x, y, z)
Definition: fem_tools.c:863
#define diffN_MBTETQ3y(x, y, z)
Definition: fem_tools.c:850
#define diffN_MBTETQ8z(x, y, z)
Definition: fem_tools.c:866
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
#define diffN_MBTETQ2x(x, y, z)
Definition: fem_tools.c:846
#define diffN_MBTETQ8x(x, y, z)
Definition: fem_tools.c:864
#define diffN_MBTETQ1x(x, y, z)
Definition: fem_tools.c:843
#define diffN_MBTETQ5y(x, y, z)
Definition: fem_tools.c:856
#define diffN_MBTETQ9x(x, y, z)
Definition: fem_tools.c:867
#define diffN_MBTETQ0x(x, y, z)
Definition: fem_tools.c:840
#define diffN_MBTETQ3z(x, y, z)
Definition: fem_tools.c:851
#define diffN_MBTETQ7y(x, y, z)
Definition: fem_tools.c:862
#define diffN_MBTETQ4z(x, y, z)
Definition: fem_tools.c:854
#define diffN_MBTETQ7x(x, y, z)
Definition: fem_tools.c:861
#define diffN_MBTETQ0y(x, y, z)
Definition: fem_tools.c:841
#define diffN_MBTETQ8y(x, y, z)
Definition: fem_tools.c:865
#define diffN_MBTETQ5z(x, y, z)
Definition: fem_tools.c:857
#define diffN_MBTETQ6z(x, y, z)
Definition: fem_tools.c:860
#define diffN_MBTETQ5x(x, y, z)
Definition: fem_tools.c:855
#define diffN_MBTETQ0z(x, y, z)
Definition: fem_tools.c:842
#define diffN_MBTETQ1z(x, y, z)
Definition: fem_tools.c:845
#define diffN_MBTETQ9y(x, y, z)
Definition: fem_tools.c:868
#define diffN_MBTETQ3x(x, y, z)
Definition: fem_tools.c:849
#define diffN_MBTETQ1y(x, y, z)
Definition: fem_tools.c:844
#define diffN_MBTETQ9z(x, y, z)
Definition: fem_tools.c:869
#define diffN_MBTETQ4x(x, y, z)
Definition: fem_tools.c:852

◆ ShapeDiffMBTRI()

PetscErrorCode ShapeDiffMBTRI ( double diffN)

calculate derivatives of shape functions

Definition at line 206 of file fem_tools.c.

206  {
208  diffN[0] = diffN_MBTRI0x;
209  diffN[1] = diffN_MBTRI0y;
210  diffN[2] = diffN_MBTRI1x;
211  diffN[3] = diffN_MBTRI1y;
212  diffN[4] = diffN_MBTRI2x;
213  diffN[5] = diffN_MBTRI2y;
215 }
#define diffN_MBTRI1y
derivative of triangle shape function
Definition: fem_tools.h:60
#define diffN_MBTRI0x
derivative of triangle shape function
Definition: fem_tools.h:57
#define diffN_MBTRI1x
derivative of triangle shape function
Definition: fem_tools.h:59
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define diffN_MBTRI2x
derivative of triangle shape function
Definition: fem_tools.h:61
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
#define diffN_MBTRI2y
derivative of triangle shape function
Definition: fem_tools.h:62
#define diffN_MBTRI0y
derivative of triangle shape function
Definition: fem_tools.h:58

◆ ShapeDiffMBTRIQ()

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

Definition at line 807 of file fem_tools.c.

808  {
810  int ii = 0;
811  for (; ii < G_DIM; ii++) {
812  double x = X[ii], y = Y[ii];
813  diffN[12 * ii + 0] = diffN_MBTRIQ0x(x, y);
814  diffN[12 * ii + 1] = diffN_MBTRIQ0y(x, y);
815  diffN[12 * ii + 2] = diffN_MBTRIQ1x(x, y);
816  diffN[12 * ii + 3] = diffN_MBTRIQ1y(x, y);
817  diffN[12 * ii + 4] = diffN_MBTRIQ2x(x, y);
818  diffN[12 * ii + 5] = diffN_MBTRIQ2y(x, y);
819  diffN[12 * ii + 6] = diffN_MBTRIQ3x(x, y);
820  diffN[12 * ii + 7] = diffN_MBTRIQ3y(x, y);
821  diffN[12 * ii + 8] = diffN_MBTRIQ4x(x, y);
822  diffN[12 * ii + 9] = diffN_MBTRIQ4y(x, y);
823  diffN[12 * ii + 10] = diffN_MBTRIQ5x(x, y);
824  diffN[12 * ii + 11] = diffN_MBTRIQ5y(x, y);
825  }
827 }
#define diffN_MBTRIQ3y(x, y)
Definition: fem_tools.h:98
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define diffN_MBTRIQ4y(x, y)
Definition: fem_tools.h:100
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
#define diffN_MBTRIQ4x(x, y)
Definition: fem_tools.h:99
#define diffN_MBTRIQ1y(x, y)
Definition: fem_tools.h:94
#define diffN_MBTRIQ5x(x, y)
Definition: fem_tools.h:101
#define diffN_MBTRIQ1x(x, y)
Definition: fem_tools.h:93
#define diffN_MBTRIQ2x(x, y)
Definition: fem_tools.h:95
#define diffN_MBTRIQ5y(x, y)
Definition: fem_tools.h:102
#define diffN_MBTRIQ0x(x, y)
Definition: fem_tools.h:91
#define diffN_MBTRIQ3x(x, y)
Definition: fem_tools.h:97
#define diffN_MBTRIQ2y(x, y)
Definition: fem_tools.h:96
#define diffN_MBTRIQ0y(x, y)
Definition: fem_tools.h:92

◆ ShapeFaceBaseMBTRI()

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

Definition at line 216 of file fem_tools.c.

217  {
219 
220  double diffX_ksi[3];
221  double diffX_eta[3];
222  int ii = 0;
223  for (; ii < 3; ii++) {
224  diffX_ksi[ii] = cblas_ddot(3, &coords[ii], 3, &diffN[0], 2);
225  diffX_eta[ii] = cblas_ddot(3, &coords[ii], 3, &diffN[1], 2);
226  }
227  if (s1 != NULL) {
228  cblas_dcopy(3, diffX_ksi, 1, s1, 1);
229  }
230  if (s2 != NULL) {
231  cblas_dcopy(3, diffX_eta, 1, s2, 1);
232  }
233  double Spin_diffX_ksi[9];
234  ierr = Spin(Spin_diffX_ksi, diffX_ksi);
235  CHKERRQ(ierr);
236  cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1., Spin_diffX_ksi, 3,
237  diffX_eta, 1, 0., normal, 1);
239 }
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
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
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: cblas_ddot.c:12
PetscErrorCode Spin(double *spinOmega, double *vecOmega)
calculate spin matrix from vector
Definition: fem_tools.c:558
CHKERRQ(ierr)
static PetscErrorCode ierr
Definition: fem_tools.c:32

◆ ShapeFaceDiffNormalMBTRI()

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

calculate derivative of normal in respect to nodal positions

Definition at line 249 of file fem_tools.c.

250  {
252  // N = Spin(dX/dksi)*dX/deta = -Spin(dX/deta)*dX/dksi
253 
254  double diffX_ksi[3];
255  double diffX_eta[3];
256  int ii = 0;
257  for (; ii < 3; ii++) {
258  diffX_ksi[ii] = cblas_ddot(3, &coords[ii], 3, &diffN[0], 2);
259  diffX_eta[ii] = cblas_ddot(3, &coords[ii], 3, &diffN[1], 2);
260  }
261  double Spin_diffX_ksi[9];
262  ierr = Spin(Spin_diffX_ksi, diffX_ksi);
263  CHKERRQ(ierr);
264  double Spin_diffX_eta[9];
265  ierr = Spin(Spin_diffX_eta, diffX_eta);
266  CHKERRQ(ierr);
267  double B_ksi[3 * 9];
268  bzero(B_ksi, 3 * 9 * sizeof(double));
269  double B_eta[3 * 9];
270  bzero(B_eta, 3 * 9 * sizeof(double));
271  // B_ksi[] = [
272  // diffN[2*0+0], 0, 0, diffN[2*1+0], 0, 0, diffN[2*2+0], 0,
273  // 0
274  // 0, diffN[2*0+0], 0, 0, diffN[2*1+0], 0, 0, diffN[2*2+0],
275  // 0
276  // 0, 0, diffM[2*0+0], 0, 0, diffN[2*1+0], 0, 0,
277  // diffN[2*2+0]
278  //]
279  // B_eta[] = [
280  // diffN[2*0+1], 0, 0, diffN[2*1+1], 0, 0, diffN[2*2+1], 0,
281  // 0
282  // 0, diffN[2*0+1], 0, 0, diffN[2*1+1], 0, 0, diffN[2*2+1],
283  // 0
284  // 0, 0, diffM[2*0+1], 0, 0, diffN[2*1+1], 0, 0,
285  // diffN[2*2+1]
286  //]
287  ii = 0;
288  for (; ii < 3; ii++) {
289  cblas_dcopy(3, &diffN[0], 2, &B_ksi[ii * 9 + ii], 3);
290  cblas_dcopy(3, &diffN[1], 2, &B_eta[ii * 9 + ii], 3);
291  }
293  Spin_diffX_ksi, 3, B_eta, 9, 0., diff_normal, 9);
295  Spin_diffX_eta, 3, B_ksi, 9, 1., diff_normal, 9);
297 }
#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
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: cblas_ddot.c:12
PetscErrorCode Spin(double *spinOmega, double *vecOmega)
calculate spin matrix from vector
Definition: fem_tools.c:558
CHKERRQ(ierr)
void cblas_dgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_TRANSPOSE TransB, const int M, const int N, const int K, const double alpha, const double *A, const int lda, const double *B, const int ldb, const double beta, double *C, const int ldc)
Definition: cblas_dgemm.c:12
static PetscErrorCode ierr
Definition: fem_tools.c:32

◆ 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 241 of file fem_tools.c.

242  {
244  ierr = ShapeFaceBaseMBTRI(diffN, coords, normal, NULL, NULL);
245  CHKERRQ(ierr);
247 }
#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
CHKERRQ(ierr)
PetscErrorCode ShapeFaceBaseMBTRI(double *diffN, const double *coords, double *normal, double *s1, double *s2)
Definition: fem_tools.c:216
static PetscErrorCode ierr
Definition: fem_tools.c:32

◆ ShapeFaceNormalMBTRI_complex()

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

Definition at line 464 of file fem_tools.c.

466  {
468  double complex diffX_x, diffX_y, diffX_z;
469  double complex diffY_x, diffY_y, diffY_z;
470  diffX_x = diffX_y = diffX_z = 0.;
471  diffY_x = diffY_y = diffY_z = 0.;
472  int ii;
473  for (ii = 0; ii < 3; ii++) {
474  diffX_x +=
475  (xcoords[3 * ii + 0].r + I * xcoords[3 * ii + 0].i) * diffN[2 * ii + 0];
476  diffX_y +=
477  (xcoords[3 * ii + 1].r + I * xcoords[3 * ii + 1].i) * diffN[2 * ii + 0];
478  diffX_z +=
479  (xcoords[3 * ii + 2].r + I * xcoords[3 * ii + 2].i) * diffN[2 * ii + 0];
480  diffY_x +=
481  (xcoords[3 * ii + 0].r + I * xcoords[3 * ii + 0].i) * diffN[2 * ii + 1];
482  diffY_y +=
483  (xcoords[3 * ii + 1].r + I * xcoords[3 * ii + 1].i) * diffN[2 * ii + 1];
484  diffY_z +=
485  (xcoords[3 * ii + 2].r + I * xcoords[3 * ii + 2].i) * diffN[2 * ii + 1];
486  }
487  double complex tmp;
488  tmp = diffX_y * diffY_z - diffX_z * diffY_y;
489  xnormal[0].r = creal(tmp);
490  xnormal[0].i = cimag(tmp);
491  tmp = diffX_z * diffY_x - diffX_x * diffY_z;
492  xnormal[1].r = creal(tmp);
493  xnormal[1].i = cimag(tmp);
494  tmp = diffX_x * diffY_y - diffX_y * diffY_x;
495  xnormal[2].r = creal(tmp);
496  xnormal[2].i = cimag(tmp);
498 }
#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
__CLPK_doublereal i
Definition: lapack_wrap.h:47
__CLPK_doublereal r
Definition: lapack_wrap.h:47

◆ ShapeInvJacVolume()

PetscErrorCode ShapeInvJacVolume ( double jac)

Definition at line 51 of file fem_tools.c.

51  {
53  __CLPK_integer ipiv[4];
54  __CLPK_doublereal work[3];
55  __CLPK_integer lwork = 3;
56  __CLPK_integer info;
57  info = lapack_dgetrf(3, 3, jac, 3, ipiv);
58  if (info != 0)
59  SETERRQ1(PETSC_COMM_SELF, 1, "info = %d", info);
60  info = lapack_dgetri(3, jac, 3, ipiv, work, lwork);
61  if (info != 0)
62  SETERRQ1(PETSC_COMM_SELF, 1, "info = %d", info);
64 }
#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
static __CLPK_integer lapack_dgetri(__CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *work, __CLPK_integer lwork)
Definition: lapack_wrap.h:197
static __CLPK_integer lapack_dgetrf(__CLPK_integer m, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv)
Definition: lapack_wrap.h:169
long int __CLPK_integer
Definition: lapack_wrap.h:35

◆ ShapeJacMBTET()

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

calculate jacobian

Definition at line 300 of file fem_tools.c.

300  {
302  int ii, jj, kk;
303  bzero(jac, sizeof(double) * 9);
304  for (ii = 0; ii < 4; ii++) // shape func.
305  for (jj = 0; jj < 3; jj++) // space
306  for (kk = 0; kk < 3; kk++) // derivative of shape func.
307  jac[jj * 3 + kk] += diffN[ii * 3 + kk] * coords[ii * 3 + jj];
309 }
#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

◆ ShapeJacMBTETQ()

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

Definition at line 979 of file fem_tools.c.

980  {
982  int ii, jj, kk;
983  bzero(Jac, sizeof(double) * 9);
984  for (ii = 0; ii < 10; ii++) // shape func.
985  for (jj = 0; jj < 3; jj++) // space
986  for (kk = 0; kk < 3; kk++) // direvative of shape func.
987  Jac[jj * 3 + kk] += diffN[ii * 3 + kk] * coords[ii * 3 + jj];
989 }
#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

◆ ShapeMBEDGE()

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

Definition at line 773 of file fem_tools.c.

773  {
775  int ii = 0;
776  for (; ii < DIM; ii++) {
777  double x = G_X[ii];
778  N[2 * ii + 0] = N_MBEDGE0(x);
779  N[2 * ii + 1] = N_MBEDGE1(x);
780  }
782 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define N_MBEDGE1(x)
edge shape function
Definition: fem_tools.h:80
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
#define N_MBEDGE0(x)
edge shape function
Definition: fem_tools.h:79
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 318 of file fem_tools.c.

319  {
321  int ii = 0;
322  for (; ii < DIM; ii++) {
323  double x = G_X[ii], y = G_Y[ii], z = G_Z[ii];
324  N[4 * ii + 0] = N_MBTET0(x, y, z);
325  N[4 * ii + 1] = N_MBTET1(x, y, z);
326  N[4 * ii + 2] = N_MBTET2(x, y, z);
327  N[4 * ii + 3] = N_MBTET3(x, y, z);
328  }
330 }
#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 N_MBTET3(x, y, z)
tetrahedral shape function
Definition: fem_tools.h:39
#define N_MBTET0(x, y, z)
tetrahedral shape function
Definition: fem_tools.h:36
#define N_MBTET2(x, y, z)
tetrahedral shape function
Definition: fem_tools.h:38
#define N_MBTET1(x, y, z)
tetrahedral shape function
Definition: fem_tools.h:37
const int N
Definition: speed_test.cpp:3

◆ 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 347 of file fem_tools.c.

350  {
352  double A[3 * 3];
353  int IPIV[3];
354  // COL MAJOR
355  // X
356  A[0 + 3 * 0] =
357  cblas_ddot(4, &diffN[0 * 3 + 0], 3, &elem_coords[0 * 3 + 0], 3);
358  A[0 + 3 * 1] =
359  cblas_ddot(4, &diffN[0 * 3 + 1], 3, &elem_coords[0 * 3 + 0], 3);
360  A[0 + 3 * 2] =
361  cblas_ddot(4, &diffN[0 * 3 + 2], 3, &elem_coords[0 * 3 + 0], 3);
362  loc_coords[0] =
363  glob_coords[0] - cblas_ddot(4, &N[0], 1, &elem_coords[0 * 3 + 0], 3);
364  // printf("A\n[ %3.2f %3.2f %3.2f ] %3.2f \n",A[0*3],A[1*3],A[2*3],R[0]);
365  // Y
366  A[1 + 3 * 0] =
367  cblas_ddot(4, &diffN[0 * 3 + 0], 3, &elem_coords[0 * 3 + 1], 3);
368  A[1 + 3 * 1] =
369  cblas_ddot(4, &diffN[0 * 3 + 1], 3, &elem_coords[0 * 3 + 1], 3);
370  A[1 + 3 * 2] =
371  cblas_ddot(4, &diffN[0 * 3 + 2], 3, &elem_coords[0 * 3 + 1], 3);
372  loc_coords[1] =
373  glob_coords[1] - cblas_ddot(4, &N[0], 1, &elem_coords[0 * 3 + 1], 3);
374  // printf("[ %3.2f %3.2f %3.2f ] %3.2f \n",A[1+3*0],A[1+3*1],A[1+3*2],R[1]);
375  // Z
376  A[2 + 3 * 0] =
377  cblas_ddot(4, &diffN[0 * 3 + 0], 3, &elem_coords[0 * 3 + 2], 3);
378  A[2 + 3 * 1] =
379  cblas_ddot(4, &diffN[0 * 3 + 1], 3, &elem_coords[0 * 3 + 2], 3);
380  A[2 + 3 * 2] =
381  cblas_ddot(4, &diffN[0 * 3 + 2], 3, &elem_coords[0 * 3 + 2], 3);
382  loc_coords[2] =
383  glob_coords[2] - cblas_ddot(4, &N[0], 1, &elem_coords[0 * 3 + 2], 3);
384  // printf("[ %3.2f %3.2f %3.2f ] %3.2f \n",A[2+3*0],A[2+3*1],A[2+3*2],R[1]);
385  int info =
386  lapack_dgesv(3, 1, &A[0], 3, (__CLPK_integer *)IPIV, loc_coords, 3);
387  if (info != 0)
388  SETERRQ1(PETSC_COMM_SELF, 1, "info == %d", info);
390 }
#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
const int N
Definition: speed_test.cpp:3
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition: lapack_wrap.h:188
long int __CLPK_integer
Definition: lapack_wrap.h:35

◆ ShapeMBTETQ()

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

Definition at line 870 of file fem_tools.c.

871  {
873  N[0] = N_MBTETQ0(x, y, z);
874  N[1] = N_MBTETQ1(x, y, z);
875  N[2] = N_MBTETQ2(x, y, z);
876  N[3] = N_MBTETQ3(x, y, z);
877  N[4] = N_MBTETQ4(x, y, z);
878  N[5] = N_MBTETQ5(x, y, z);
879  N[6] = N_MBTETQ6(x, y, z);
880  N[7] = N_MBTETQ7(x, y, z);
881  N[8] = N_MBTETQ8(x, y, z);
882  N[9] = N_MBTETQ9(x, y, z);
884 }
#define N_MBTETQ5(x, y, z)
Definition: fem_tools.c:835
#define N_MBTETQ8(x, y, z)
Definition: fem_tools.c:838
#define N_MBTETQ0(x, y, z)
Definition: fem_tools.c:830
#define N_MBTETQ2(x, y, z)
Definition: fem_tools.c:832
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define N_MBTETQ1(x, y, z)
Definition: fem_tools.c:831
#define N_MBTETQ7(x, y, z)
Definition: fem_tools.c:837
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
#define N_MBTETQ6(x, y, z)
Definition: fem_tools.c:836
#define N_MBTETQ3(x, y, z)
Definition: fem_tools.c:833
#define N_MBTETQ9(x, y, z)
Definition: fem_tools.c:839
const int N
Definition: speed_test.cpp:3
#define N_MBTETQ4(x, y, z)
Definition: fem_tools.c:834

◆ 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 991 of file fem_tools.c.

993  {
995 
996  double Jac[9];
997  int ii = 0;
998  for (; ii < G_DIM; ii++) {
999  ierr = ShapeJacMBTETQ(&diffN[30 * ii], coords, Jac);
1000  CHKERRQ(ierr);
1001  detJac_at_Gauss_Points[ii] = ShapeDetJacVolume(Jac);
1002  }
1004 }
#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
CHKERRQ(ierr)
PetscErrorCode ShapeJacMBTETQ(const double *diffN, const double *coords, double *Jac)
Definition: fem_tools.c:979
static PetscErrorCode ierr
Definition: fem_tools.c:32
double ShapeDetJacVolume(double *jac)
determined of jacobian
Definition: fem_tools.c:34

◆ ShapeMBTETQ_GAUSS()

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

Definition at line 920 of file fem_tools.c.

921  {
923  int ii = 0;
924  for (; ii < G_DIM; ii++) {
925  double x = X[ii], y = Y[ii], z = Z[ii];
926  N[10 * ii + 0] = N_MBTETQ0(x, y, z);
927  N[10 * ii + 1] = N_MBTETQ1(x, y, z);
928  N[10 * ii + 2] = N_MBTETQ2(x, y, z);
929  N[10 * ii + 3] = N_MBTETQ3(x, y, z);
930  N[10 * ii + 4] = N_MBTETQ4(x, y, z);
931  N[10 * ii + 5] = N_MBTETQ5(x, y, z);
932  N[10 * ii + 6] = N_MBTETQ6(x, y, z);
933  N[10 * ii + 7] = N_MBTETQ7(x, y, z);
934  N[10 * ii + 8] = N_MBTETQ8(x, y, z);
935  N[10 * ii + 9] = N_MBTETQ9(x, y, z);
936  }
938 }
#define N_MBTETQ5(x, y, z)
Definition: fem_tools.c:835
#define N_MBTETQ8(x, y, z)
Definition: fem_tools.c:838
#define N_MBTETQ0(x, y, z)
Definition: fem_tools.c:830
#define N_MBTETQ2(x, y, z)
Definition: fem_tools.c:832
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define N_MBTETQ1(x, y, z)
Definition: fem_tools.c:831
#define N_MBTETQ7(x, y, z)
Definition: fem_tools.c:837
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
#define N_MBTETQ6(x, y, z)
Definition: fem_tools.c:836
#define N_MBTETQ3(x, y, z)
Definition: fem_tools.c:833
#define N_MBTETQ9(x, y, z)
Definition: fem_tools.c:839
const int N
Definition: speed_test.cpp:3
#define N_MBTETQ4(x, y, z)
Definition: fem_tools.c:834

◆ 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 1019 of file fem_tools.c.

1022  {
1024  double A[3 * 3];
1025  double R[3];
1026  int IPIV[3];
1027  float NORM_dR = 1000.;
1028  float NORM_R0;
1029  ShapeMBTETQ(N, 0.1, 0.1, 0.1);
1030  ShapeDiffMBTETQ(diffN, 0.1, 0.1, 0.1);
1031  R[0] = glob_coords[0] - cblas_ddot(10, &N[0], 1, &elem_coords[0], 3);
1032  R[1] = glob_coords[1] - cblas_ddot(10, &N[0], 1, &elem_coords[1], 3);
1033  R[2] = glob_coords[2] - cblas_ddot(10, &N[0], 1, &elem_coords[2], 3);
1034  NORM_R0 = cblas_dnrm2(3, &R[0], 1);
1035  while ((NORM_dR / NORM_R0) > eps) {
1036  // COL MAJOR
1037  // X
1038  A[0 + 3 * 0] = cblas_ddot(10, &diffN[0 * 3 + 0], 3, &elem_coords[0], 3);
1039  A[0 + 3 * 1] = cblas_ddot(10, &diffN[0 * 3 + 1], 3, &elem_coords[0], 3);
1040  A[0 + 3 * 2] = cblas_ddot(10, &diffN[0 * 3 + 2], 3, &elem_coords[0], 3);
1041  R[0] = glob_coords[0] - cblas_ddot(10, &N[0], 1, &elem_coords[0], 3);
1042  // Y
1043  A[1 + 3 * 0] = cblas_ddot(10, &diffN[0 * 3 + 0], 3, &elem_coords[1], 3);
1044  A[1 + 3 * 1] = cblas_ddot(10, &diffN[0 * 3 + 1], 3, &elem_coords[1], 3);
1045  A[1 + 3 * 2] = cblas_ddot(10, &diffN[0 * 3 + 2], 3, &elem_coords[1], 3);
1046  R[1] = glob_coords[1] - cblas_ddot(10, &N[0], 1, &elem_coords[1], 3);
1047  // Z
1048  A[2 + 3 * 0] =
1049  cblas_ddot(10, &diffN[0 * 3 + 0], 3, &elem_coords[0 * 3 + 2], 3);
1050  A[2 + 3 * 1] =
1051  cblas_ddot(10, &diffN[0 * 3 + 1], 3, &elem_coords[0 * 3 + 2], 3);
1052  A[2 + 3 * 2] =
1053  cblas_ddot(10, &diffN[0 * 3 + 2], 3, &elem_coords[0 * 3 + 2], 3);
1054  R[2] = glob_coords[2] - cblas_ddot(10, &N[0], 1, &elem_coords[2], 3);
1055  int info = lapack_dgesv(3, 1, &A[0], 3, (__CLPK_integer *)IPIV, R, 3);
1056  assert(info == 0);
1057  NOT_USED(info);
1058  cblas_daxpy(3, 1., R, 1, loc_coords, 1);
1059  NORM_dR = cblas_dnrm2(3, &R[0], 1);
1060  ShapeMBTETQ(N, loc_coords[0], loc_coords[1], loc_coords[2]);
1061  ShapeDiffMBTETQ(diffN, loc_coords[0], loc_coords[1], loc_coords[2]);
1062  }
1064 }
PetscErrorCode ShapeMBTETQ(double *N, const double x, const double y, const double z)
Definition: fem_tools.c:870
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
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
PetscErrorCode ShapeDiffMBTETQ(double *diffN, const double x, const double y, const double z)
Definition: fem_tools.c:885
double cblas_dnrm2(const int N, const double *X, const int incX)
Definition: cblas_dnrm2.c:12
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: cblas_ddot.c:12
static const double eps
const int N
Definition: speed_test.cpp:3
#define NOT_USED(x)
Definition: definitions.h:303
static __CLPK_integer lapack_dgesv(__CLPK_integer n, __CLPK_integer nrhs, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_integer *ipiv, __CLPK_doublereal *b, __CLPK_integer ldb)
Definition: lapack_wrap.h:188
long int __CLPK_integer
Definition: lapack_wrap.h:35

◆ 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 194 of file fem_tools.c.

195  {
197  int ii = 0;
198  for (; ii < G_DIM; ii++) {
199  double x = X[ii], y = Y[ii];
200  N[3 * ii + 0] = N_MBTRI0(x, y);
201  N[3 * ii + 1] = N_MBTRI1(x, y);
202  N[3 * ii + 2] = N_MBTRI2(x, y);
203  }
205 }
#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 N_MBTRI2(x, y)
triangle shape function
Definition: fem_tools.h:56
const int N
Definition: speed_test.cpp:3
#define N_MBTRI1(x, y)
triangle shape function
Definition: fem_tools.h:55
#define N_MBTRI0(x, y)
triangle shape function
Definition: fem_tools.h:54

◆ 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 392 of file fem_tools.c.

395  {
397  double A[2 * 2];
398 
399  // 1st and 2nd element of matrix A
400  A[0] = cblas_ddot(3, &diffN[0], 2, &elem_coords[0], 2); // dot product
401  A[1] = cblas_ddot(3, &diffN[1], 2, &elem_coords[0], 2);
402  loc_coords[0] = glob_coords[0] - cblas_ddot(3, &N[0], 1, &elem_coords[0], 2);
403 
404  // 3rd and 4th element of matrix A
405  A[2] = cblas_ddot(3, &diffN[0], 2, &elem_coords[1], 2);
406  A[3] = cblas_ddot(3, &diffN[1], 2, &elem_coords[1], 2);
407  loc_coords[1] = glob_coords[1] - cblas_ddot(3, &N[0], 1, &elem_coords[1], 2);
408 
409  // calculate directly the solution (as the size of matrix is only 2x2)
410  double invA[2 * 2], detA;
411  detA = A[0] * A[3] - A[1] * A[2];
412  detA = 1.0 / detA;
413  invA[0] = A[3] * detA;
414  invA[1] = -1.0 * A[1] * detA;
415  invA[2] = -1.0 * A[2] * detA;
416  invA[3] = A[0] * detA;
417 
418  double loc_coords_new[2];
419  loc_coords_new[0] = invA[0] * loc_coords[0] + invA[1] * loc_coords[1];
420  loc_coords_new[1] = invA[2] * loc_coords[0] + invA[3] * loc_coords[1];
421 
422  loc_coords[0] = loc_coords_new[0];
423  loc_coords[1] = loc_coords_new[1];
425 }
#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
const int N
Definition: speed_test.cpp:3

◆ ShapeMBTRIQ()

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

Definition at line 792 of file fem_tools.c.

793  {
795  int ii = 0;
796  for (; ii < G_DIM; ii++) {
797  double x = X[ii], y = Y[ii];
798  N[6 * ii + 0] = N_MBTRIQ0(x, y);
799  N[6 * ii + 1] = N_MBTRIQ1(x, y);
800  N[6 * ii + 2] = N_MBTRIQ2(x, y);
801  N[6 * ii + 3] = N_MBTRIQ3(x, y);
802  N[6 * ii + 4] = N_MBTRIQ4(x, y);
803  N[6 * ii + 5] = N_MBTRIQ5(x, y);
804  }
806 }
#define N_MBTRIQ0(x, y)
Definition: fem_tools.h:85
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
#define N_MBTRIQ2(x, y)
Definition: fem_tools.h:87
#define N_MBTRIQ5(x, y)
Definition: fem_tools.h:90
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
#define N_MBTRIQ4(x, y)
Definition: fem_tools.h:89
#define N_MBTRIQ1(x, y)
Definition: fem_tools.h:86
#define N_MBTRIQ3(x, y)
Definition: fem_tools.h:88
const int N
Definition: speed_test.cpp:3

◆ ShapeVolumeMBTET()

double ShapeVolumeMBTET ( double diffN,
const double coords 
)

calculate TET volume

Definition at line 310 of file fem_tools.c.

310  {
311  double Jac[9];
312  ShapeJacMBTET(diffN, coords, Jac);
313  double detJac = ShapeDetJacVolume(Jac);
314  // printf("detJac = +%6.4e\n",detJac);
315  // print_mat(Jac,3,3);
316  return detJac * G_TET_W1[0] / 6.;
317 }
PetscErrorCode ShapeJacMBTET(double *diffN, const double *coords, double *jac)
calculate jacobian
Definition: fem_tools.c:300
static const double G_TET_W1[]
Definition: fem_tools.h:495
double ShapeDetJacVolume(double *jac)
determined of jacobian
Definition: fem_tools.c:34

◆ ShapeVolumeMBTETQ()

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

Definition at line 1005 of file fem_tools.c.

1006  {
1007 
1008  int ii = 0;
1009  double vol = 0;
1010  double detJac_at_Gauss_Points[G_DIM];
1011  ierr = ShapeMBTETQ_detJac_at_Gauss_Points(detJac_at_Gauss_Points, diffN,
1012  coords, G_DIM);
1013  CHKERRQ(ierr);
1014  for (; ii < G_DIM; ii++) {
1015  vol += G_TET_W[ii] * (detJac_at_Gauss_Points[ii]) / 6;
1016  }
1017  return vol;
1018 }
CHKERRQ(ierr)
PetscErrorCode ShapeMBTETQ_detJac_at_Gauss_Points(double *detJac_at_Gauss_Points, const double *diffN, const double *coords, int G_DIM)
Definition: fem_tools.c:991
static PetscErrorCode ierr
Definition: fem_tools.c:32

◆ Spin()

PetscErrorCode Spin ( double spinOmega,
double vecOmega 
)

calculate spin matrix from vector

Definition at line 558 of file fem_tools.c.

558  {
560  bzero(spinOmega, 9 * sizeof(double));
561  spinOmega[0 * 3 + 1] = -vecOmega[2];
562  spinOmega[0 * 3 + 2] = +vecOmega[1];
563  spinOmega[1 * 3 + 0] = +vecOmega[2];
564  spinOmega[1 * 3 + 2] = -vecOmega[0];
565  spinOmega[2 * 3 + 0] = -vecOmega[1];
566  spinOmega[2 * 3 + 1] = +vecOmega[0];
568 }
#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

Variable Documentation

◆ ierr

PetscErrorCode ierr
static

Definition at line 32 of file fem_tools.c.