v0.13.2
Loading...
Searching...
No Matches
FaceElementForcesAndSourcesCoreOnParent.cpp
Go to the documentation of this file.
1/** \file FaceElementForcesAndSourcesCoreOnParent.cpp
2
3\brief Implementation of face element
4
5*/
6
7
8
9namespace MoFEM {
10
12 return -1;
13};
14
18
19 auto ref_fe = refinePtrFE;
20 if (ref_fe == nullptr)
21 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
22 "Pointer to face element is not set");
23
24 const auto ref_ent = ref_fe->getFEEntityHandle();
25 const auto ent = getFEEntityHandle();
26
27 const auto &ref_gauss_pts = ref_fe->gaussPts;
28 const auto nb_integration_points = ref_gauss_pts.size2();
29
30 auto set_integration_pts_for_tri = [&]() {
32 auto get_coords = [&](const auto ent) {
33 int num_nodes;
34 const EntityHandle *conn;
35 CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
36 std::array<double, 9> node_coords;
37 CHKERR mField.get_moab().get_coords(conn, num_nodes, node_coords.data());
38 return node_coords;
39 };
40
41 auto ref_node_coords = get_coords(ref_ent);
42 auto node_coords = get_coords(ent);
43
44 MatrixDouble ref_shapes(nb_integration_points, 3);
45 CHKERR Tools::shapeFunMBTRI<1>(&ref_shapes(0, 0), &ref_gauss_pts(0, 0),
46 &ref_gauss_pts(1, 0), nb_integration_points);
47
48#ifndef NDEBUG
49 if (ref_shapes.size1() * ref_shapes.size2() !=
50 nb_integration_points * ref_node_coords.size() / 3)
51 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
52 "Wrong number of shape functions");
53#endif
54
55 auto get_glob_coords = [&]() {
56 MatrixDouble glob_coords(nb_integration_points, 3);
57
59 &glob_coords(0, 0), &glob_coords(0, 1), &glob_coords(0, 2)};
60 auto t_shape =
62
63 for (auto gg = 0; gg != nb_integration_points; ++gg) {
65
66 auto t_ref_node_coords =
68 &ref_node_coords[0], &ref_node_coords[1], &ref_node_coords[2]};
69
70 t_glob_coords(i) = 0;
71 for (int nn = 0; nn != ref_node_coords.size() / 3; ++nn) {
72 t_glob_coords(i) += t_shape * t_ref_node_coords(i);
73 ++t_ref_node_coords;
74 ++t_shape;
75 }
76
77 ++t_glob_coords;
78 }
79 return glob_coords;
80 };
81
82 auto glob_coords = get_glob_coords();
83 MatrixDouble local_coords(nb_integration_points, 2);
84
86 node_coords.data(), &glob_coords(0, 0), nb_integration_points,
87 &local_coords(0, 0));
88
89 gaussPts.resize(3, nb_integration_points, false);
91 &gaussPts(0, 0), &gaussPts(1, 0), &gaussPts(2, 0)};
92 auto t_local_coords = FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2>{
93 &local_coords(0, 0), &local_coords(0, 1)};
94
95 auto scale_quadarture = [&]() {
97 FTensor::Tensor1<double, 3> t_ref_normal;
98 CHKERR Tools::getTriNormal(node_coords.data(), &t_normal(0));
99 CHKERR Tools::getTriNormal(ref_node_coords.data(), &t_ref_normal(0));
101 const double scale = std::sqrt((t_ref_normal(J) * t_ref_normal(J)) /
102 (t_normal(J) * t_normal(J)));
103 return scale;
104 };
105
106 const auto sq = scale_quadarture();
107
108 for (auto gg = 0; gg != nb_integration_points; ++gg) {
110 t_gauss_pts(i) = t_local_coords(i);
111 t_gauss_pts(2) = sq * ref_gauss_pts(2, gg);
112 ++t_gauss_pts;
113 ++t_local_coords;
114 }
116 };
117
118 const auto type = numeredEntFiniteElementPtr->getEntType();
119
120 switch (type) {
121 case MBTRI:
122 CHKERR set_integration_pts_for_tri();
123 break;
124 default:
125 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
126 "Element type not implemented: %d", type);
127 }
128
130}
131
132} // namespace MoFEM
static Index< 'J', 3 > J
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
FTensor::Index< 'i', SPACE_DIM > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
double scale
Definition: plastic.cpp:97
virtual moab::Interface & get_moab()=0
EntityHandle getFEEntityHandle() const
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
MatrixDouble gaussPts
Matrix of integration points.
ForcesAndSourcesCore * refinePtrFE
Element to integrate parent or child.
static MoFEMErrorCode getLocalCoordinatesOnReferenceTriNodeTri(const double *elem_coords, const double *glob_coords, const int nb_nodes, double *local_coords)
Get the Local Coordinates On Reference Four Node Tet object.
Definition: Tools.cpp:142
static MoFEMErrorCode getTriNormal(const double *coords, double *normal)
Get the Tri Normal objectGet triangle normal.
Definition: Tools.cpp:350