v0.13.2
Loading...
Searching...
No Matches
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[] 
)

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}
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:90
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
help
static char help[]
Definition: edge_and_bubble_shape_functions_on_quad.cpp:10
sum_matrix
static double sum_matrix(MatrixDouble &m)
Definition: edge_and_bubble_shape_functions_on_quad.cpp:12
diffN_MBQUAD2y
#define diffN_MBQUAD2y(x)
Definition: fem_tools.h:66
diffN_MBQUAD1x
#define diffN_MBQUAD1x(y)
Definition: fem_tools.h:63
N_MBQUAD3
#define N_MBQUAD3(x, y)
quad shape function
Definition: fem_tools.h:60
diffN_MBQUAD0x
#define diffN_MBQUAD0x(y)
Definition: fem_tools.h:61
diffN_MBQUAD1y
#define diffN_MBQUAD1y(x)
Definition: fem_tools.h:64
diffN_MBQUAD3y
#define diffN_MBQUAD3y(x)
Definition: fem_tools.h:68
diffN_MBQUAD0y
#define diffN_MBQUAD0y(x)
Definition: fem_tools.h:62
N_MBQUAD0
#define N_MBQUAD0(x, y)
quad shape function
Definition: fem_tools.h:57
diffN_MBQUAD3x
#define diffN_MBQUAD3x(y)
Definition: fem_tools.h:67
diffN_MBQUAD2x
#define diffN_MBQUAD2x(y)
Definition: fem_tools.h:65
N_MBQUAD2
#define N_MBQUAD2(x, y)
quad shape function
Definition: fem_tools.h:59
N_MBQUAD1
#define N_MBQUAD1(x, y)
quad shape function
Definition: fem_tools.h:58
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
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
NBEDGE_H1
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
Definition: h1_hdiv_hcurl_l2.h:55
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
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
N
const int N
Definition: speed_test.cpp:3
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
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
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112

◆ 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}
m
FTensor::Index< 'm', SPACE_DIM > m
Definition: fluid_structure_eigenproblem.cpp:65

Variable Documentation

◆ help

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

Definition at line 10 of file edge_and_bubble_shape_functions_on_quad.cpp.

Generated by Doxygen 1.9.5 and hosted at University of Glasgow