v0.13.1
FaceElementForcesAndSourcesCore.cpp
Go to the documentation of this file.
1/** \file FaceElementForcesAndSourcesCore.cpp
2
3\brief Implementation of face element
4
5*/
6
7/* This file is part of MoFEM.
8 * MoFEM is free software: you can redistribute it and/or modify it under
9 * the terms of the GNU Lesser General Public License as published by the
10 * Free Software Foundation, either version 3 of the License, or (at your
11 * option) any later version.
12 *
13 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 * License for more details.
17 *
18 * You should have received a copy of the GNU Lesser General Public
19 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20
21namespace MoFEM {
22
24 Interface &m_field)
25 : ForcesAndSourcesCore(m_field),
26 meshPositionsFieldName("MESH_NODE_POSITIONS") {}
27
31
32 auto type = numeredEntFiniteElementPtr->getEntType();
33
37
38 auto get_ftensor_from_vec_3d = [](VectorDouble &v) {
40 &v[2]);
41 };
42
43 auto get_ftensor_n_diff = [&]() {
44 const auto &m = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
46 &m(0, 1));
47 };
48
49 auto get_ftensor_from_mat_3d = [](MatrixDouble &m) {
51 &m(0, 0), &m(0, 1), &m(0, 2));
52 };
53
54 if (type == MBTRI) {
55
56 const size_t nb_gauss_pts = gaussPts.size2();
57 normalsAtGaussPts.resize(nb_gauss_pts, 3);
58 tangentOneAtGaussPts.resize(nb_gauss_pts, 3);
59 tangentTwoAtGaussPts.resize(nb_gauss_pts, 3);
60
61 auto t_tan1 = get_ftensor_from_mat_3d(tangentOneAtGaussPts);
62 auto t_tan2 = get_ftensor_from_mat_3d(tangentTwoAtGaussPts);
63 auto t_normal = get_ftensor_from_mat_3d(normalsAtGaussPts);
64
65 auto t_n =
68 &tangentOne[2]);
70 &tangentTwo[2]);
71
72 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
73 t_normal(i) = t_n(i);
74 t_tan1(i) = t_t1(i);
75 t_tan2(i) = t_t2(i);
76 ++t_tan1;
77 ++t_tan2;
78 ++t_normal;
79 }
80
81 } else if (type == MBQUAD) {
82
83 EntityHandle ent = numeredEntFiniteElementPtr->getEnt();
84 CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
85 coords.resize(num_nodes * 3, false);
86 CHKERR mField.get_moab().get_coords(conn, num_nodes,
87 &*coords.data().begin());
88
89 const size_t nb_gauss_pts = gaussPts.size2();
90 normalsAtGaussPts.resize(nb_gauss_pts, 3);
91 tangentOneAtGaussPts.resize(nb_gauss_pts, 3);
92 tangentTwoAtGaussPts.resize(nb_gauss_pts, 3);
93 normalsAtGaussPts.clear();
96
97 auto t_t1 = get_ftensor_from_mat_3d(tangentOneAtGaussPts);
98 auto t_t2 = get_ftensor_from_mat_3d(tangentTwoAtGaussPts);
99 auto t_normal = get_ftensor_from_mat_3d(normalsAtGaussPts);
100
103
104 auto t_diff = get_ftensor_n_diff();
105 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
106 auto t_coords = get_ftensor_from_vec_3d(coords);
107 for (int nn = 0; nn != num_nodes; ++nn) {
108 t_t1(i) += t_coords(i) * t_diff(N0);
109 t_t2(i) += t_coords(i) * t_diff(N1);
110 ++t_diff;
111 ++t_coords;
112 }
113 t_normal(j) = FTensor::levi_civita(i, j, k) * t_t1(k) * t_t2(i);
114
115 ++t_t1;
116 ++t_t2;
117 ++t_normal;
118 }
119 } else {
120 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
121 "Element type not implemented");
122 }
123
125}
126
129
130 EntityHandle ent = numeredEntFiniteElementPtr->getEnt();
131 CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
132 coords.resize(num_nodes * 3, false);
133 CHKERR mField.get_moab().get_coords(conn, num_nodes, &*coords.data().begin());
134 nOrmal.resize(3, false);
135 tangentOne.resize(3, false);
136 tangentTwo.resize(3, false);
137
138 auto calc_normal = [&](const double *diff_ptr) {
141 &coords[0], &coords[1], &coords[2]);
143 &nOrmal[0], &nOrmal[1], &nOrmal[2]);
145 &tangentOne[0], &tangentOne[1], &tangentOne[2]);
147 &tangentTwo[0], &tangentTwo[1], &tangentTwo[2]);
149 diff_ptr, &diff_ptr[1]);
150
154
157 t_t1(i) = 0;
158 t_t2(i) = 0;
159
160 for (int nn = 0; nn != num_nodes; ++nn) {
161 t_t1(i) += t_coords(i) * t_diff(N0);
162 t_t2(i) += t_coords(i) * t_diff(N1);
163 ++t_coords;
164 ++t_diff;
165 }
166 t_normal(j) = FTensor::levi_civita(i, j, k) * t_t1(k) * t_t2(i);
167 aRea = sqrt(t_normal(i) * t_normal(i));
169 };
170
171 const double *diff_ptr;
172 switch (numeredEntFiniteElementPtr->getEntType()) {
173 case MBTRI:
174 diff_ptr = Tools::diffShapeFunMBTRI.data();
175 CHKERR calc_normal(diff_ptr);
176 // FIXME: Normal should be divided not the area for triangle!!
177 aRea /= 2;
178 break;
179 case MBQUAD:
180 diff_ptr = Tools::diffShapeFunMBQUADAtCenter.data();
181 CHKERR calc_normal(diff_ptr);
182 break;
183 default:
184 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
185 "Element type not implemented");
186 }
187
189}
190
193 // Set integration points
194 int order_data = getMaxDataOrder();
195 int order_row = getMaxRowOrder();
196 int order_col = getMaxColOrder();
197 int rule = getRule(order_row, order_col, order_data);
198
199 auto set_integration_pts_for_tri = [&]() {
201 if (rule < QUAD_2D_TABLE_SIZE) {
202 if (QUAD_2D_TABLE[rule]->dim != 2) {
203 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
204 }
205 if (QUAD_2D_TABLE[rule]->order < rule) {
206 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
207 "wrong order %d != %d", QUAD_2D_TABLE[rule]->order, rule);
208 }
209 const size_t nb_gauss_pts = QUAD_2D_TABLE[rule]->npoints;
210 gaussPts.resize(3, nb_gauss_pts, false);
211 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[1], 3,
212 &gaussPts(0, 0), 1);
213 cblas_dcopy(nb_gauss_pts, &QUAD_2D_TABLE[rule]->points[2], 3,
214 &gaussPts(1, 0), 1);
215 cblas_dcopy(nb_gauss_pts, QUAD_2D_TABLE[rule]->weights, 1,
216 &gaussPts(2, 0), 1);
217 dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 3,
218 false);
219 double *shape_ptr =
220 &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
221 cblas_dcopy(3 * nb_gauss_pts, QUAD_2D_TABLE[rule]->points, 1, shape_ptr,
222 1);
223 dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(3, 2, false);
224 std::copy(
226 dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).data().begin());
227 } else {
228 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
229 "rule > quadrature order %d < %d", rule, QUAD_2D_TABLE_SIZE);
230 }
232 };
233
234 auto calc_base_for_quad = [&]() {
236 const size_t nb_gauss_pts = gaussPts.size2();
237 auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
238 auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
239 base.resize(nb_gauss_pts, 4, false);
240 diff_base.resize(nb_gauss_pts, 8, false);
241 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
242 const double ksi = gaussPts(0, gg);
243 const double zeta = gaussPts(1, gg);
244 base(gg, 0) = N_MBQUAD0(ksi, zeta);
245 base(gg, 1) = N_MBQUAD1(ksi, zeta);
246 base(gg, 2) = N_MBQUAD2(ksi, zeta);
247 base(gg, 3) = N_MBQUAD3(ksi, zeta);
248 diff_base(gg, 0) = diffN_MBQUAD0x(zeta);
249 diff_base(gg, 1) = diffN_MBQUAD0y(ksi);
250 diff_base(gg, 2) = diffN_MBQUAD1x(zeta);
251 diff_base(gg, 3) = diffN_MBQUAD1y(ksi);
252 diff_base(gg, 4) = diffN_MBQUAD2x(zeta);
253 diff_base(gg, 5) = diffN_MBQUAD2y(ksi);
254 diff_base(gg, 6) = diffN_MBQUAD3x(zeta);
255 diff_base(gg, 7) = diffN_MBQUAD3y(ksi);
256 }
258 };
259
260 const auto type = numeredEntFiniteElementPtr->getEntType();
261
262 if (rule >= 0) {
263 switch (type) {
264 case MBTRI:
265 CHKERR set_integration_pts_for_tri();
266 break;
267 case MBQUAD:
269 rule);
270 CHKERR calc_base_for_quad();
271 break;
272 default:
273 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
274 "Element type not implemented: %d", type);
275 }
276
277 } else {
278 // If rule is negative, set user defined integration points
279 CHKERR setGaussPts(order_row, order_col, order_data);
280 const size_t nb_gauss_pts = gaussPts.size2();
281 auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
282 auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
283 if (nb_gauss_pts) {
284 switch (type) {
285 case MBTRI:
286 base.resize(nb_gauss_pts, 3, false);
287 diff_base.resize(3, 2, false);
288 CHKERR ShapeMBTRI(&*base.data().begin(), &gaussPts(0, 0),
289 &gaussPts(1, 0), nb_gauss_pts);
290 std::copy(
292 dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).data().begin());
293 break;
294 case MBQUAD:
295 CHKERR calc_base_for_quad();
296 break;
297 default:
298 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
299 "Element type not implemented: %d", type);
300 }
301 }
302 }
304}
305
309 // Get spaces order/base and sense of entities.
310
312
313 // H1
314 if (dataH1.spacesOnEntities[MBEDGE].test(H1)) {
315 CHKERR getEntitySense<MBEDGE>(dataH1);
316 CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
317 }
318 if (dataH1.spacesOnEntities[MBTRI].test(H1)) {
319 CHKERR getEntitySense<MBTRI>(dataH1);
320 CHKERR getEntityDataOrder<MBTRI>(dataH1, H1);
321 }
322 if (dataH1.spacesOnEntities[MBQUAD].test(H1)) {
323 CHKERR getEntitySense<MBQUAD>(dataH1);
324 CHKERR getEntityDataOrder<MBQUAD>(dataH1, H1);
325 }
326
327 // Hcurl
328 if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
329 CHKERR getEntitySense<MBEDGE>(dataHcurl);
330 CHKERR getEntityDataOrder<MBEDGE>(dataHcurl, HCURL);
331 dataHcurl.spacesOnEntities[MBEDGE].set(HCURL);
332 }
333 if (dataH1.spacesOnEntities[MBTRI].test(HCURL)) {
334 CHKERR getEntitySense<MBTRI>(dataHcurl);
335 CHKERR getEntityDataOrder<MBTRI>(dataHcurl, HCURL);
336 dataHcurl.spacesOnEntities[MBTRI].set(HCURL);
337 }
338 if (dataH1.spacesOnEntities[MBQUAD].test(HCURL)) {
339 CHKERR getEntitySense<MBQUAD>(dataHcurl);
340 CHKERR getEntityDataOrder<MBQUAD>(dataHcurl, HCURL);
341 dataHcurl.spacesOnEntities[MBQUAD].set(HCURL);
342 }
343
344 // Hdiv
345 if (dataH1.spacesOnEntities[MBTRI].test(HDIV)) {
346 CHKERR getEntitySense<MBTRI>(dataHdiv);
347 CHKERR getEntityDataOrder<MBTRI>(dataHdiv, HDIV);
348 dataHdiv.spacesOnEntities[MBTRI].set(HDIV);
349 }
350 if (dataH1.spacesOnEntities[MBQUAD].test(HDIV)) {
351 CHKERR getEntitySense<MBQUAD>(dataHdiv);
352 CHKERR getEntityDataOrder<MBQUAD>(dataHdiv, HDIV);
353 dataHdiv.spacesOnEntities[MBQUAD].set(HDIV);
354 }
355
356 // L2
357 if (dataH1.spacesOnEntities[MBTRI].test(L2)) {
358 CHKERR getEntitySense<MBTRI>(dataL2);
359 CHKERR getEntityDataOrder<MBTRI>(dataL2, L2);
360 dataL2.spacesOnEntities[MBTRI].set(L2);
361 }
362 if (dataH1.spacesOnEntities[MBQUAD].test(L2)) {
363 CHKERR getEntitySense<MBQUAD>(dataL2);
364 CHKERR getEntityDataOrder<MBQUAD>(dataL2, L2);
365 dataL2.spacesOnEntities[MBQUAD].set(L2);
366 }
367
369}
370
374
375 const size_t nb_nodes =
376 dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).size2();
377 double *shape_functions =
378 &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
379 const size_t nb_gauss_pts = gaussPts.size2();
380 coordsAtGaussPts.resize(nb_gauss_pts, 3, false);
381 for (int gg = 0; gg != nb_gauss_pts; ++gg)
382 for (int dd = 0; dd != 3; ++dd)
383 coordsAtGaussPts(gg, dd) = cblas_ddot(
384 nb_nodes, &shape_functions[nb_nodes * gg], 1, &coords[dd], 3);
385
387}
388
392 if (!(ptrFE = dynamic_cast<FaceElementForcesAndSourcesCore *>(ptr)))
393 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
394 "User operator and finite element do not work together");
396}
397
400
401 const auto type = numeredEntFiniteElementPtr->getEntType();
403 switch (type) {
404 case MBTRI:
406 boost::shared_ptr<BaseFunction>(new TriPolynomialBase());
407 break;
408 case MBQUAD:
410 boost::shared_ptr<BaseFunction>(new QuadPolynomialBase());
411 break;
412 default:
414 }
417 }
418
419 // Calculate normal and tangent vectors for face geometry
422
424 if (gaussPts.size2() == 0)
426
431
432 // Iterate over operators
434
436}
437
440 const string fe_name, VolumeElementForcesAndSourcesCoreOnSide &fe_method) {
441 return loopSide(fe_name, &fe_method, 3);
442}
443
444} // namespace MoFEM
static Number< 1 > N1
static Number< 0 > N0
@ NOBASE
Definition: definitions.h:72
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ H1
continuous field
Definition: definitions.h:98
@ HCURL
field with continuous tangents
Definition: definitions.h:99
@ HDIV
field with continuous normal traction
Definition: definitions.h:100
#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
Definition: definitions.h:44
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:45
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
const int dim
#define diffN_MBQUAD2y(x)
Definition: fem_tools.h:76
#define diffN_MBQUAD1x(y)
Definition: fem_tools.h:73
#define N_MBQUAD3(x, y)
quad shape function
Definition: fem_tools.h:70
#define diffN_MBQUAD0x(y)
Definition: fem_tools.h:71
#define diffN_MBQUAD1y(x)
Definition: fem_tools.h:74
#define diffN_MBQUAD3y(x)
Definition: fem_tools.h:78
#define diffN_MBQUAD0y(x)
Definition: fem_tools.h:72
#define N_MBQUAD0(x, y)
quad shape function
Definition: fem_tools.h:67
#define diffN_MBQUAD3x(y)
Definition: fem_tools.h:77
#define diffN_MBQUAD2x(y)
Definition: fem_tools.h:75
#define N_MBQUAD2(x, y)
quad shape function
Definition: fem_tools.h:69
#define N_MBQUAD1(x, y)
quad shape function
Definition: fem_tools.h:68
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition: fem_tools.c:194
FTensor::Index< 'i', SPACE_DIM > i
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
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
#define QUAD_2D_TABLE_SIZE
Definition: quad.h:174
static QUAD *const QUAD_2D_TABLE[]
Definition: quad.h:175
FTensor::Index< 'm', 3 > m
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
MoFEMErrorCode loopSideVolumes(const string fe_name, VolumeElementForcesAndSourcesCoreOnSide &fe_method)
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:591
static constexpr std::array< double, 8 > diffShapeFunMBQUADAtCenter
Definition: Tools.hpp:221
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:113
Calculate base functions on triangle.
Base volume element used to integrate on skeleton.
int npoints
Definition: quad.h:29