v0.13.1
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 36 of file edge_and_bubble_shape_functions_on_quad.cpp.

36 {
37
38 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
39 try {
40 moab::Core mb_instance;
41 moab::Interface &moab = mb_instance;
42
43 MoFEM::Core core(moab);
44
45 PetscBool flg = PETSC_TRUE;
46 char mesh_file_name[255];
47#if PETSC_VERSION_GE(3, 6, 4)
48 CHKERR PetscOptionsGetString(PETSC_NULL, "", "-my_file", mesh_file_name,
49 255, &flg);
50#else
51 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
52 mesh_file_name, 255, &flg);
53#endif
54 if (flg != PETSC_TRUE) {
55 SETERRQ(PETSC_COMM_SELF, 1, "Error -my_file (mesh file needed)");
56 }
57
58 const char *option;
59 option = "";
60 CHKERR moab.load_file(mesh_file_name, 0, option);
61 Range verts;
62 CHKERR moab.get_entities_by_type(0, MBVERTEX, verts, true);
63
64 MatrixDouble coords(verts.size(), 3);
65 MatrixDouble N(verts.size(), 4);
66 MatrixDouble diffN(verts.size(), 8);
67 CHKERR moab.get_coords(verts, &coords(0, 0));
68 for (int i = 0; i < verts.size(); ++i) {
69 N(i, 0) = N_MBQUAD0(coords(i, 0), coords(i, 1));
70 N(i, 1) = N_MBQUAD1(coords(i, 0), coords(i, 1));
71 N(i, 2) = N_MBQUAD2(coords(i, 0), coords(i, 1));
72 N(i, 3) = N_MBQUAD3(coords(i, 0), coords(i, 1));
73
74 diffN(i, 0) = diffN_MBQUAD0x(coords(i, 1));
75 diffN(i, 1) = diffN_MBQUAD0y(coords(i, 0));
76 diffN(i, 2) = diffN_MBQUAD1x(coords(i, 1));
77 diffN(i, 3) = diffN_MBQUAD1y(coords(i, 0));
78 diffN(i, 4) = diffN_MBQUAD2x(coords(i, 1));
79 diffN(i, 5) = diffN_MBQUAD2y(coords(i, 0));
80 diffN(i, 6) = diffN_MBQUAD3x(coords(i, 1));
81 diffN(i, 7) = diffN_MBQUAD3y(coords(i, 0));
82 }
83
84 /* BUBBLES */
85 {
86 std::array<int, 4> faces_nodes = {0, 1, 2, 3};
87 int p = 7;
88 int P = NBFACEQUAD_H1(p);
89 MatrixDouble quad_bubbles(verts.size(), P);
90 MatrixDouble quad_diff_bubbles(verts.size(), P * 2);
91 constexpr double eps = 1e-4;
92 double quad_bubbles_sum = 3.275491e+01;
93 double quad_diff_bubbles_sum = 5.838131e+02;
94
95 CHKERR H1_QuadShapeFunctions_MBQUAD(
96 faces_nodes.data(), p, &*N.data().begin(), &*diffN.data().begin(),
97 &*quad_bubbles.data().begin(), &*quad_diff_bubbles.data().begin(),
98 verts.size(), Legendre_polynomials);
99
100 if (std::abs(quad_bubbles_sum - sum_matrix(quad_bubbles)) > eps) {
101 SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
102 "wrong result = %8.6e", sum_matrix(quad_bubbles));
103 }
104 if (std::abs(quad_diff_bubbles_sum - sum_matrix(quad_diff_bubbles)) >
105 eps) {
106 SETERRQ1(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
107 "wrong result = %8.6e", sum_matrix(quad_diff_bubbles));
108 }
109 }
110
111 /* EDGES */
112 {
113 int sense[] = {1, -1, 1, -1};
114 int p[] = {3, 4, 5, 6};
115 int P[4];
116 int ee = 0;
117 for (; ee < 4; ee++) {
118 P[ee] = NBEDGE_H1(p[ee]);
119 }
120
121 double *quad_edges_ptr[4];
122 std::array<MatrixDouble, 4> quad_edges;
123 for (auto ee : {0, 1, 2, 3}) {
124 quad_edges[ee] = MatrixDouble(verts.size(), P[ee]);
125 quad_edges_ptr[ee] = &quad_edges[ee](0, 0);
126 }
127
128 double *quad_diff_edges_ptr[4];
129 std::array<MatrixDouble, 4> quad_diff_edges;
130 for (auto ee : {0, 1, 2, 3}) {
131 quad_diff_edges[ee] = MatrixDouble(verts.size(), P[ee] * 2);
132 quad_diff_edges_ptr[ee] = &quad_diff_edges[ee](0, 0);
133 }
134
135 double eps = 1e-4;
136 double quad_edges_sum[] = {1.237500e+01, 1.535050e+01, 1.781890e+01,
137 2.015678e+01};
138 double quad_diff_edges_sum[] = {8.470000e+01, 1.192510e+02, 1.580128e+02,
139 1.976378e+02};
140
141 CHKERR H1_EdgeShapeFunctions_MBQUAD(
142 sense, p, &*N.data().begin(), &*diffN.data().begin(), quad_edges_ptr,
143 quad_diff_edges_ptr, verts.size(), Legendre_polynomials);
144
145 for (int ee = 0; ee < 4; ++ee) {
146 if (std::abs(quad_edges_sum[ee] - sum_matrix(quad_edges[ee])) > eps)
147 SETERRQ2(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
148 "edge %d wrong result = %8.6e", ee,
149 sum_matrix(quad_edges[ee]));
150
151 if (std::abs(quad_diff_edges_sum[ee] -
152 sum_matrix(quad_diff_edges[ee])) > eps)
153 SETERRQ2(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
154 "edge %d wrong result = %8.6e", ee,
155 sum_matrix(quad_diff_edges[ee]));
156 }
157 }
158 }
159 CATCH_ERRORS;
160
161 CHKERR MoFEM::Core::Finalize();
162}
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:90
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:25
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:53
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
help
static char help[]
Definition: edge_and_bubble_shape_functions_on_quad.cpp:24
sum_matrix
static double sum_matrix(MatrixDouble &m)
Definition: edge_and_bubble_shape_functions_on_quad.cpp:26
diffN_MBQUAD2y
#define diffN_MBQUAD2y(x)
Definition: fem_tools.h:76
diffN_MBQUAD1x
#define diffN_MBQUAD1x(y)
Definition: fem_tools.h:73
N_MBQUAD3
#define N_MBQUAD3(x, y)
quad shape function
Definition: fem_tools.h:70
diffN_MBQUAD0x
#define diffN_MBQUAD0x(y)
Definition: fem_tools.h:71
diffN_MBQUAD1y
#define diffN_MBQUAD1y(x)
Definition: fem_tools.h:74
diffN_MBQUAD3y
#define diffN_MBQUAD3y(x)
Definition: fem_tools.h:78
diffN_MBQUAD0y
#define diffN_MBQUAD0y(x)
Definition: fem_tools.h:72
N_MBQUAD0
#define N_MBQUAD0(x, y)
quad shape function
Definition: fem_tools.h:67
diffN_MBQUAD3x
#define diffN_MBQUAD3x(y)
Definition: fem_tools.h:77
diffN_MBQUAD2x
#define diffN_MBQUAD2x(y)
Definition: fem_tools.h:75
N_MBQUAD2
#define N_MBQUAD2(x, y)
quad shape function
Definition: fem_tools.h:69
N_MBQUAD1
#define N_MBQUAD1(x, y)
quad shape function
Definition: fem_tools.h:68
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:25
NBFACEQUAD_H1
#define NBFACEQUAD_H1(P)
Number of base functions on quad for H1 space.
Definition: h1_hdiv_hcurl_l2.h:75
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:969
NBEDGE_H1
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
Definition: h1_hdiv_hcurl_l2.h:65
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:1101
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:39
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:35
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1096
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:183
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
N
const int N
Definition: speed_test.cpp:3
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:92
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:85
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125

◆ 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 26 of file edge_and_bubble_shape_functions_on_quad.cpp.

26 {
27 double s = 0;
28 for (unsigned int ii = 0; ii < m.size1(); ii++) {
29 for (unsigned int jj = 0; jj < m.size2(); jj++) {
30 s += std::abs(m(ii, jj));
31 }
32 }
33 return s;
34}
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:95

Variable Documentation

◆ help

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

Definition at line 24 of file edge_and_bubble_shape_functions_on_quad.cpp.

Generated on Tue Jul 5 2022 22:52:23 for MoFEM by Doxygen 1.9.4 and hosted at University of Glasgow