v0.14.0
Functions | Variables
edge_and_bubble_shape_functions_on_quad.cpp File Reference
#include <MoFEM.hpp>

Go to the source code of this file.

Functions

static double sum_matrix (MatrixDouble &m)
 
int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)
Examples
edge_and_bubble_shape_functions_on_quad.cpp.

Definition at line 22 of file edge_and_bubble_shape_functions_on_quad.cpp.

22  {
23 
24  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
25  try {
26  moab::Core mb_instance;
27  moab::Interface &moab = mb_instance;
28 
29  MoFEM::Core core(moab);
30 
31  PetscBool flg = PETSC_TRUE;
32  char mesh_file_name[255];
33 #if PETSC_VERSION_GE(3, 6, 4)
34  CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
35  255, &flg);
36 #else
37  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
38  mesh_file_name, 255, &flg);
39 #endif
40  if (flg != PETSC_TRUE) {
41  SETERRQ(PETSC_COMM_SELF, 1, "Error -my_file (mesh file needed)");
42  }
43 
44  const char *option;
45  option = "";
46  CHKERR moab.load_file(mesh_file_name, 0, option);
47  Range verts;
48  CHKERR moab.get_entities_by_type(0, MBVERTEX, verts, true);
49 
50  MatrixDouble coords(verts.size(), 3);
51  MatrixDouble N(verts.size(), 4);
52  MatrixDouble diffN(verts.size(), 8);
53  CHKERR moab.get_coords(verts, &coords(0, 0));
54  for (int i = 0; i < verts.size(); ++i) {
55  N(i, 0) = N_MBQUAD0(coords(i, 0), coords(i, 1));
56  N(i, 1) = N_MBQUAD1(coords(i, 0), coords(i, 1));
57  N(i, 2) = N_MBQUAD2(coords(i, 0), coords(i, 1));
58  N(i, 3) = N_MBQUAD3(coords(i, 0), coords(i, 1));
59 
60  diffN(i, 0) = diffN_MBQUAD0x(coords(i, 1));
61  diffN(i, 1) = diffN_MBQUAD0y(coords(i, 0));
62  diffN(i, 2) = diffN_MBQUAD1x(coords(i, 1));
63  diffN(i, 3) = diffN_MBQUAD1y(coords(i, 0));
64  diffN(i, 4) = diffN_MBQUAD2x(coords(i, 1));
65  diffN(i, 5) = diffN_MBQUAD2y(coords(i, 0));
66  diffN(i, 6) = diffN_MBQUAD3x(coords(i, 1));
67  diffN(i, 7) = diffN_MBQUAD3y(coords(i, 0));
68  }
69 
70  /* BUBBLES */
71  {
72  std::array<int, 4> faces_nodes = {0, 1, 2, 3};
73  int p = 7;
74  int P = NBFACEQUAD_H1(p);
75  MatrixDouble quad_bubbles(verts.size(), P);
76  MatrixDouble quad_diff_bubbles(verts.size(), P * 2);
77  constexpr double eps = 1e-4;
78  double quad_bubbles_sum = 3.275491e+01;
79  double quad_diff_bubbles_sum = 5.838131e+02;
80 
81  CHKERR H1_QuadShapeFunctions_MBQUAD(
82  faces_nodes.data(), p, &*N.data().begin(), &*diffN.data().begin(),
83  &*quad_bubbles.data().begin(), &*quad_diff_bubbles.data().begin(),
84  verts.size(), Legendre_polynomials);
85 
86  if (std::abs(quad_bubbles_sum - sum_matrix(quad_bubbles)) > eps) {
87  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
88  "wrong result = %8.6e", sum_matrix(quad_bubbles));
89  }
90  if (std::abs(quad_diff_bubbles_sum - sum_matrix(quad_diff_bubbles)) >
91  eps) {
92  SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
93  "wrong result = %8.6e", sum_matrix(quad_diff_bubbles));
94  }
95  }
96 
97  /* EDGES */
98  {
99  int sense[] = {1, -1, 1, -1};
100  int p[] = {3, 4, 5, 6};
101  int P[4];
102  int ee = 0;
103  for (; ee < 4; ee++) {
104  P[ee] = NBEDGE_H1(p[ee]);
105  }
106 
107  double *quad_edges_ptr[4];
108  std::array<MatrixDouble, 4> quad_edges;
109  for (auto ee : {0, 1, 2, 3}) {
110  quad_edges[ee] = MatrixDouble(verts.size(), P[ee]);
111  quad_edges_ptr[ee] = &quad_edges[ee](0, 0);
112  }
113 
114  double *quad_diff_edges_ptr[4];
115  std::array<MatrixDouble, 4> quad_diff_edges;
116  for (auto ee : {0, 1, 2, 3}) {
117  quad_diff_edges[ee] = MatrixDouble(verts.size(), P[ee] * 2);
118  quad_diff_edges_ptr[ee] = &quad_diff_edges[ee](0, 0);
119  }
120 
121  double eps = 1e-4;
122  double quad_edges_sum[] = {1.237500e+01, 1.535050e+01, 1.781890e+01,
123  2.015678e+01};
124  double quad_diff_edges_sum[] = {8.470000e+01, 1.192510e+02, 1.580128e+02,
125  1.976378e+02};
126 
127  CHKERR H1_EdgeShapeFunctions_MBQUAD(
128  sense, p, &*N.data().begin(), &*diffN.data().begin(), quad_edges_ptr,
129  quad_diff_edges_ptr, verts.size(), Legendre_polynomials);
130 
131  for (int ee = 0; ee < 4; ++ee) {
132  if (std::abs(quad_edges_sum[ee] - sum_matrix(quad_edges[ee])) > eps)
133  SETERRQ2(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
134  "edge %d wrong result = %8.6e", ee,
135  sum_matrix(quad_edges[ee]));
136 
137  if (std::abs(quad_diff_edges_sum[ee] -
138  sum_matrix(quad_diff_edges[ee])) > eps)
139  SETERRQ2(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
140  "edge %d wrong result = %8.6e", ee,
141  sum_matrix(quad_diff_edges[ee]));
142  }
143  }
144  }
145  CATCH_ERRORS;
146 
147  CHKERR MoFEM::Core::Finalize();
148 }

◆ sum_matrix()

static double sum_matrix ( MatrixDouble &  m)
static
Examples
edge_and_bubble_shape_functions_on_quad.cpp, and prism_elements_from_surface.cpp.

Definition at line 12 of file edge_and_bubble_shape_functions_on_quad.cpp.

12  {
13  double s = 0;
14  for (unsigned int ii = 0; ii < m.size1(); ii++) {
15  for (unsigned int jj = 0; jj < m.size2(); jj++) {
16  s += std::abs(m(ii, jj));
17  }
18  }
19  return s;
20 }

Variable Documentation

◆ help

char help[] = "...\n\n"
static
Examples
edge_and_bubble_shape_functions_on_quad.cpp.

Definition at line 10 of file edge_and_bubble_shape_functions_on_quad.cpp.

MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
NBEDGE_H1
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
Definition: h1_hdiv_hcurl_l2.h:55
diffN_MBQUAD1x
#define diffN_MBQUAD1x(y)
Definition: fem_tools.h:63
H1_EdgeShapeFunctions_MBQUAD
PetscErrorCode H1_EdgeShapeFunctions_MBQUAD(int *sense, int *p, double *N, double *diffN, double *edgeN[4], double *diff_edgeN[4], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: h1.c:1091
diffN_MBQUAD1y
#define diffN_MBQUAD1y(x)
Definition: fem_tools.h:64
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
H1_QuadShapeFunctions_MBQUAD
PetscErrorCode H1_QuadShapeFunctions_MBQUAD(int *faces_nodes, int p, double *N, double *diffN, double *faceN, double *diff_faceN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: h1.c:959
diffN_MBQUAD2x
#define diffN_MBQUAD2x(y)
Definition: fem_tools.h:65
help
static char help[]
Definition: edge_and_bubble_shape_functions_on_quad.cpp:10
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
N_MBQUAD2
#define N_MBQUAD2(x, y)
quad shape function
Definition: fem_tools.h:59
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2002
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
diffN_MBQUAD3x
#define diffN_MBQUAD3x(y)
Definition: fem_tools.h:67
EshelbianPlasticity::P
@ P
Definition: EshelbianContact.cpp:193
N_MBQUAD1
#define N_MBQUAD1(x, y)
quad shape function
Definition: fem_tools.h:58
diffN_MBQUAD2y
#define diffN_MBQUAD2y(x)
Definition: fem_tools.h:66
Legendre_polynomials
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Legendre approximation basis.
Definition: base_functions.c:15
NBFACEQUAD_H1
#define NBFACEQUAD_H1(P)
Number of base functions on quad for H1 space.
Definition: h1_hdiv_hcurl_l2.h:65
diffN_MBQUAD0x
#define diffN_MBQUAD0x(y)
Definition: fem_tools.h:61
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
diffN_MBQUAD0y
#define diffN_MBQUAD0y(x)
Definition: fem_tools.h:62
diffN_MBQUAD3y
#define diffN_MBQUAD3y(x)
Definition: fem_tools.h:68
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
N
const int N
Definition: speed_test.cpp:3
Range
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1140
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
N_MBQUAD0
#define N_MBQUAD0(x, y)
quad shape function
Definition: fem_tools.h:57
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
N_MBQUAD3
#define N_MBQUAD3(x, y)
quad shape function
Definition: fem_tools.h:60
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
sum_matrix
static double sum_matrix(MatrixDouble &m)
Definition: edge_and_bubble_shape_functions_on_quad.cpp:12
Generated by Doxygen 1.8.17 and hosted at University of Glasgow