v0.15.0
Loading...
Searching...
No Matches
Classes | Public Member Functions | Public Attributes | Private Attributes | List of all members
MoFEM::DGSetIntegrationAtFrontFace Struct Reference
Collaboration diagram for MoFEM::DGSetIntegrationAtFrontFace:
[legend]

Classes

struct  Fe
 

Public Member Functions

 DGSetIntegrationAtFrontFace (boost::shared_ptr< Range > edges_ptr, ForcesAndSourcesCore::RuleHookFun fun_rule)
 
MoFEMErrorCode operator() (ForcesAndSourcesCore *fe_raw_ptr, int order_row, int order_col, int order_data)
 

Public Attributes

ForcesAndSourcesCore::RuleHookFun funRule
 

Private Attributes

boost::shared_ptr< Range > edgesPtr
 

Detailed Description

Definition at line 21 of file DGProjection.cpp.

Constructor & Destructor Documentation

◆ DGSetIntegrationAtFrontFace()

MoFEM::DGSetIntegrationAtFrontFace::DGSetIntegrationAtFrontFace ( boost::shared_ptr< Range >  edges_ptr,
ForcesAndSourcesCore::RuleHookFun  fun_rule 
)
inline

Definition at line 27 of file DGProjection.cpp.

29 : edgesPtr(edges_ptr), funRule(fun_rule) {};
MoFEM::DGSetIntegrationAtFrontFace::edgesPtr
boost::shared_ptr< Range > edgesPtr
Definition DGProjection.cpp:236
MoFEM::DGSetIntegrationAtFrontFace::funRule
ForcesAndSourcesCore::RuleHookFun funRule
Definition DGProjection.cpp:23

Member Function Documentation

◆ operator()()

MoFEMErrorCode MoFEM::DGSetIntegrationAtFrontFace::operator() ( ForcesAndSourcesCore *  fe_raw_ptr,
int  order_row,
int  order_col,
int  order_data 
)
inline

Definition at line 31 of file DGProjection.cpp.

32 {
33 MoFEMFunctionBegin;
34
35 constexpr bool debug = false;
36
37 constexpr int numNodes = 3;
38 constexpr int numEdges = 3;
39 constexpr int refinementLevels = 1;
40
41 auto &m_field = fe_raw_ptr->mField;
42 auto fe_ptr = static_cast<Fe *>(fe_raw_ptr);
43 auto fe_handle = fe_ptr->getFEEntityHandle();
44 auto type = type_from_handle(fe_handle);
45 if (type != MBTRI) {
46 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
47 "only MBTRI is implemented");
48 }
49
50
51
52 auto set_base_quadrature = [&]() {
53 MoFEMFunctionBegin;
54 int rule = funRule(order_row, order_col, order_data);
55 if (rule < QUAD_2D_TABLE_SIZE) {
56 if (QUAD_2D_TABLE[rule]->dim != 2) {
57 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
58 }
59 if (QUAD_2D_TABLE[rule]->order < rule) {
60 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
61 "wrong order %d != %d", QUAD_2D_TABLE[rule]->order, rule);
62 }
63 const size_t nb_gauss_pts = QUAD_2D_TABLE[rule]->npoints;
64 fe_ptr->gaussPts.resize(3, nb_gauss_pts, false);
65 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[1], 3,
66 &fe_ptr->gaussPts(0, 0), 1);
67 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[2], 3,
68 &fe_ptr->gaussPts(1, 0), 1);
69 cblas_dcopy(nb_gauss_pts, QUAD_2D_TABLE[rule]->weights, 1,
70 &fe_ptr->gaussPts(2, 0), 1);
71 auto &data = fe_ptr->dataOnElement[H1];
72 data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 3,
73 false);
74 double *shape_ptr =
75 &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
76 cblas_dcopy(3 * nb_gauss_pts, QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
77 1);
78 data->dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(3, 2, false);
79 std::copy(
80 Tools::diffShapeFunMBTRI.begin(), Tools::diffShapeFunMBTRI.end(),
81 data->dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).data().begin());
82
83 } else {
84 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
85 "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
86 }
87 MoFEMFunctionReturn(0);
88 };
89
90 CHKERR set_base_quadrature();
91
92 if (edgesPtr) {
93
94 auto get_edges = [&]() {
95 std::bitset<numEdges> edges;
96 for (int ee = 0; ee != numEdges; ee++) {
97 EntityHandle edge;
98 CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
99 if (edgesPtr->find(edge) != edgesPtr->end()) {
100 edges.set(ee);
101 } else {
102 edges.reset(ee);
103 }
104 }
105 return edges;
106 };
107
108 auto set_gauss_pts = [&](auto &ref_gauss_pts) {
109 MoFEMFunctionBegin;
110 fe_ptr->gaussPts.swap(ref_gauss_pts);
111 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
112 auto &data = fe_ptr->dataOnElement[H1];
113 data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4);
114 auto *shape_ptr =
115 data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().data();
116 CHKERR Tools::shapeFunMBTRI(shape_ptr, &fe_ptr->gaussPts(0, 0),
117 &fe_ptr->gaussPts(1, 0), nb_gauss_pts);
118 MoFEMFunctionReturn(0);
119 };
120
121 auto refine_quadrature = [&]() {
122 MoFEMFunctionBegin;
123
124 const int max_level = refinementLevels;
125
126 moab::Core moab_ref;
127 double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0};
128 EntityHandle nodes[numNodes];
129 for (int nn = 0; nn != numNodes; nn++)
130 CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
131 EntityHandle tri;
132 CHKERR moab_ref.create_element(MBTRI, nodes, numNodes, tri);
133 MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
134 MoFEM::Interface &m_field_ref = mofem_ref_core;
135 {
136 Range tris(tri, tri);
137 Range edges;
138 CHKERR m_field_ref.get_moab().get_adjacencies(tris, 1, true, edges,
139 moab::Interface::UNION);
140 CHKERR m_field_ref.getInterface<BitRefManager>()->setBitRefLevel(
141 tris, BitRefLevel().set(0), false, VERBOSE);
142 }
143
144 auto edges = get_edges();
145
146 EntityHandle meshset;
147 CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
148 for (int ee = 0; ee != numEdges; ee++) {
149 if (edges[ee]) {
150 EntityHandle ent;
151 CHKERR moab_ref.side_element(tri, 1, ee, ent);
152 CHKERR moab_ref.add_entities(meshset, &ent, 1);
153 }
154 }
155
156 // refine mesh
157 auto *m_ref = m_field_ref.getInterface<MeshRefinement>();
158 for (int ll = 0; ll != max_level; ll++) {
159 Range ref_edges;
160 CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ref_edges,
161 true);
162 CHKERR m_field_ref.getInterface<BitRefManager>()
163 ->filterEntitiesByRefLevel(BitRefLevel().set(ll),
164 BitRefLevel().set(), ref_edges);
165
166 Range tris;
167 CHKERR m_field_ref.getInterface<BitRefManager>()
168 ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
169 BitRefLevel().set(), MBTRI, tris);
170 CHKERR m_ref->addVerticesInTheMiddleOfEdges(
171 ref_edges, BitRefLevel().set(ll + 1));
172 CHKERR m_ref->refineTris(tris, BitRefLevel().set(ll + 1));
173 CHKERR m_field_ref.getInterface<BitRefManager>()
174 ->updateMeshsetByEntitiesChildren(
175 meshset, BitRefLevel().set(ll + 1), meshset, MBEDGE, true);
176 }
177
178 // get ref coords
179 Range tris;
180 CHKERR m_field_ref.getInterface<BitRefManager>()
181 ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
182 BitRefLevel().set(), MBTRI, tris);
183
184 if (debug) {
185 CHKERR dg_save_range(moab_ref, "ref_tris.vtk", tris);
186 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "debug");
187 }
188
189 MatrixDouble ref_coords(tris.size(), 9, false);
190 int tt = 0;
191 for (Range::iterator tit = tris.begin(); tit != tris.end();
192 tit++, tt++) {
193 int num_nodes;
194 const EntityHandle *conn;
195 CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
196 CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
197 }
198
199 auto &data = fe_ptr->dataOnElement[H1];
200 const size_t nb_gauss_pts = fe_ptr->gaussPts.size2();
201 MatrixDouble ref_gauss_pts(3, nb_gauss_pts * ref_coords.size1());
202 MatrixDouble &shape_n = data->dataOnEntities[MBVERTEX][0].getN(NOBASE);
203 int gg = 0;
204 for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
205 double *tri_coords = &ref_coords(tt, 0);
206 FTensor::Tensor1<double, 3> t_normal;
207 CHKERR Tools::getTriNormal(tri_coords, &t_normal(0));
208 auto det = t_normal.l2();
209 for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
210 for (int dd = 0; dd != 2; dd++) {
211 ref_gauss_pts(dd, gg) = shape_n(ggg, 0) * tri_coords[3 * 0 + dd] +
212 shape_n(ggg, 1) * tri_coords[3 * 1 + dd] +
213 shape_n(ggg, 2) * tri_coords[3 * 2 + dd];
214 }
215 ref_gauss_pts(2, gg) = fe_ptr->gaussPts(2, ggg) * det;
216 }
217 }
218
219 MoFEMFunctionReturn(0);
220 };
221
222 CHKERR refine_quadrature();
223 }
224
225 MoFEMFunctionReturn(0);
226 }
VERBOSE
@ VERBOSE
Definition definitions.h:222
NOBASE
@ NOBASE
Definition definitions.h:59
H1
@ H1
continuous field
Definition definitions.h:85
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
order
constexpr int order
Order displacement.
Definition dg_projection.cpp:38
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition ddTensor0.hpp:33
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition Templates.hpp:1992
MoFEM::debug
static const bool debug
Definition ConstrainMatrixCtx.cpp:9
MoFEM::dg_save_range
static auto dg_save_range(moab::Interface &moab, const std::string name, const Range r, std::vector< Tag > tags={})
Definition DGProjection.cpp:7
QUAD_2D_TABLE_SIZE
#define QUAD_2D_TABLE_SIZE
Definition quad.h:174
QUAD_3D_TABLE_SIZE
#define QUAD_3D_TABLE_SIZE
Definition quad.h:186
QUAD_2D_TABLE
static QUAD *const QUAD_2D_TABLE[]
Definition quad.h:175
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::Tools::diffShapeFunMBTRI
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition Tools.hpp:104
MoFEM::Tools::getTriNormal
static MoFEMErrorCode getTriNormal(const double *coords, double *normal, double *d_normal=nullptr)
Get the Tri Normal objectGet triangle normal.
Definition Tools.cpp:353
MoFEM::Tools::shapeFunMBTRI
static MoFEMErrorCode shapeFunMBTRI(double *shape, const double *ksi, const double *eta, const int nb)
Calculate shape functions on triangle.
Definition Tools.hpp:728
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93
QUAD_::npoints
int npoints
Definition quad.h:29

Member Data Documentation

◆ edgesPtr

boost::shared_ptr<Range> MoFEM::DGSetIntegrationAtFrontFace::edgesPtr
private

Definition at line 236 of file DGProjection.cpp.

◆ funRule

ForcesAndSourcesCore::RuleHookFun MoFEM::DGSetIntegrationAtFrontFace::funRule
Initial value:
= [](int, int, int p) {
return 2 * (p + 1);
}

Definition at line 23 of file DGProjection.cpp.

23 {
24 return 2 * (p + 1);
25 };
The documentation for this struct was generated from the following file:
Generated by Doxygen 1.9.8 and hosted at University of Glasgow