v0.13.1
edge_and_bubble_shape_functions_on_quad.cpp
Go to the documentation of this file.
1/** \file edge_and_bubble_shape_functions_on_quad.cpp
2 \example edge_and_bubble_shape_functions_on_quad.cpp
3 \brief Testing edge and bubble shape functions on a quad face
4
5*/
6
7/* This file is part of MoFEM.
8 * MoFEM is free software: you can redistribute it and/or modify it under
9 * the terms of the GNU Lesser General Public License as published by the
10 * Free Software Foundation, either version 3 of the License, or (at your
11 * option) any later version.
12 *
13 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 * License for more details.
17 *
18 * You should have received a copy of the GNU Lesser General Public
19 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20
21#include <MoFEM.hpp>
22
23using namespace MoFEM;
24static char help[] = "...\n\n";
25
26static double sum_matrix(MatrixDouble &m) {
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}
35
36int main(int argc, char *argv[]) {
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
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
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 }
160
162}
static Index< 'p', 3 > p
static const double eps
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:53
#define CHKERR
Inline error check.
Definition: definitions.h:548
int main(int argc, char *argv[])
static double sum_matrix(MatrixDouble &m)
#define diffN_MBQUAD2y(x)
Definition: fem_tools.h:76
#define diffN_MBQUAD1x(y)
Definition: fem_tools.h:73
#define N_MBQUAD3(x, y)
quad shape function
Definition: fem_tools.h:70
#define diffN_MBQUAD0x(y)
Definition: fem_tools.h:71
#define diffN_MBQUAD1y(x)
Definition: fem_tools.h:74
#define diffN_MBQUAD3y(x)
Definition: fem_tools.h:78
#define diffN_MBQUAD0y(x)
Definition: fem_tools.h:72
#define N_MBQUAD0(x, y)
quad shape function
Definition: fem_tools.h:67
#define diffN_MBQUAD3x(y)
Definition: fem_tools.h:77
#define diffN_MBQUAD2x(y)
Definition: fem_tools.h:75
#define N_MBQUAD2(x, y)
quad shape function
Definition: fem_tools.h:69
#define N_MBQUAD1(x, y)
quad shape function
Definition: fem_tools.h:68
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Legendre approximation basis.
#define NBFACEQUAD_H1(P)
Number of base functions on quad for H1 space.
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
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
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
FTensor::Index< 'i', SPACE_DIM > i
char mesh_file_name[255]
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
CoreTmp< 0 > Core
Definition: Core.hpp:1096
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
FTensor::Index< 'm', 3 > m
const int N
Definition: speed_test.cpp:3
Core (interface) class.
Definition: Core.hpp:92
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125