v0.13.2
Loading...
Searching...
No Matches
FaceElementForcesAndSourcesCore.cpp
Go to the documentation of this file.
1/** \file FaceElementForcesAndSourcesCore.cpp
2
3\brief Implementation of face element
4
5*/
6
7namespace MoFEM {
8
10 Interface &m_field)
11 : ForcesAndSourcesCore(m_field),
12 meshPositionsFieldName("MESH_NODE_POSITIONS") {}
13
17
18 auto type = numeredEntFiniteElementPtr->getEntType();
19
23
24 auto get_ftensor_from_vec_3d = [](VectorDouble &v) {
26 &v[2]);
27 };
28
29 auto get_ftensor_n_diff = [&]() {
30 const auto &m = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
32 &m(0, 1));
33 };
34
35 auto get_ftensor_from_mat_3d = [](MatrixDouble &m) {
37 &m(0, 0), &m(0, 1), &m(0, 2));
38 };
39
40 if (type == MBTRI) {
41
42 const size_t nb_gauss_pts = gaussPts.size2();
43 normalsAtGaussPts.resize(nb_gauss_pts, 3);
44 tangentOneAtGaussPts.resize(nb_gauss_pts, 3);
45 tangentTwoAtGaussPts.resize(nb_gauss_pts, 3);
46
47 auto t_tan1 = get_ftensor_from_mat_3d(tangentOneAtGaussPts);
48 auto t_tan2 = get_ftensor_from_mat_3d(tangentTwoAtGaussPts);
49 auto t_normal = get_ftensor_from_mat_3d(normalsAtGaussPts);
50
51 auto t_n =
54 &tangentOne[2]);
56 &tangentTwo[2]);
57
58 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
59 t_normal(i) = t_n(i);
60 t_tan1(i) = t_t1(i);
61 t_tan2(i) = t_t2(i);
62 ++t_tan1;
63 ++t_tan2;
64 ++t_normal;
65 }
66
67 } else if (type == MBQUAD) {
68
70 CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
71 coords.resize(num_nodes * 3, false);
72 CHKERR mField.get_moab().get_coords(conn, num_nodes,
73 &*coords.data().begin());
74
75 const size_t nb_gauss_pts = gaussPts.size2();
76 normalsAtGaussPts.resize(nb_gauss_pts, 3);
77 tangentOneAtGaussPts.resize(nb_gauss_pts, 3);
78 tangentTwoAtGaussPts.resize(nb_gauss_pts, 3);
79 normalsAtGaussPts.clear();
82
83 auto t_t1 = get_ftensor_from_mat_3d(tangentOneAtGaussPts);
84 auto t_t2 = get_ftensor_from_mat_3d(tangentTwoAtGaussPts);
85 auto t_normal = get_ftensor_from_mat_3d(normalsAtGaussPts);
86
89
90 auto t_diff = get_ftensor_n_diff();
91 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
92 auto t_coords = get_ftensor_from_vec_3d(coords);
93 for (int nn = 0; nn != num_nodes; ++nn) {
94 t_t1(i) += t_coords(i) * t_diff(N0);
95 t_t2(i) += t_coords(i) * t_diff(N1);
96 ++t_diff;
97 ++t_coords;
98 }
99 t_normal(j) = FTensor::levi_civita(i, j, k) * t_t1(k) * t_t2(i);
100
101 ++t_t1;
102 ++t_t2;
103 ++t_normal;
104 }
105 } else {
106 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
107 "Element type not implemented");
108 }
109
111}
112
115
117 CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
118 coords.resize(num_nodes * 3, false);
119 CHKERR mField.get_moab().get_coords(conn, num_nodes, &*coords.data().begin());
120 nOrmal.resize(3, false);
121 tangentOne.resize(3, false);
122 tangentTwo.resize(3, false);
123
124 auto calc_normal = [&](const double *diff_ptr) {
127 &coords[0], &coords[1], &coords[2]);
129 &nOrmal[0], &nOrmal[1], &nOrmal[2]);
131 &tangentOne[0], &tangentOne[1], &tangentOne[2]);
133 &tangentTwo[0], &tangentTwo[1], &tangentTwo[2]);
135 diff_ptr, &diff_ptr[1]);
136
140
143 t_t1(i) = 0;
144 t_t2(i) = 0;
145
146 for (int nn = 0; nn != num_nodes; ++nn) {
147 t_t1(i) += t_coords(i) * t_diff(N0);
148 t_t2(i) += t_coords(i) * t_diff(N1);
149 ++t_coords;
150 ++t_diff;
151 }
152 t_normal(j) = FTensor::levi_civita(i, j, k) * t_t1(k) * t_t2(i);
153 aRea = sqrt(t_normal(i) * t_normal(i));
155 };
156
157 const double *diff_ptr;
158 switch (numeredEntFiniteElementPtr->getEntType()) {
159 case MBTRI:
160 diff_ptr = Tools::diffShapeFunMBTRI.data();
161 CHKERR calc_normal(diff_ptr);
162 // FIXME: Normal should be divided not the area for triangle!!
163 aRea /= 2;
164 break;
165 case MBQUAD:
166 diff_ptr = Tools::diffShapeFunMBQUADAtCenter.data();
167 CHKERR calc_normal(diff_ptr);
168 break;
169 default:
170 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
171 "Element type not implemented");
172 }
173
175}
176
179 // Set integration points
180 int order_data = getMaxDataOrder();
181 int order_row = getMaxRowOrder();
182 int order_col = getMaxColOrder();
183 int rule = getRule(order_row, order_col, order_data);
184
185 auto set_integration_pts_for_tri = [&]() {
187 if (rule < QUAD_2D_TABLE_SIZE) {
188 if (QUAD_2D_TABLE[rule]->dim != 2) {
189 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
190 }
191 if (QUAD_2D_TABLE[rule]->order < rule) {
192 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
193 "wrong order %d != %d", QUAD_2D_TABLE[rule]->order, rule);
194 }
195 const size_t nb_gauss_pts = QUAD_2D_TABLE[rule]->npoints;
196 gaussPts.resize(3, nb_gauss_pts, false);
197 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[1], 3,
198 &gaussPts(0, 0), 1);
199 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[2], 3,
200 &gaussPts(1, 0), 1);
201 cblas_dcopy(nb_gauss_pts, QUAD_2D_TABLE[rule]->weights, 1,
202 &gaussPts(2, 0), 1);
203 dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 3,
204 false);
205 double *shape_ptr =
206 &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
207 cblas_dcopy(3 * nb_gauss_pts, QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
208 1);
209 dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(3, 2, false);
210 std::copy(
212 dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).data().begin());
213 } else {
214 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
215 "rule > quadrature order %d < %d", rule, QUAD_2D_TABLE_SIZE);
216 }
218 };
219
220 auto calc_base_for_tri = [&]() {
222 const size_t nb_gauss_pts = gaussPts.size2();
223 auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
224 auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
225 base.resize(nb_gauss_pts, 3, false);
226 diff_base.resize(3, 2, false);
227 CHKERR ShapeMBTRI(&*base.data().begin(), &gaussPts(0, 0), &gaussPts(1, 0),
228 nb_gauss_pts);
229 std::copy(
231 dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).data().begin());
233 };
234
235 auto calc_base_for_quad = [&]() {
237 const size_t nb_gauss_pts = gaussPts.size2();
238 auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
239 auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
240 base.resize(nb_gauss_pts, 4, false);
241 diff_base.resize(nb_gauss_pts, 8, false);
242 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
243 const double ksi = gaussPts(0, gg);
244 const double zeta = gaussPts(1, gg);
245 base(gg, 0) = N_MBQUAD0(ksi, zeta);
246 base(gg, 1) = N_MBQUAD1(ksi, zeta);
247 base(gg, 2) = N_MBQUAD2(ksi, zeta);
248 base(gg, 3) = N_MBQUAD3(ksi, zeta);
249 diff_base(gg, 0) = diffN_MBQUAD0x(zeta);
250 diff_base(gg, 1) = diffN_MBQUAD0y(ksi);
251 diff_base(gg, 2) = diffN_MBQUAD1x(zeta);
252 diff_base(gg, 3) = diffN_MBQUAD1y(ksi);
253 diff_base(gg, 4) = diffN_MBQUAD2x(zeta);
254 diff_base(gg, 5) = diffN_MBQUAD2y(ksi);
255 diff_base(gg, 6) = diffN_MBQUAD3x(zeta);
256 diff_base(gg, 7) = diffN_MBQUAD3y(ksi);
257 }
259 };
260
261 const auto type = numeredEntFiniteElementPtr->getEntType();
262
263 if (rule >= 0) {
264 switch (type) {
265 case MBTRI:
266 CHKERR set_integration_pts_for_tri();
267 break;
268 case MBQUAD:
270 rule);
271 CHKERR calc_base_for_quad();
272 break;
273 default:
274 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
275 "Element type not implemented: %d", type);
276 }
277
278 } else {
279 // If rule is negative, set user defined integration points
280 CHKERR setGaussPts(order_row, order_col, order_data);
281 const size_t nb_gauss_pts = gaussPts.size2();
282 if (nb_gauss_pts) {
283 switch (type) {
284 case MBTRI:
285 CHKERR calc_base_for_tri();
286 break;
287 case MBQUAD:
288 CHKERR calc_base_for_quad();
289 break;
290 default:
291 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
292 "Element type not implemented: %d", type);
293 }
294 }
295 }
297}
298
302 // Get spaces order/base and sense of entities.
303
305
306 // H1
307 if (dataH1.spacesOnEntities[MBEDGE].test(H1)) {
308 CHKERR getEntitySense<MBEDGE>(dataH1);
309 CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
310 }
311 if (dataH1.spacesOnEntities[MBTRI].test(H1)) {
312 CHKERR getEntitySense<MBTRI>(dataH1);
313 CHKERR getEntityDataOrder<MBTRI>(dataH1, H1);
314 }
315 if (dataH1.spacesOnEntities[MBQUAD].test(H1)) {
316 CHKERR getEntitySense<MBQUAD>(dataH1);
317 CHKERR getEntityDataOrder<MBQUAD>(dataH1, H1);
318 }
319
320 // Hcurl
321 if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
322 CHKERR getEntitySense<MBEDGE>(dataHcurl);
323 CHKERR getEntityDataOrder<MBEDGE>(dataHcurl, HCURL);
324 dataHcurl.spacesOnEntities[MBEDGE].set(HCURL);
325 }
326 if (dataH1.spacesOnEntities[MBTRI].test(HCURL)) {
327 CHKERR getEntitySense<MBTRI>(dataHcurl);
328 CHKERR getEntityDataOrder<MBTRI>(dataHcurl, HCURL);
329 dataHcurl.spacesOnEntities[MBTRI].set(HCURL);
330 }
331 if (dataH1.spacesOnEntities[MBQUAD].test(HCURL)) {
332 CHKERR getEntitySense<MBQUAD>(dataHcurl);
333 CHKERR getEntityDataOrder<MBQUAD>(dataHcurl, HCURL);
334 dataHcurl.spacesOnEntities[MBQUAD].set(HCURL);
335 }
336
337 // Hdiv
338 if (dataH1.spacesOnEntities[MBTRI].test(HDIV)) {
339 CHKERR getEntitySense<MBTRI>(dataHdiv);
340 CHKERR getEntityDataOrder<MBTRI>(dataHdiv, HDIV);
341 dataHdiv.spacesOnEntities[MBTRI].set(HDIV);
342 }
343 if (dataH1.spacesOnEntities[MBQUAD].test(HDIV)) {
344 CHKERR getEntitySense<MBQUAD>(dataHdiv);
345 CHKERR getEntityDataOrder<MBQUAD>(dataHdiv, HDIV);
346 dataHdiv.spacesOnEntities[MBQUAD].set(HDIV);
347 }
348
349 // L2
350 if (dataH1.spacesOnEntities[MBTRI].test(L2)) {
351 CHKERR getEntitySense<MBTRI>(dataL2);
352 CHKERR getEntityDataOrder<MBTRI>(dataL2, L2);
353 dataL2.spacesOnEntities[MBTRI].set(L2);
354 }
355 if (dataH1.spacesOnEntities[MBQUAD].test(L2)) {
356 CHKERR getEntitySense<MBQUAD>(dataL2);
357 CHKERR getEntityDataOrder<MBQUAD>(dataL2, L2);
358 dataL2.spacesOnEntities[MBQUAD].set(L2);
359 }
360
362}
363
367
368 const size_t nb_nodes =
369 dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).size2();
370 double *shape_functions =
371 &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
372 const size_t nb_gauss_pts = gaussPts.size2();
373 coordsAtGaussPts.resize(nb_gauss_pts, 3, false);
374 for (int gg = 0; gg != nb_gauss_pts; ++gg)
375 for (int dd = 0; dd != 3; ++dd)
376 coordsAtGaussPts(gg, dd) = cblas_ddot(
377 nb_nodes, &shape_functions[nb_nodes * gg], 1, &coords[dd], 3);
378
380}
381
382MoFEMErrorCode FaceElementForcesAndSourcesCore::UserDataOperator::setPtrFE(
385 if (!(ptrFE = dynamic_cast<FaceElementForcesAndSourcesCore *>(ptr)))
386 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
387 "User operator and finite element do not work together");
389}
390
393
394 const auto type = numeredEntFiniteElementPtr->getEntType();
395 if (type != lastEvaluatedElementEntityType) {
396 switch (type) {
397 case MBTRI:
399 boost::shared_ptr<BaseFunction>(new TriPolynomialBase());
400 break;
401 case MBQUAD:
403 boost::shared_ptr<BaseFunction>(new QuadPolynomialBase());
404 break;
405 default:
407 }
410 }
411
412 // Calculate normal and tangent vectors for face geometry
415
417 if (gaussPts.size2() == 0)
419
424
425 // Iterate over operators
427
429}
430
432FaceElementForcesAndSourcesCore::UserDataOperator::loopSideVolumes(
433 const string fe_name, VolumeElementForcesAndSourcesCoreOnSide &fe_method) {
434 return loopSide(fe_name, &fe_method, 3);
435}
436
437} // namespace MoFEM
static Number< 1 > N1
static Number< 0 > N0
@ NOBASE
Definition: definitions.h:59
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#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
const int dim
#define diffN_MBQUAD2y(x)
Definition: fem_tools.h:66
#define diffN_MBQUAD1x(y)
Definition: fem_tools.h:63
#define N_MBQUAD3(x, y)
quad shape function
Definition: fem_tools.h:60
#define diffN_MBQUAD0x(y)
Definition: fem_tools.h:61
#define diffN_MBQUAD1y(x)
Definition: fem_tools.h:64
#define diffN_MBQUAD3y(x)
Definition: fem_tools.h:68
#define diffN_MBQUAD0y(x)
Definition: fem_tools.h:62
#define N_MBQUAD0(x, y)
quad shape function
Definition: fem_tools.h:57
#define diffN_MBQUAD3x(y)
Definition: fem_tools.h:67
#define diffN_MBQUAD2x(y)
Definition: fem_tools.h:65
#define N_MBQUAD2(x, y)
quad shape function
Definition: fem_tools.h:59
#define N_MBQUAD1(x, y)
quad shape function
Definition: fem_tools.h:58
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition: fem_tools.c:182
FTensor::Index< 'm', SPACE_DIM > m
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
double zeta
Definition: plastic.cpp:104
#define QUAD_2D_TABLE_SIZE
Definition: quad.h:174
static QUAD *const QUAD_2D_TABLE[]
Definition: quad.h:175
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
virtual MoFEMErrorCode calculateAreaAndNormal()
Calculate element area and normal of the face.
MoFEMErrorCode operator()()
function is run for every finite element
virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement()
Determine approximation space and order of base functions.
virtual MoFEMErrorCode calculateAreaAndNormalAtIntegrationPts()
Calculate element area and normal of the face at integration points.
virtual MoFEMErrorCode setIntegrationPts()
Set integration points.
virtual MoFEMErrorCode calculateCoordinatesAtGaussPts()
Calculate coordinate at integration points.
structure to get information form mofem into EntitiesFieldData
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
int getMaxRowOrder() const
Get max order of approximation for field in rows.
MoFEMErrorCode getSpacesAndBaseOnEntities(EntitiesFieldData &data) const
Get field approximation space and base on entities.
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
MatrixDouble coordsAtGaussPts
coordinated at gauss points
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
MoFEMErrorCode calHierarchicalBaseFunctionsOnElement()
Calculate base functions.
MatrixDouble gaussPts
Matrix of integration points.
MoFEMErrorCode calBernsteinBezierBaseFunctionsOnElement()
Calculate Bernstein-Bezier base.
int getMaxColOrder() const
Get max order of approximation for field in columns.
EntityType lastEvaluatedElementEntityType
Last evaluated type of element entity.
int getMaxDataOrder() const
Get max order of approximation for data fields.
MoFEMErrorCode createDataOnElement(EntityType type)
Create a entity data on element object.
Calculate base functions on triangle.
static MoFEMErrorCode outerProductOfEdgeIntegrationPtsForQuad(MatrixDouble &pts, const int edge0, const int edge1)
Definition: Tools.cpp:579
static constexpr std::array< double, 8 > diffShapeFunMBQUADAtCenter
Definition: Tools.hpp:212
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:104
Calculate base functions on triangle.
Base volume element used to integrate on skeleton.
int npoints
Definition: quad.h:29