v0.14.0
gauss_points_on_outer_product.cpp

Testing gauss points coordinates and weights for a quad face

/** \file gauss_points_on_outer_product.cpp
\example gauss_points_on_outer_product.cpp
\brief Testing gauss points coordinates and weights for a quad face
*/
#include <MoFEM.hpp>
#ifdef __cplusplus
extern "C" {
#endif
#include <cblas.h>
#include <quad.h>
#ifdef __cplusplus
}
#endif
using namespace MoFEM;
static char help[] = "...\n\n";
int main(int argc, char *argv[]) {
MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
try {
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
MoFEM::Core core(moab);
auto test_quad = [&]() {
MatrixDouble pts_quad;
int rule_ksi = 6;
int rule_eta = 8;
rule_eta);
int nb_gauss_pts = pts_quad.size2();
double sum_coords = 0, sum_gauss_pts = 0;
for (int i = 0; i < nb_gauss_pts; ++i) {
sum_coords += pts_quad(0, i) + pts_quad(1, i);
sum_gauss_pts += pts_quad(2, i);
}
double eps = 1e-8;
if (fabs(20.0 - sum_coords) > eps) {
SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"wrong result %3.4e", sum_coords);
}
if (fabs(1.0 - sum_gauss_pts) > eps) {
SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"wrong result %3.4e", sum_gauss_pts);
}
};
auto test_hex = [&]() {
MatrixDouble pts_quad;
int rule_ksi = 4;
int rule_eta = 5;
int rule_zet = 6;
rule_eta, rule_zet);
int nb_gauss_pts = pts_quad.size2();
double sum_coords = 0, sum_gauss_pts = 0;
for (int i = 0; i < nb_gauss_pts; ++i) {
sum_coords += pts_quad(0, i) + pts_quad(1, i) + pts_quad(2, i);
sum_gauss_pts += pts_quad(3, i);
}
double eps = 1e-8;
if (fabs(54.0 - sum_coords) > eps) {
SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"wrong result %3.4e", sum_coords);
}
if (fabs(1.0 - sum_gauss_pts) > eps) {
SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"wrong result %3.4e", sum_gauss_pts);
}
};
auto test_refined_triangle = [&]() {
auto refine = Tools::refineTriangle(2);
auto new_gauss_pts = Tools::refineTriangleIntegrationPts(4, refine);
int new_nb_gauss_pts = new_gauss_pts.size2();
double sum_coords = 0, sum_gauss_pts = 0;
for (int i = 0; i < new_nb_gauss_pts; ++i) {
sum_coords += new_gauss_pts(0, i) + new_gauss_pts(1, i);
sum_gauss_pts += new_gauss_pts(2, i);
}
constexpr double eps = 1e-8;
MOFEM_LOG("WORLD", Sev::verbose)
<< "sum_gauss_pts " << sum_gauss_pts << endl;
MOFEM_LOG("WORLD", Sev::verbose) << "sum_coords " << sum_coords << endl;
if (fabs(64.0 - sum_coords) > eps) {
SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"wrong result %3.4e", sum_coords);
}
if (fabs(1.0 - sum_gauss_pts) > eps) {
SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"wrong result %3.4e", sum_gauss_pts);
}
constexpr bool debug = false;
if (debug) {
auto [nodes, triangles, level_index] = refine;
auto get_coords = [&](auto t) {
auto [nodes, triangles, level_index] = refine;
std::array<double, 9> ele_coords;
for (auto n : {0, 1, 2}) {
for (auto i : {0, 1}) {
ele_coords[3 * n + i] = nodes[2 * triangles[6 * t + n] + i];
}
ele_coords[3 * n + 2] = 0;
}
return ele_coords;
};
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
for (auto t = level_index[level_index.size() - 2];
t != level_index.back(); ++t) {
std::array<EntityHandle, 3> conn;
auto ele_coords = get_coords(t);
for (auto n : {0, 1, 2}) {
CHKERR moab.create_vertex(&ele_coords[3 * n], conn[n]);
MOFEM_LOG("WORLD", Sev::noisy)
<< ele_coords[3 * n] << " " << ele_coords[3 * n + 1] << " "
<< ele_coords[3 * n + 2];
}
MOFEM_LOG("WORLD", Sev::noisy)
<< triangles[6 * t + 0] << " " << triangles[6 * t + 1] << " "
<< triangles[6 * t + 2] << std::endl;
CHKERR moab.create_element(MBTRI, conn.data(), 3, tri);
}
for (int i = 0; i < new_nb_gauss_pts; ++i) {
std::array<double, 3> coords{new_gauss_pts(0, i), new_gauss_pts(1, i),
0};
CHKERR moab.create_vertex(coords.data(), node);
}
CHKERR moab.write_file("gauss_refine.vtk", "VTK", "");
}
};
CHKERR test_quad();
CHKERR test_hex();
CHKERR test_refined_triangle();
}
}
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:1094
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