v0.14.0
gauss_points_on_outer_product.cpp
Go to the documentation of this file.
1 /** \file gauss_points_on_outer_product.cpp
2  \example gauss_points_on_outer_product.cpp
3  \brief Testing gauss points coordinates and weights for a quad face
4 
5 */
6 
7 #include <MoFEM.hpp>
8 
9 #ifdef __cplusplus
10 extern "C" {
11 #endif
12 #include <cblas.h>
13 #include <quad.h>
14 #ifdef __cplusplus
15 }
16 #endif
17 
18 using namespace MoFEM;
19 static char help[] = "...\n\n";
20 
21 int main(int argc, char *argv[]) {
22 
23  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
24  try {
25  moab::Core mb_instance;
26  moab::Interface &moab = mb_instance;
27  MoFEM::Core core(moab);
28 
29  auto test_quad = [&]() {
31  MatrixDouble pts_quad;
32  int rule_ksi = 6;
33  int rule_eta = 8;
35  rule_eta);
36  int nb_gauss_pts = pts_quad.size2();
37  double sum_coords = 0, sum_gauss_pts = 0;
38  for (int i = 0; i < nb_gauss_pts; ++i) {
39  sum_coords += pts_quad(0, i) + pts_quad(1, i);
40  sum_gauss_pts += pts_quad(2, i);
41  }
42  double eps = 1e-8;
43  if (fabs(20.0 - sum_coords) > eps) {
44  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
45  "wrong result %3.4e", sum_coords);
46  }
47  if (fabs(1.0 - sum_gauss_pts) > eps) {
48  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
49  "wrong result %3.4e", sum_gauss_pts);
50  }
52  };
53 
54  auto test_hex = [&]() {
56  MatrixDouble pts_quad;
57  int rule_ksi = 4;
58  int rule_eta = 5;
59  int rule_zet = 6;
61  rule_eta, rule_zet);
62  int nb_gauss_pts = pts_quad.size2();
63  double sum_coords = 0, sum_gauss_pts = 0;
64  for (int i = 0; i < nb_gauss_pts; ++i) {
65  sum_coords += pts_quad(0, i) + pts_quad(1, i) + pts_quad(2, i);
66  sum_gauss_pts += pts_quad(3, i);
67  }
68  double eps = 1e-8;
69  if (fabs(54.0 - sum_coords) > eps) {
70  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
71  "wrong result %3.4e", sum_coords);
72  }
73  if (fabs(1.0 - sum_gauss_pts) > eps) {
74  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
75  "wrong result %3.4e", sum_gauss_pts);
76  }
78  };
79 
80  auto test_refined_triangle = [&]() {
82 
83  auto refine = Tools::refineTriangle(2);
84  auto new_gauss_pts = Tools::refineTriangleIntegrationPts(4, refine);
85 
86  int new_nb_gauss_pts = new_gauss_pts.size2();
87  double sum_coords = 0, sum_gauss_pts = 0;
88  for (int i = 0; i < new_nb_gauss_pts; ++i) {
89  sum_coords += new_gauss_pts(0, i) + new_gauss_pts(1, i);
90  sum_gauss_pts += new_gauss_pts(2, i);
91  }
92  constexpr double eps = 1e-8;
93 
94  MOFEM_LOG("WORLD", Sev::verbose)
95  << "sum_gauss_pts " << sum_gauss_pts << endl;
96  MOFEM_LOG("WORLD", Sev::verbose) << "sum_coords " << sum_coords << endl;
97  if (fabs(64.0 - sum_coords) > eps) {
98  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
99  "wrong result %3.4e", sum_coords);
100  }
101  if (fabs(1.0 - sum_gauss_pts) > eps) {
102  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
103  "wrong result %3.4e", sum_gauss_pts);
104  }
105 
106  constexpr bool debug = false;
107 
108  if (debug) {
109 
110  auto [nodes, triangles, level_index] = refine;
111 
112  auto get_coords = [&](auto t) {
113  auto [nodes, triangles, level_index] = refine;
114  std::array<double, 9> ele_coords;
115  for (auto n : {0, 1, 2}) {
116  for (auto i : {0, 1}) {
117  ele_coords[3 * n + i] = nodes[2 * triangles[6 * t + n] + i];
118  }
119  ele_coords[3 * n + 2] = 0;
120  }
121  return ele_coords;
122  };
123 
124  moab::Core mb_instance;
125  moab::Interface &moab = mb_instance;
126 
127  for (auto t = level_index[level_index.size() - 2];
128  t != level_index.back(); ++t) {
129  std::array<EntityHandle, 3> conn;
130  auto ele_coords = get_coords(t);
131 
132  for (auto n : {0, 1, 2}) {
133  CHKERR moab.create_vertex(&ele_coords[3 * n], conn[n]);
134  MOFEM_LOG("WORLD", Sev::noisy)
135  << ele_coords[3 * n] << " " << ele_coords[3 * n + 1] << " "
136  << ele_coords[3 * n + 2];
137  }
138 
139  MOFEM_LOG("WORLD", Sev::noisy)
140  << triangles[6 * t + 0] << " " << triangles[6 * t + 1] << " "
141  << triangles[6 * t + 2] << std::endl;
142 
143  EntityHandle tri;
144  CHKERR moab.create_element(MBTRI, conn.data(), 3, tri);
145  }
146 
147  for (int i = 0; i < new_nb_gauss_pts; ++i) {
148  std::array<double, 3> coords{new_gauss_pts(0, i), new_gauss_pts(1, i),
149  0};
150  EntityHandle node;
151  CHKERR moab.create_vertex(coords.data(), node);
152  }
153 
154  CHKERR moab.write_file("gauss_refine.vtk", "VTK", "");
155  }
156 
158  };
159 
160  CHKERR test_quad();
161  CHKERR test_hex();
162  CHKERR test_refined_triangle();
163  }
164  CATCH_ERRORS;
165 
167 }
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
EntityHandle
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM.hpp
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::Tools::outerProductOfEdgeIntegrationPtsForHex
static MoFEMErrorCode outerProductOfEdgeIntegrationPtsForHex(MatrixDouble &pts, const int edge0, const int edge1, const int edge2)
Definition: Tools.cpp:659
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::Tools::refineTriangleIntegrationPts
static MatrixDouble refineTriangleIntegrationPts(MatrixDouble pts, RefineTrianglesReturn refined)
generate integration points for refined triangle mesh for last level
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
t
constexpr double t
plate stiffness
Definition: plate.cpp:59
debug
static const bool debug
Definition: dm_create_subdm.cpp:12
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
convert.n
n
Definition: convert.py:82
help
static char help[]
Definition: gauss_points_on_outer_product.cpp:19
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_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
MoFEM::Tools::refineTriangle
static RefineTrianglesReturn refineTriangle(int nb_levels)
create uniform triangle mesh of refined elements
Definition: Tools.cpp:721
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::Tools::outerProductOfEdgeIntegrationPtsForQuad
static MoFEMErrorCode outerProductOfEdgeIntegrationPtsForQuad(MatrixDouble &pts, const int edge0, const int edge1)
Definition: Tools.cpp:610
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
main
int main(int argc, char *argv[])
Definition: gauss_points_on_outer_product.cpp:21
quad.h