v0.13.2
Loading...
Searching...
No Matches
TetPolynomialBase.cpp
Go to the documentation of this file.
1/** \file TetPolynomialBase.cpp
2\brief Implementation of hierarchical bases on tetrahedral
3
4A l2, h1, h-div and h-curl spaces are implemented.
5
6*/
7
8using namespace MoFEM;
9
10MoFEMErrorCode
11TetPolynomialBase::query_interface(boost::typeindex::type_index type_index,
12 UnknownInterface **iface) const {
13
14 MoFEMFunctionBeginHot;
15 *iface = const_cast<TetPolynomialBase *>(this);
16 MoFEMFunctionReturnHot(0);
17}
18
19MoFEMErrorCode TetPolynomialBase::getValueH1(MatrixDouble &pts) {
20 MoFEMFunctionBegin;
21
22 switch (cTx->bAse) {
23 case AINSWORTH_LEGENDRE_BASE:
24 case AINSWORTH_LOBATTO_BASE:
25 CHKERR getValueH1AinsworthBase(pts);
26 break;
27 case AINSWORTH_BERNSTEIN_BEZIER_BASE:
28 CHKERR getValueH1BernsteinBezierBase(pts);
29 break;
30 default:
31 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
32 }
33
34 MoFEMFunctionReturn(0);
35}
36
37MoFEMErrorCode TetPolynomialBase::getValueH1AinsworthBase(MatrixDouble &pts) {
38 MoFEMFunctionBegin;
39
40 EntitiesFieldData &data = cTx->dAta;
41 const FieldApproximationBase base = cTx->bAse;
42 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
43 double *diffL, const int dim) =
44 cTx->basePolynomialsType0;
45
46 int nb_gauss_pts = pts.size2();
47
48 int sense[6], order[6];
49 if (data.spacesOnEntities[MBEDGE].test(H1)) {
50 // edges
51 if (data.dataOnEntities[MBEDGE].size() != 6) {
52 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
53 }
54 double *h1_edge_n[6], *diff_h1_egde_n[6];
55 for (int ee = 0; ee != 6; ++ee) {
56 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
57 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
58 "data inconsistency");
59 }
60 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
61 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
62 int nb_dofs = NBEDGE_H1(data.dataOnEntities[MBEDGE][ee].getOrder());
63 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
64 false);
65 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
66 3 * nb_dofs, false);
67 h1_edge_n[ee] =
68 &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
69 diff_h1_egde_n[ee] =
70 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
71 }
72 CHKERR H1_EdgeShapeFunctions_MBTET(
73 sense, order,
74 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
75 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
76 h1_edge_n, diff_h1_egde_n, nb_gauss_pts, base_polynomials);
77 } else {
78 for (int ee = 0; ee != 6; ++ee) {
79 data.dataOnEntities[MBEDGE][ee].getN(base).resize(0, 0, false);
80 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(0, 0, false);
81 }
82 }
83
84 if (data.spacesOnEntities[MBTRI].test(H1)) {
85 // faces
86 if (data.dataOnEntities[MBTRI].size() != 4) {
87 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
88 }
89 double *h1_face_n[4], *diff_h1_face_n[4];
90 for (int ff = 0; ff != 4; ++ff) {
91 if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
92 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
93 "data inconsistency");
94 }
95 int nb_dofs = NBFACETRI_H1(data.dataOnEntities[MBTRI][ff].getOrder());
96 order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
97 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, nb_dofs,
98 false);
99 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
100 3 * nb_dofs, false);
101 h1_face_n[ff] =
102 &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
103 diff_h1_face_n[ff] =
104 &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
105 }
106 if (data.facesNodes.size1() != 4) {
107 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
108 }
109 if (data.facesNodes.size2() != 3) {
110 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
111 }
112 CHKERR H1_FaceShapeFunctions_MBTET(
113 &*data.facesNodes.data().begin(), order,
114 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
115 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
116 h1_face_n, diff_h1_face_n, nb_gauss_pts, base_polynomials);
117
118 } else {
119 for (int ff = 0; ff != 4; ++ff) {
120 data.dataOnEntities[MBTRI][ff].getN(base).resize(0, false);
121 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(0, 0, false);
122 }
123 }
124
125 if (data.spacesOnEntities[MBTET].test(H1)) {
126 // volume
127 int order = data.dataOnEntities[MBTET][0].getOrder();
128 int nb_vol_dofs = NBVOLUMETET_H1(order);
129 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, nb_vol_dofs,
130 false);
131 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
132 3 * nb_vol_dofs, false);
133 CHKERR H1_VolumeShapeFunctions_MBTET(
134 data.dataOnEntities[MBTET][0].getOrder(),
135 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
136 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
137 &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
138 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
139 nb_gauss_pts, base_polynomials);
140 } else {
141 data.dataOnEntities[MBTET][0].getN(base).resize(0, 0, false);
142 data.dataOnEntities[MBTET][0].getDiffN(base).resize(0, 0, false);
143 }
144
145 MoFEMFunctionReturn(0);
146}
147
148MoFEMErrorCode
149TetPolynomialBase::getValueH1BernsteinBezierBase(MatrixDouble &pts) {
150 MoFEMFunctionBegin;
151
152 EntitiesFieldData &data = cTx->dAta;
153 const std::string field_name = cTx->fieldName;
154 const int nb_gauss_pts = pts.size2();
155
156 if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
157 (unsigned int)nb_gauss_pts)
158 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
159 "Base functions or nodes has wrong number of integration points "
160 "for base %s",
161 ApproximationBaseNames[NOBASE]);
162 auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
163
164 auto get_alpha = [field_name](auto &data) -> MatrixInt & {
165 auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
166 if (!ptr)
167 ptr.reset(new MatrixInt());
168 return *ptr;
169 };
170
171 auto get_base = [field_name](auto &data) -> MatrixDouble & {
172 auto &ptr = data.getBBNSharedPtr(field_name);
173 if (!ptr)
174 ptr.reset(new MatrixDouble());
175 return *ptr;
176 };
177
178 auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
179 auto &ptr = data.getBBDiffNSharedPtr(field_name);
180 if (!ptr)
181 ptr.reset(new MatrixDouble());
182 return *ptr;
183 };
184
185 auto get_alpha_by_name_ptr =
186 [](auto &data,
187 const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
188 return data.getBBAlphaIndicesSharedPtr(field_name);
189 };
190
191 auto get_base_by_name_ptr =
192 [](auto &data,
193 const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
194 return data.getBBNSharedPtr(field_name);
195 };
196
197 auto get_diff_base_by_name_ptr =
198 [](auto &data,
199 const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
200 return data.getBBDiffNSharedPtr(field_name);
201 };
202
203 auto get_alpha_by_order_ptr =
204 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
205 return data.getBBAlphaIndicesByOrderSharedPtr(o);
206 };
207
208 auto get_base_by_order_ptr =
209 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
210 return data.getBBNByOrderSharedPtr(o);
211 };
212
213 auto get_diff_base_by_order_ptr =
214 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
215 return data.getBBDiffNByOrderSharedPtr(o);
216 };
217
218 auto &vert_ent_data = data.dataOnEntities[MBVERTEX][0];
219 auto &vertex_alpha = get_alpha(vert_ent_data);
220 vertex_alpha.resize(4, 4, false);
221 vertex_alpha.clear();
222 for (int n = 0; n != 4; ++n)
223 vertex_alpha(n, n) = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n];
224
225 auto &vert_get_n = get_base(vert_ent_data);
226 auto &vert_get_diff_n = get_diff_base(vert_ent_data);
227 vert_get_n.resize(nb_gauss_pts, 4, false);
228 vert_get_diff_n.resize(nb_gauss_pts, 12, false);
229 CHKERR BernsteinBezier::baseFunctionsTet(
230 1, lambda.size1(), vertex_alpha.size1(), &vertex_alpha(0, 0),
231 &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &vert_get_n(0, 0),
232 &vert_get_diff_n(0, 0));
233 for (int n = 0; n != 4; ++n) {
234 const double f = boost::math::factorial<double>(
235 data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n]);
236 for (int g = 0; g != nb_gauss_pts; ++g) {
237 vert_get_n(g, n) *= f;
238 for (int d = 0; d != 3; ++d)
239 vert_get_diff_n(g, 3 * n + d) *= f;
240 }
241 }
242
243 // edges
244 if (data.spacesOnEntities[MBEDGE].test(H1)) {
245 if (data.dataOnEntities[MBEDGE].size() != 6)
246 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
247 "Wrong size of ent data");
248
249 constexpr int edges_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
250 {0, 3}, {1, 3}, {2, 3}};
251 for (int ee = 0; ee != 6; ++ee) {
252 auto &ent_data = data.dataOnEntities[MBEDGE][ee];
253 const int sense = ent_data.getSense();
254 if (sense == 0)
255 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
256 "Sense of the edge unknown");
257 const int order = ent_data.getOrder();
258 const int nb_dofs = NBEDGE_H1(order);
259
260 if (nb_dofs) {
261 if (get_alpha_by_order_ptr(ent_data, order)) {
262 get_alpha_by_name_ptr(ent_data, field_name) =
263 get_alpha_by_order_ptr(ent_data, order);
264 get_base_by_name_ptr(ent_data, field_name) =
265 get_base_by_order_ptr(ent_data, order);
266 get_diff_base_by_name_ptr(ent_data, field_name) =
267 get_diff_base_by_order_ptr(ent_data, order);
268 } else {
269 auto &get_n = get_base(ent_data);
270 auto &get_diff_n = get_diff_base(ent_data);
271 get_n.resize(nb_gauss_pts, nb_dofs, false);
272 get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
273
274 auto &edge_alpha = get_alpha(data.dataOnEntities[MBEDGE][ee]);
275 edge_alpha.resize(nb_dofs, 4, false);
276 CHKERR BernsteinBezier::generateIndicesEdgeTet(ee, order,
277 &edge_alpha(0, 0));
278 if (sense == -1) {
279 for (int i = 0; i != edge_alpha.size1(); ++i) {
280 int a = edge_alpha(i, edges_nodes[ee][0]);
281 edge_alpha(i, edges_nodes[ee][0]) =
282 edge_alpha(i, edges_nodes[ee][1]);
283 edge_alpha(i, edges_nodes[ee][1]) = a;
284 }
285 }
286 CHKERR BernsteinBezier::baseFunctionsTet(
287 order, lambda.size1(), edge_alpha.size1(), &edge_alpha(0, 0),
288 &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
289 &get_diff_n(0, 0));
290
291 get_alpha_by_order_ptr(ent_data, order) =
292 get_alpha_by_name_ptr(ent_data, field_name);
293 get_base_by_order_ptr(ent_data, order) =
294 get_base_by_name_ptr(ent_data, field_name);
295 get_diff_base_by_order_ptr(ent_data, order) =
296 get_diff_base_by_name_ptr(ent_data, field_name);
297 }
298 }
299 }
300 } else {
301 for (int ee = 0; ee != 6; ++ee) {
302 auto &ent_data = data.dataOnEntities[MBEDGE][ee];
303 ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
304 auto &get_n = get_base(ent_data);
305 auto &get_diff_n = get_diff_base(ent_data);
306 get_n.resize(nb_gauss_pts, 0, false);
307 get_diff_n.resize(nb_gauss_pts, 0, false);
308 }
309 }
310
311 // face
312 if (data.spacesOnEntities[MBTRI].test(H1)) {
313 if (data.dataOnEntities[MBTRI].size() != 4)
314 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
315 "Wrong size of ent data");
316 if (data.facesNodes.size1() != 4)
317 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
318 if (data.facesNodes.size2() != 3)
319 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
320
321 for (int ff = 0; ff != 4; ++ff) {
322 auto &ent_data = data.dataOnEntities[MBTRI][ff];
323 const int order = ent_data.getOrder();
324 const int nb_dofs = NBFACETRI_H1(order);
325
326 if (nb_dofs) {
327 if (get_alpha_by_order_ptr(ent_data, order)) {
328 get_alpha_by_name_ptr(ent_data, field_name) =
329 get_alpha_by_order_ptr(ent_data, order);
330 get_base_by_name_ptr(ent_data, field_name) =
331 get_base_by_order_ptr(ent_data, order);
332 get_diff_base_by_name_ptr(ent_data, field_name) =
333 get_diff_base_by_order_ptr(ent_data, order);
334 } else {
335
336 auto &get_n = get_base(ent_data);
337 auto &get_diff_n = get_diff_base(ent_data);
338 get_n.resize(nb_gauss_pts, nb_dofs, false);
339 get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
340
341 auto &face_alpha = get_alpha(ent_data);
342 face_alpha.resize(nb_dofs, 4, false);
343
344 CHKERR BernsteinBezier::generateIndicesTriTet(ff, order,
345 &face_alpha(0, 0));
346 senseFaceAlpha.resize(face_alpha.size1(), face_alpha.size2(), false);
347 senseFaceAlpha.clear();
348 constexpr int tri_nodes[4][3] = {
349 {0, 1, 3}, {1, 2, 3}, {0, 2, 3}, {0, 1, 2}};
350 for (int d = 0; d != nb_dofs; ++d)
351 for (int n = 0; n != 3; ++n)
352 senseFaceAlpha(d, data.facesNodes(ff, n)) =
353 face_alpha(d, tri_nodes[ff][n]);
354 face_alpha.swap(senseFaceAlpha);
355 CHKERR BernsteinBezier::baseFunctionsTet(
356 order, lambda.size1(), face_alpha.size1(), &face_alpha(0, 0),
357 &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
358 &get_diff_n(0, 0));
359
360 get_alpha_by_order_ptr(ent_data, order) =
361 get_alpha_by_name_ptr(ent_data, field_name);
362 get_base_by_order_ptr(ent_data, order) =
363 get_base_by_name_ptr(ent_data, field_name);
364 get_diff_base_by_order_ptr(ent_data, order) =
365 get_diff_base_by_name_ptr(ent_data, field_name);
366 }
367 }
368 }
369 } else {
370 for (int ff = 0; ff != 4; ++ff) {
371 auto &ent_data = data.dataOnEntities[MBTRI][ff];
372 ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
373 auto &get_n = get_base(ent_data);
374 auto &get_diff_n = get_diff_base(ent_data);
375 get_n.resize(nb_gauss_pts, 0, false);
376 get_diff_n.resize(nb_gauss_pts, 0, false);
377 }
378 }
379
380 if (data.spacesOnEntities[MBTET].test(H1)) {
381 if (data.dataOnEntities[MBTET].size() != 1)
382 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
383 "Wrong size ent of ent data");
384
385 auto &ent_data = data.dataOnEntities[MBTET][0];
386 const int order = ent_data.getOrder();
387 const int nb_dofs = NBVOLUMETET_H1(order);
388 if (get_alpha_by_order_ptr(ent_data, order)) {
389 get_alpha_by_name_ptr(ent_data, field_name) =
390 get_alpha_by_order_ptr(ent_data, order);
391 get_base_by_name_ptr(ent_data, field_name) =
392 get_base_by_order_ptr(ent_data, order);
393 get_diff_base_by_name_ptr(ent_data, field_name) =
394 get_diff_base_by_order_ptr(ent_data, order);
395 } else {
396
397 auto &get_n = get_base(ent_data);
398 auto &get_diff_n = get_diff_base(ent_data);
399 get_n.resize(nb_gauss_pts, nb_dofs, false);
400 get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
401 if (nb_dofs) {
402 auto &tet_alpha = get_alpha(ent_data);
403 tet_alpha.resize(nb_dofs, 4, false);
404
405 CHKERR BernsteinBezier::generateIndicesTetTet(order, &tet_alpha(0, 0));
406 CHKERR BernsteinBezier::baseFunctionsTet(
407 order, lambda.size1(), tet_alpha.size1(), &tet_alpha(0, 0),
408 &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
409 &get_diff_n(0, 0));
410
411 get_alpha_by_order_ptr(ent_data, order) =
412 get_alpha_by_name_ptr(ent_data, field_name);
413 get_base_by_order_ptr(ent_data, order) =
414 get_base_by_name_ptr(ent_data, field_name);
415 get_diff_base_by_order_ptr(ent_data, order) =
416 get_diff_base_by_name_ptr(ent_data, field_name);
417 }
418 }
419 } else {
420 auto &ent_data = data.dataOnEntities[MBTET][0];
421 ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
422 auto &get_n = get_base(ent_data);
423 auto &get_diff_n = get_diff_base(ent_data);
424 get_n.resize(nb_gauss_pts, 0, false);
425 get_diff_n.resize(nb_gauss_pts, 0, false);
426 }
427
428 MoFEMFunctionReturn(0);
429}
430
431MoFEMErrorCode TetPolynomialBase::getValueL2(MatrixDouble &pts) {
432 MoFEMFunctionBegin;
433
434 switch (cTx->bAse) {
435 case AINSWORTH_LEGENDRE_BASE:
436 case AINSWORTH_LOBATTO_BASE:
437 CHKERR getValueL2AinsworthBase(pts);
438 break;
439 case AINSWORTH_BERNSTEIN_BEZIER_BASE:
440 CHKERR getValueL2BernsteinBezierBase(pts);
441 break;
442 default:
443 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
444 }
445
446 MoFEMFunctionReturn(0);
447}
448
449MoFEMErrorCode TetPolynomialBase::getValueL2AinsworthBase(MatrixDouble &pts) {
450 MoFEMFunctionBegin;
451
452 EntitiesFieldData &data = cTx->dAta;
453 const FieldApproximationBase base = cTx->bAse;
454 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
455 double *diffL, const int dim) =
456 cTx->basePolynomialsType0;
457
458 int nb_gauss_pts = pts.size2();
459
460 data.dataOnEntities[MBTET][0].getN(base).resize(
461 nb_gauss_pts,
462 NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getOrder()), false);
463 data.dataOnEntities[MBTET][0].getDiffN(base).resize(
464 nb_gauss_pts,
465 3 * NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getOrder()), false);
466
467 CHKERR L2_Ainsworth_ShapeFunctions_MBTET(
468 data.dataOnEntities[MBTET][0].getOrder(),
469 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
470 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
471 &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
472 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
473 nb_gauss_pts, base_polynomials);
474
475 MoFEMFunctionReturn(0);
476}
477
478MoFEMErrorCode
479TetPolynomialBase::getValueL2BernsteinBezierBase(MatrixDouble &pts) {
480 MoFEMFunctionBegin;
481
482 EntitiesFieldData &data = cTx->dAta;
483 const std::string field_name = cTx->fieldName;
484 const int nb_gauss_pts = pts.size2();
485
486 if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
487 (unsigned int)nb_gauss_pts)
488 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
489 "Base functions or nodes has wrong number of integration points "
490 "for base %s",
491 ApproximationBaseNames[NOBASE]);
492 auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
493
494 auto get_alpha = [field_name](auto &data) -> MatrixInt & {
495 auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
496 if (!ptr)
497 ptr.reset(new MatrixInt());
498 return *ptr;
499 };
500
501 auto get_base = [field_name](auto &data) -> MatrixDouble & {
502 auto &ptr = data.getBBNSharedPtr(field_name);
503 if (!ptr)
504 ptr.reset(new MatrixDouble());
505 return *ptr;
506 };
507
508 auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
509 auto &ptr = data.getBBDiffNSharedPtr(field_name);
510 if (!ptr)
511 ptr.reset(new MatrixDouble());
512 return *ptr;
513 };
514
515 auto get_alpha_by_name_ptr =
516 [](auto &data,
517 const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
518 return data.getBBAlphaIndicesSharedPtr(field_name);
519 };
520
521 auto get_base_by_name_ptr =
522 [](auto &data,
523 const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
524 return data.getBBNSharedPtr(field_name);
525 };
526
527 auto get_diff_base_by_name_ptr =
528 [](auto &data,
529 const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
530 return data.getBBDiffNSharedPtr(field_name);
531 };
532
533 auto get_alpha_by_order_ptr =
534 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
535 return data.getBBAlphaIndicesByOrderSharedPtr(o);
536 };
537
538 auto get_base_by_order_ptr =
539 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
540 return data.getBBNByOrderSharedPtr(o);
541 };
542
543 auto get_diff_base_by_order_ptr =
544 [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
545 return data.getBBDiffNByOrderSharedPtr(o);
546 };
547
548 if (data.spacesOnEntities[MBTET].test(L2)) {
549 if (data.dataOnEntities[MBTET].size() != 1)
550 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
551 "Wrong size ent of ent data");
552
553 auto &ent_data = data.dataOnEntities[MBTET][0];
554 const int order = ent_data.getOrder();
555 const int nb_dofs = NBVOLUMETET_L2(order);
556
557 if (get_alpha_by_order_ptr(ent_data, order)) {
558 get_alpha_by_name_ptr(ent_data, field_name) =
559 get_alpha_by_order_ptr(ent_data, order);
560 get_base_by_name_ptr(ent_data, field_name) =
561 get_base_by_order_ptr(ent_data, order);
562 get_diff_base_by_name_ptr(ent_data, field_name) =
563 get_diff_base_by_order_ptr(ent_data, order);
564 } else {
565
566 auto &get_n = get_base(ent_data);
567 auto &get_diff_n = get_diff_base(ent_data);
568 get_n.resize(nb_gauss_pts, nb_dofs, false);
569 get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
570
571 if (nb_dofs) {
572
573 if (order == 0) {
574
575 if (nb_dofs != 1)
576 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
577 "Inconsistent number of DOFs");
578
579 auto &tri_alpha = get_alpha(ent_data);
580 tri_alpha.clear();
581 get_n(0, 0) = 1;
582 get_diff_n.clear();
583
584 } else {
585
586 if (nb_dofs != 4 + 6 * NBEDGE_H1(order) + 4 * NBFACETRI_H1(order) +
587 NBVOLUMETET_H1(order) ||
588 nb_dofs != NBVOLUMETET_L2(order))
589 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
590 "Inconsistent number of DOFs");
591
592 auto &tet_alpha = get_alpha(ent_data);
593 tet_alpha.resize(nb_dofs, 4, false);
594
595 CHKERR BernsteinBezier::generateIndicesVertexTet(order,
596 &tet_alpha(0, 0));
597 if (order > 1) {
598 std::array<int, 6> edge_n{order, order, order, order, order, order};
599 std::array<int *, 6> tet_edge_ptr{
600 &tet_alpha(4, 0),
601 &tet_alpha(4 + 1 * NBEDGE_H1(order), 0),
602 &tet_alpha(4 + 2 * NBEDGE_H1(order), 0),
603 &tet_alpha(4 + 3 * NBEDGE_H1(order), 0),
604 &tet_alpha(4 + 4 * NBEDGE_H1(order), 0),
605 &tet_alpha(4 + 5 * NBEDGE_H1(order), 0)};
606 CHKERR BernsteinBezier::generateIndicesEdgeTet(edge_n.data(),
607 tet_edge_ptr.data());
608 if (order > 2) {
609 std::array<int, 6> face_n{order, order, order, order};
610 std::array<int *, 6> tet_face_ptr{
611 &tet_alpha(4 + 6 * NBEDGE_H1(order), 0),
612 &tet_alpha(4 + 6 * NBEDGE_H1(order) + 1 * NBFACETRI_H1(order),
613 0),
614 &tet_alpha(4 + 6 * NBEDGE_H1(order) + 2 * NBFACETRI_H1(order),
615 0),
616 &tet_alpha(4 + 6 * NBEDGE_H1(order) + 3 * NBFACETRI_H1(order),
617 0),
618 };
619 CHKERR BernsteinBezier::generateIndicesTriTet(
620 face_n.data(), tet_face_ptr.data());
621 if (order > 3)
622 CHKERR BernsteinBezier::generateIndicesTetTet(
623 order,
624 &tet_alpha(
625 4 + 6 * NBEDGE_H1(order) + 4 * NBFACETRI_H1(order), 0));
626 }
627 }
628
629 CHKERR BernsteinBezier::baseFunctionsTet(
630 order, lambda.size1(), tet_alpha.size1(), &tet_alpha(0, 0),
631 &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
632 &get_diff_n(0, 0));
633
634 get_alpha_by_order_ptr(ent_data, order) =
635 get_alpha_by_name_ptr(ent_data, field_name);
636 get_base_by_order_ptr(ent_data, order) =
637 get_base_by_name_ptr(ent_data, field_name);
638 get_diff_base_by_order_ptr(ent_data, order) =
639 get_diff_base_by_name_ptr(ent_data, field_name);
640 }
641 }
642 }
643 } else {
644 auto &ent_data = data.dataOnEntities[MBTET][0];
645 ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
646 auto &get_n = get_base(ent_data);
647 auto &get_diff_n = get_diff_base(ent_data);
648 get_n.resize(nb_gauss_pts, 0, false);
649 get_diff_n.resize(nb_gauss_pts, 0, false);
650 }
651
652 MoFEMFunctionReturn(0);
653}
654
655MoFEMErrorCode TetPolynomialBase::getValueHdivAinsworthBase(MatrixDouble &pts) {
656 MoFEMFunctionBegin;
657
658 EntitiesFieldData &data = cTx->dAta;
659 const FieldApproximationBase base = cTx->bAse;
660 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
661 double *diffL, const int dim) =
662 cTx->basePolynomialsType0;
663
664 int nb_gauss_pts = pts.size2();
665
666 // face shape functions
667
668 double *phi_f_e[4][3];
669 double *phi_f[4];
670 double *diff_phi_f_e[4][3];
671 double *diff_phi_f[4];
672
673 N_face_edge.resize(4, 3, false);
674 N_face_bubble.resize(4, false);
675 diffN_face_edge.resize(4, 3, false);
676 diffN_face_bubble.resize(4, false);
677
678 int faces_order[4];
679 for (int ff = 0; ff != 4; ++ff) {
680 if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
681 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
682 }
683 faces_order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
684 // three edges on face
685 for (int ee = 0; ee < 3; ee++) {
686 N_face_edge(ff, ee).resize(
687 nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(faces_order[ff]),
688 false);
689 diffN_face_edge(ff, ee).resize(
690 nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_EDGE_HDIV(faces_order[ff]),
691 false);
692 phi_f_e[ff][ee] = &*N_face_edge(ff, ee).data().begin();
693 diff_phi_f_e[ff][ee] = &*diffN_face_edge(ff, ee).data().begin();
694 }
695 N_face_bubble[ff].resize(nb_gauss_pts,
696 3 * NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]),
697 false);
698 diffN_face_bubble[ff].resize(
699 nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]),
700 false);
701 phi_f[ff] = &*(N_face_bubble[ff].data().begin());
702 diff_phi_f[ff] = &*(diffN_face_bubble[ff].data().begin());
703 }
704
705 CHKERR Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET(
706 &data.facesNodes(0, 0), faces_order,
707 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
708 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(), phi_f_e,
709 diff_phi_f_e, nb_gauss_pts, base_polynomials);
710
711 CHKERR Hdiv_Ainsworth_FaceBubbleShapeFunctions(
712 &data.facesNodes(0, 0), faces_order,
713 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
714 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(), phi_f,
715 diff_phi_f, nb_gauss_pts, base_polynomials);
716
717 // volume shape functions
718
719 double *phi_v_e[6];
720 double *phi_v_f[4];
721 double *phi_v;
722 double *diff_phi_v_e[6];
723 double *diff_phi_v_f[4];
724 double *diff_phi_v;
725
726 int volume_order = data.dataOnEntities[MBTET][0].getOrder();
727
728 N_volume_edge.resize(6, false);
729 diffN_volume_edge.resize(6, false);
730 for (int ee = 0; ee != 6; ++ee) {
731 N_volume_edge[ee].resize(
732 nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_EDGE_HDIV(volume_order), false);
733 diffN_volume_edge[ee].resize(
734 nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_EDGE_HDIV(volume_order), false);
735 phi_v_e[ee] = &*(N_volume_edge[ee].data().begin());
736 diff_phi_v_e[ee] = &*(diffN_volume_edge[ee].data().begin());
737 }
738 CHKERR Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET(
739 volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
740 &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_v_e,
741 diff_phi_v_e, nb_gauss_pts, base_polynomials);
742
743 N_volume_face.resize(4, false);
744 diffN_volume_face.resize(4, false);
745 for (int ff = 0; ff != 4; ++ff) {
746 N_volume_face[ff].resize(
747 nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order), false);
748 diffN_volume_face[ff].resize(
749 nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order), false);
750 phi_v_f[ff] = &*(N_volume_face[ff].data().begin());
751 diff_phi_v_f[ff] = &*(diffN_volume_face[ff].data().begin());
752 }
753 CHKERR Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET(
754 volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
755 &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_v_f,
756 diff_phi_v_f, nb_gauss_pts, base_polynomials);
757
758 N_volume_bubble.resize(
759 nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order), false);
760 diffN_volume_bubble.resize(
761 nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order), false);
762 phi_v = &*(N_volume_bubble.data().begin());
763 diff_phi_v = &*(diffN_volume_bubble.data().begin());
764 CHKERR Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET(
765 volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
766 &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_v, diff_phi_v,
767 nb_gauss_pts, base_polynomials);
768
769 // Set shape functions into data structure Shape functions hast to be put
770 // in arrays in order which guarantee hierarchical series of degrees of
771 // freedom, i.e. in other words dofs form sub-entities has to be group
772 // by order.
773
774 FTensor::Index<'i', 3> i;
775 FTensor::Index<'j', 3> j;
776
777 // faces
778 if (data.dataOnEntities[MBTRI].size() != 4) {
779 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
780 }
781 for (int ff = 0; ff != 4; ff++) {
782 data.dataOnEntities[MBTRI][ff].getN(base).resize(
783 nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_HDIV(faces_order[ff]), false);
784 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
785 nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_HDIV(faces_order[ff]), false);
786 if (NBFACETRI_AINSWORTH_HDIV(faces_order[ff]) == 0)
787 continue;
788 // face
789 double *base_ptr =
790 &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
791 FTensor::Tensor1<double *, 3> t_base(base_ptr, &base_ptr[HVEC1],
792 &base_ptr[HVEC2], 3);
793 double *diff_base_ptr =
794 &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
795 FTensor::Tensor2<double *, 3, 3> t_diff_base(
796 &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
797 &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
798 &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
799 &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
800 &diff_base_ptr[HVEC2_2], 9);
801 // face-face
802 boost::shared_ptr<FTensor::Tensor1<double *, 3>> t_base_f;
803 boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>> t_diff_base_f;
804 if (NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]) > 0) {
805 base_ptr = phi_f[ff];
806 t_base_f = boost::shared_ptr<FTensor::Tensor1<double *, 3>>(
807 new FTensor::Tensor1<double *, 3>(base_ptr, &base_ptr[HVEC1],
808 &base_ptr[HVEC2], 3));
809 diff_base_ptr = diff_phi_f[ff];
810 t_diff_base_f = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>>(
811 new FTensor::Tensor2<double *, 3, 3>(
812 &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
813 &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
814 &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
815 &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
816 &diff_base_ptr[HVEC2_2], 9));
817 }
818 // edge-face
819 base_ptr = phi_f_e[ff][0];
820 FTensor::Tensor1<double *, 3> t_base_f_e0(base_ptr, &base_ptr[HVEC1],
821 &base_ptr[HVEC2], 3);
822 diff_base_ptr = diff_phi_f_e[ff][0];
823 FTensor::Tensor2<double *, 3, 3> t_diff_base_f_e0(
824 &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
825 &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
826 &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
827 &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
828 &diff_base_ptr[HVEC2_2], 9);
829 base_ptr = phi_f_e[ff][1];
830 FTensor::Tensor1<double *, 3> t_base_f_e1(base_ptr, &base_ptr[HVEC1],
831 &base_ptr[HVEC2], 3);
832 diff_base_ptr = diff_phi_f_e[ff][1];
833 FTensor::Tensor2<double *, 3, 3> t_diff_base_f_e1(
834 &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
835 &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
836 &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
837 &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
838 &diff_base_ptr[HVEC2_2], 9);
839 base_ptr = phi_f_e[ff][2];
840 FTensor::Tensor1<double *, 3> t_base_f_e2(base_ptr, &base_ptr[HVEC1],
841 &base_ptr[HVEC2], 3);
842 diff_base_ptr = diff_phi_f_e[ff][2];
843 FTensor::Tensor2<double *, 3, 3> t_diff_base_f_e2(
844 &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
845 &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
846 &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
847 &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
848 &diff_base_ptr[HVEC2_2], 9);
849 for (int gg = 0; gg != nb_gauss_pts; gg++) {
850 for (int oo = 0; oo != faces_order[ff]; oo++) {
851 for (int dd = NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
852 dd != NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
853 t_base(i) = t_base_f_e0(i);
854 ++t_base;
855 ++t_base_f_e0;
856 t_diff_base(i, j) = t_diff_base_f_e0(i, j);
857 ++t_diff_base;
858 ++t_diff_base_f_e0;
859 t_base(i) = t_base_f_e1(i);
860 ++t_base;
861 ++t_base_f_e1;
862 t_diff_base(i, j) = t_diff_base_f_e1(i, j);
863 ++t_diff_base;
864 ++t_diff_base_f_e1;
865 t_base(i) = t_base_f_e2(i);
866 ++t_base;
867 ++t_base_f_e2;
868 t_diff_base(i, j) = t_diff_base_f_e2(i, j);
869 ++t_diff_base;
870 ++t_diff_base_f_e2;
871 }
872 for (int dd = NBFACETRI_AINSWORTH_FACE_HDIV(oo);
873 dd != NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
874 t_base(i) = (*t_base_f)(i);
875 ++t_base;
876 ++(*t_base_f);
877 t_diff_base(i, j) = (*t_diff_base_f)(i, j);
878 ++t_diff_base;
879 ++(*t_diff_base_f);
880 }
881 }
882 }
883 }
884
885 // volume
886 data.dataOnEntities[MBTET][0].getN(base).resize(
887 nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_HDIV(volume_order), false);
888 data.dataOnEntities[MBTET][0].getDiffN(base).resize(
889 nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_HDIV(volume_order), false);
890 if (NBVOLUMETET_AINSWORTH_HDIV(volume_order) > 0) {
891 double *base_ptr =
892 &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
893 FTensor::Tensor1<double *, 3> t_base(base_ptr, &base_ptr[HVEC1],
894 &base_ptr[HVEC2], 3);
895 double *diff_base_ptr =
896 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
897 FTensor::Tensor2<double *, 3, 3> t_diff_base(
898 &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
899 &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
900 &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
901 &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
902 &diff_base_ptr[HVEC2_2], 9);
903 // edges
904 std::vector<FTensor::Tensor1<double *, 3>> t_base_v_e;
905 t_base_v_e.reserve(6);
906 std::vector<FTensor::Tensor2<double *, 3, 3>> t_diff_base_v_e;
907 t_diff_base_v_e.reserve(6);
908 for (int ee = 0; ee != 6; ee++) {
909 base_ptr = phi_v_e[ee];
910 diff_base_ptr = diff_phi_v_e[ee];
911 t_base_v_e.push_back(FTensor::Tensor1<double *, 3>(
912 base_ptr, &base_ptr[HVEC1], &base_ptr[HVEC2], 3));
913 t_diff_base_v_e.push_back(FTensor::Tensor2<double *, 3, 3>(
914 &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
915 &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
916 &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
917 &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
918 &diff_base_ptr[HVEC2_2], 9));
919 }
920 // faces
921 std::vector<FTensor::Tensor1<double *, 3>> t_base_v_f;
922 t_base_v_f.reserve(4);
923 std::vector<FTensor::Tensor2<double *, 3, 3>> t_diff_base_v_f;
924 t_diff_base_v_f.reserve(4);
925 if (NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order) > 0) {
926 for (int ff = 0; ff != 4; ff++) {
927 base_ptr = phi_v_f[ff];
928 diff_base_ptr = diff_phi_v_f[ff];
929 t_base_v_f.push_back(FTensor::Tensor1<double *, 3>(
930 base_ptr, &base_ptr[HVEC1], &base_ptr[HVEC2], 3));
931 t_diff_base_v_f.push_back(FTensor::Tensor2<double *, 3, 3>(
932 &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
933 &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
934 &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
935 &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
936 &diff_base_ptr[HVEC2_2], 9));
937 }
938 }
939 boost::shared_ptr<FTensor::Tensor1<double *, 3>> t_base_v;
940 boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>> t_diff_base_v;
941 if (NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order) > 0) {
942 base_ptr = phi_v;
943 t_base_v = boost::shared_ptr<FTensor::Tensor1<double *, 3>>(
944 new FTensor::Tensor1<double *, 3>(base_ptr, &base_ptr[HVEC1],
945 &base_ptr[HVEC2], 3));
946 diff_base_ptr = diff_phi_v;
947 t_diff_base_v = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>>(
948 new FTensor::Tensor2<double *, 3, 3>(
949 &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
950 &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
951 &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
952 &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
953 &diff_base_ptr[HVEC2_2], 9));
954 }
955 for (int gg = 0; gg != nb_gauss_pts; gg++) {
956 for (int oo = 0; oo < volume_order; oo++) {
957 for (int dd = NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo);
958 dd < NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
959 for (int ee = 0; ee < 6; ee++) {
960 t_base(i) = t_base_v_e[ee](i);
961 ++t_base;
962 ++t_base_v_e[ee];
963 t_diff_base(i, j) = t_diff_base_v_e[ee](i, j);
964 ++t_diff_base;
965 ++t_diff_base_v_e[ee];
966 }
967 }
968 for (int dd = NBVOLUMETET_AINSWORTH_FACE_HDIV(oo);
969 dd < NBVOLUMETET_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
970 for (int ff = 0; ff < 4; ff++) {
971 t_base(i) = t_base_v_f[ff](i);
972 ++t_base;
973 ++t_base_v_f[ff];
974 t_diff_base(i, j) = t_diff_base_v_f[ff](i, j);
975 ++t_diff_base;
976 ++t_diff_base_v_f[ff];
977 }
978 }
979 for (int dd = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo);
980 dd < NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo + 1); dd++) {
981 t_base(i) = (*t_base_v)(i);
982 ++t_base;
983 ++(*t_base_v);
984 t_diff_base(i, j) = (*t_diff_base_v)(i, j);
985 ++t_diff_base;
986 ++(*t_diff_base_v);
987 }
988 }
989 }
990 }
991
992 MoFEMFunctionReturn(0);
993}
994
995MoFEMErrorCode TetPolynomialBase::getValueHdivDemkowiczBase(MatrixDouble &pts) {
996 MoFEMFunctionBegin;
997
998 EntitiesFieldData &data = cTx->dAta;
999 const FieldApproximationBase base = cTx->bAse;
1000 if (base != DEMKOWICZ_JACOBI_BASE) {
1001 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1002 "This should be used only with DEMKOWICZ_JACOBI_BASE "
1003 "but base is %s",
1004 ApproximationBaseNames[base]);
1005 }
1006 int nb_gauss_pts = pts.size2();
1007
1008 int volume_order = data.dataOnEntities[MBTET][0].getOrder();
1009
1010 int p_f[4];
1011 double *phi_f[4];
1012 double *diff_phi_f[4];
1013
1014 // Calculate base function on tet faces
1015 for (int ff = 0; ff != 4; ff++) {
1016 int face_order = data.dataOnEntities[MBTRI][ff].getOrder();
1017 int order = volume_order > face_order ? volume_order : face_order;
1018 data.dataOnEntities[MBTRI][ff].getN(base).resize(
1019 nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
1020 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
1021 nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
1022 p_f[ff] = order;
1023 phi_f[ff] = &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1024 diff_phi_f[ff] =
1025 &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1026 if (NBFACETRI_DEMKOWICZ_HDIV(order) == 0)
1027 continue;
1028 CHKERR Hdiv_Demkowicz_Face_MBTET_ON_FACE(
1029 &data.facesNodes(ff, 0), order,
1030 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1031 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1032 phi_f[ff], diff_phi_f[ff], nb_gauss_pts, 4);
1033 }
1034
1035 // Calculate base functions in tet interior
1036 if (NBVOLUMETET_DEMKOWICZ_HDIV(volume_order) > 0) {
1037 data.dataOnEntities[MBTET][0].getN(base).resize(
1038 nb_gauss_pts, 3 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
1039 data.dataOnEntities[MBTET][0].getDiffN(base).resize(
1040 nb_gauss_pts, 9 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
1041 double *phi_v = &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
1042 double *diff_phi_v =
1043 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
1044 CHKERR Hdiv_Demkowicz_Interior_MBTET(
1045 volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
1046 &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), p_f, phi_f,
1047 diff_phi_f, phi_v, diff_phi_v, nb_gauss_pts);
1048 }
1049
1050 // Set size of face base correctly
1051 for (int ff = 0; ff != 4; ff++) {
1052 int face_order = data.dataOnEntities[MBTRI][ff].getOrder();
1053 data.dataOnEntities[MBTRI][ff].getN(base).resize(
1054 nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
1055 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
1056 nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
1057 }
1058
1059 MoFEMFunctionReturn(0);
1060}
1061
1062MoFEMErrorCode TetPolynomialBase::getValueHdiv(MatrixDouble &pts) {
1063 MoFEMFunctionBeginHot;
1064
1065 switch (cTx->bAse) {
1066 case AINSWORTH_LEGENDRE_BASE:
1067 case AINSWORTH_LOBATTO_BASE:
1068 return getValueHdivAinsworthBase(pts);
1069 case DEMKOWICZ_JACOBI_BASE:
1070 return getValueHdivDemkowiczBase(pts);
1071 default:
1072 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1073 }
1074
1075 MoFEMFunctionReturnHot(0);
1076}
1077
1078MoFEMErrorCode
1079TetPolynomialBase::getValueHcurlAinsworthBase(MatrixDouble &pts) {
1080 MoFEMFunctionBegin;
1081
1082 EntitiesFieldData &data = cTx->dAta;
1083 const FieldApproximationBase base = cTx->bAse;
1084 PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
1085 double *diffL, const int dim) =
1086 cTx->basePolynomialsType0;
1087
1088 int nb_gauss_pts = pts.size2();
1089
1090 // edges
1091 if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
1092 int sense[6], order[6];
1093 if (data.dataOnEntities[MBEDGE].size() != 6) {
1094 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1095 }
1096 double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
1097 for (int ee = 0; ee != 6; ee++) {
1098 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
1099 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1100 "data inconsistency");
1101 }
1102 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
1103 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
1104 int nb_dofs = NBEDGE_AINSWORTH_HCURL(
1105 data.dataOnEntities[MBEDGE][ee].getOrder());
1106 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
1107 3 * nb_dofs, false);
1108 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
1109 9 * nb_dofs, false);
1110 hcurl_edge_n[ee] =
1111 &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
1112 diff_hcurl_edge_n[ee] =
1113 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
1114 }
1115 CHKERR Hcurl_Ainsworth_EdgeBaseFunctions_MBTET(
1116 sense, order,
1117 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1118 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1119 hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
1120 } else {
1121 for (int ee = 0; ee != 6; ee++) {
1122 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
1123 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
1124 false);
1125 }
1126 }
1127
1128 // triangles
1129 if (data.spacesOnEntities[MBTRI].test(HCURL)) {
1130 int order[4];
1131 // faces
1132 if (data.dataOnEntities[MBTRI].size() != 4) {
1133 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1134 }
1135 double *hcurl_base_n[4], *diff_hcurl_base_n[4];
1136 for (int ff = 0; ff != 4; ff++) {
1137 if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
1138 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1139 "data inconsistency");
1140 }
1141 order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
1142 int nb_dofs = NBFACETRI_AINSWORTH_HCURL(order[ff]);
1143 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1144 3 * nb_dofs, false);
1145 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1146 9 * nb_dofs, false);
1147 hcurl_base_n[ff] =
1148 &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1149 diff_hcurl_base_n[ff] =
1150 &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1151 }
1152 if (data.facesNodes.size1() != 4) {
1153 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1154 }
1155 if (data.facesNodes.size2() != 3) {
1156 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1157 }
1158 CHKERR Hcurl_Ainsworth_FaceFunctions_MBTET(
1159 &*data.facesNodes.data().begin(), order,
1160 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1161 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1162 hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts, base_polynomials);
1163 } else {
1164 for (int ff = 0; ff != 4; ff++) {
1165 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
1166 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
1167 false);
1168 }
1169 }
1170
1171 if (data.spacesOnEntities[MBTET].test(HCURL)) {
1172
1173 // volume
1174 int order = data.dataOnEntities[MBTET][0].getOrder();
1175 int nb_vol_dofs = NBVOLUMETET_AINSWORTH_HCURL(order);
1176 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
1177 3 * nb_vol_dofs, false);
1178 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
1179 9 * nb_vol_dofs, false);
1180 CHKERR Hcurl_Ainsworth_VolumeFunctions_MBTET(
1181 data.dataOnEntities[MBTET][0].getOrder(),
1182 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1183 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1184 &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
1185 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
1186 nb_gauss_pts, base_polynomials);
1187
1188 } else {
1189 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
1190 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
1191 }
1192
1193 MoFEMFunctionReturn(0);
1194}
1195
1196MoFEMErrorCode
1197TetPolynomialBase::getValueHcurlDemkowiczBase(MatrixDouble &pts) {
1198 MoFEMFunctionBegin;
1199
1200 EntitiesFieldData &data = cTx->dAta;
1201 const FieldApproximationBase base = cTx->bAse;
1202 if (base != DEMKOWICZ_JACOBI_BASE) {
1203 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1204 "This should be used only with DEMKOWICZ_JACOBI_BASE "
1205 "but base is %s",
1206 ApproximationBaseNames[base]);
1207 }
1208
1209 int nb_gauss_pts = pts.size2();
1210
1211 // edges
1212 if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
1213 int sense[6], order[6];
1214 if (data.dataOnEntities[MBEDGE].size() != 6) {
1215 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1216 "wrong size of data structure, expected space for six edges "
1217 "but is %d",
1218 data.dataOnEntities[MBEDGE].size());
1219 }
1220 double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
1221 for (int ee = 0; ee != 6; ee++) {
1222 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
1223 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1224 "orintation of edges is not set");
1225 }
1226 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
1227 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
1228 int nb_dofs = NBEDGE_DEMKOWICZ_HCURL(
1229 data.dataOnEntities[MBEDGE][ee].getOrder());
1230 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
1231 3 * nb_dofs, false);
1232 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
1233 9 * nb_dofs, false);
1234 hcurl_edge_n[ee] =
1235 &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
1236 diff_hcurl_edge_n[ee] =
1237 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
1238 }
1239 CHKERR Hcurl_Demkowicz_EdgeBaseFunctions_MBTET(
1240 sense, order,
1241 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1242 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1243 hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts);
1244 } else {
1245 // No DOFs on edges, resize base function matrices, indicating that no
1246 // dofs on them.
1247 for (int ee = 0; ee != 6; ee++) {
1248 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
1249 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
1250 false);
1251 }
1252 }
1253
1254 // triangles
1255 if (data.spacesOnEntities[MBTRI].test(HCURL)) {
1256 int order[4];
1257 // faces
1258 if (data.dataOnEntities[MBTRI].size() != 4) {
1259 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1260 "data structure for storing face h-curl base have wrong size "
1261 "should be four but is %d",
1262 data.dataOnEntities[MBTRI].size());
1263 }
1264 double *hcurl_base_n[4], *diff_hcurl_base_n[4];
1265 for (int ff = 0; ff != 4; ff++) {
1266 if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
1267 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1268 "orintation of face is not set");
1269 }
1270 order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
1271 int nb_dofs = NBFACETRI_DEMKOWICZ_HCURL(order[ff]);
1272 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1273 3 * nb_dofs, false);
1274 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1275 9 * nb_dofs, false);
1276 hcurl_base_n[ff] =
1277 &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1278 diff_hcurl_base_n[ff] =
1279 &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1280 }
1281 if (data.facesNodes.size1() != 4) {
1282 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1283 "data inconsistency, should be four faces");
1284 }
1285 if (data.facesNodes.size2() != 3) {
1286 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1287 "data inconsistency, should be three nodes on face");
1288 }
1289 CHKERR Hcurl_Demkowicz_FaceBaseFunctions_MBTET(
1290 &*data.facesNodes.data().begin(), order,
1291 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1292 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1293 hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts);
1294 } else {
1295 // No DOFs on faces, resize base function matrices, indicating that no
1296 // dofs on them.
1297 for (int ff = 0; ff != 4; ff++) {
1298 data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
1299 data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
1300 false);
1301 }
1302 }
1303
1304 if (data.spacesOnEntities[MBTET].test(HCURL)) {
1305 // volume
1306 int order = data.dataOnEntities[MBTET][0].getOrder();
1307 int nb_vol_dofs = NBVOLUMETET_DEMKOWICZ_HCURL(order);
1308 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
1309 3 * nb_vol_dofs, false);
1310 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
1311 9 * nb_vol_dofs, false);
1312 CHKERR Hcurl_Demkowicz_VolumeBaseFunctions_MBTET(
1313 data.dataOnEntities[MBTET][0].getOrder(),
1314 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1315 &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1316 &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
1317 &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
1318 nb_gauss_pts);
1319 } else {
1320 // No DOFs on faces, resize base function matrices, indicating that no
1321 // dofs on them.
1322 data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
1323 data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
1324 }
1325
1326 MoFEMFunctionReturn(0);
1327}
1328
1329MoFEMErrorCode TetPolynomialBase::getValueHcurl(MatrixDouble &pts) {
1330 MoFEMFunctionBegin;
1331
1332 switch (cTx->bAse) {
1333 case AINSWORTH_LEGENDRE_BASE:
1334 case AINSWORTH_LOBATTO_BASE:
1335 CHKERR getValueHcurlAinsworthBase(pts);
1336 break;
1337 case DEMKOWICZ_JACOBI_BASE:
1338 CHKERR getValueHcurlDemkowiczBase(pts);
1339 break;
1340 default:
1341 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1342 }
1343
1344 MoFEMFunctionReturn(0);
1345}
1346
1347MoFEMErrorCode
1348TetPolynomialBase::getValue(MatrixDouble &pts,
1349 boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
1350 MoFEMFunctionBegin;
1351
1352 cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
1353
1354 int nb_gauss_pts = pts.size2();
1355 if (!nb_gauss_pts)
1356 MoFEMFunctionReturnHot(0);
1357
1358 if (pts.size1() < 3)
1359 SETERRQ(
1360 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1361 "Wrong dimension of pts, should be at least 3 rows with coordinates");
1362
1363 if (cTx->bAse != AINSWORTH_BERNSTEIN_BEZIER_BASE) {
1364 const FieldApproximationBase base = cTx->bAse;
1365 EntitiesFieldData &data = cTx->dAta;
1366 if (cTx->copyNodeBase == LASTBASE) {
1367 data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4,
1368 false);
1369 CHKERR Tools::shapeFunMBTET(
1370 &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1371 &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
1372 } else {
1373 data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
1374 data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
1375 }
1376 if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
1377 (unsigned int)nb_gauss_pts) {
1378 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1379 "Base functions or nodes has wrong number of integration points "
1380 "for base %s",
1381 ApproximationBaseNames[base]);
1382 }
1383 data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, false);
1384 std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
1385 data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
1386 }
1387
1388 switch (cTx->sPace) {
1389 case H1:
1390 CHKERR getValueH1(pts);
1391 break;
1392 case HDIV:
1393 CHKERR getValueHdiv(pts);
1394 break;
1395 case HCURL:
1396 CHKERR getValueHcurl(pts);
1397 break;
1398 case L2:
1399 CHKERR getValueL2(pts);
1400 break;
1401 default:
1402 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space");
1403 }
1404
1405 MoFEMFunctionReturn(0);
1406}
o
static Index< 'o', 3 > o
Definition: BasicFeTools.hpp:89
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:90
a
constexpr double a
Definition: approx_sphere.cpp:30
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
LASTBASE
@ LASTBASE
Definition: definitions.h:69
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
NOBASE
@ NOBASE
Definition: definitions.h:59
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
H1
@ H1
continuous field
Definition: definitions.h:85
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
HVEC1
@ HVEC1
Definition: definitions.h:186
HVEC2
@ HVEC2
Definition: definitions.h:186
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
ApproximationBaseNames
static const char *const ApproximationBaseNames[]
Definition: definitions.h:72
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
HVEC1_1
@ HVEC1_1
Definition: definitions.h:196
HVEC0_1
@ HVEC0_1
Definition: definitions.h:195
HVEC1_0
@ HVEC1_0
Definition: definitions.h:193
HVEC2_1
@ HVEC2_1
Definition: definitions.h:197
HVEC1_2
@ HVEC1_2
Definition: definitions.h:199
HVEC2_2
@ HVEC2_2
Definition: definitions.h:200
HVEC2_0
@ HVEC2_0
Definition: definitions.h:194
HVEC0_2
@ HVEC0_2
Definition: definitions.h:198
HVEC0_0
@ HVEC0_0
Definition: definitions.h:192
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
dim
const int dim
Definition: elec_phys_2D.cpp:12
n
FTensor::Index< 'n', SPACE_DIM > n
Definition: fluid_structure_eigenproblem.cpp:66
order
int order
Definition: fluid_structure_eigenproblem.cpp:84
lambda
constexpr double lambda
surface tension
Definition: free_surface.cpp:165
NBVOLUMETET_AINSWORTH_EDGE_HDIV
#define NBVOLUMETET_AINSWORTH_EDGE_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:133
NBEDGE_DEMKOWICZ_HCURL
#define NBEDGE_DEMKOWICZ_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:108
NBVOLUMETET_H1
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
Definition: h1_hdiv_hcurl_l2.h:75
NBVOLUMETET_AINSWORTH_HDIV
#define NBVOLUMETET_AINSWORTH_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:137
NBVOLUMETET_AINSWORTH_HCURL
#define NBVOLUMETET_AINSWORTH_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:105
NBFACETRI_AINSWORTH_HCURL
#define NBFACETRI_AINSWORTH_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:100
NBVOLUMETET_DEMKOWICZ_HDIV
#define NBVOLUMETET_DEMKOWICZ_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:140
H1_EdgeShapeFunctions_MBTET
PetscErrorCode H1_EdgeShapeFunctions_MBTET(int *sense, int *p, double *N, double *diffN, double *edgeN[], double *diff_edgeN[], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: h1.c:274
NBVOLUMETET_AINSWORTH_FACE_HDIV
#define NBVOLUMETET_AINSWORTH_FACE_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:134
NBEDGE_H1
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
Definition: h1_hdiv_hcurl_l2.h:55
NBFACETRI_DEMKOWICZ_HDIV
#define NBFACETRI_DEMKOWICZ_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:139
L2_Ainsworth_ShapeFunctions_MBTET
PetscErrorCode L2_Ainsworth_ShapeFunctions_MBTET(int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Get base functions on tetrahedron for L2 space.
Definition: l2.c:74
NBVOLUMETET_DEMKOWICZ_HCURL
#define NBVOLUMETET_DEMKOWICZ_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:110
NBFACETRI_DEMKOWICZ_HCURL
#define NBFACETRI_DEMKOWICZ_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:109
NBFACETRI_AINSWORTH_FACE_HDIV
#define NBFACETRI_AINSWORTH_FACE_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:131
NBVOLUMETET_AINSWORTH_VOLUME_HDIV
#define NBVOLUMETET_AINSWORTH_VOLUME_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:135
NBFACETRI_AINSWORTH_EDGE_HDIV
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:130
H1_VolumeShapeFunctions_MBTET
PetscErrorCode H1_VolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *volumeN, double *diff_volumeN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: h1.c:475
NBEDGE_AINSWORTH_HCURL
#define NBEDGE_AINSWORTH_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:97
H1_FaceShapeFunctions_MBTET
PetscErrorCode H1_FaceShapeFunctions_MBTET(int *faces_nodes, int *p, double *N, double *diffN, double *faceN[], double *diff_faceN[], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: h1.c:373
NBFACETRI_H1
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
Definition: h1_hdiv_hcurl_l2.h:60
NBFACETRI_AINSWORTH_HDIV
#define NBFACETRI_AINSWORTH_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:132
NBVOLUMETET_L2
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
Definition: h1_hdiv_hcurl_l2.h:27
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::Types::MatrixInt
UBlasMatrix< int > MatrixInt
Definition: Types.hpp:76
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::Hdiv_Ainsworth_FaceBubbleShapeFunctions
MoFEMErrorCode Hdiv_Ainsworth_FaceBubbleShapeFunctions(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f[], double *diff_phi_f[], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face bubble functions by Ainsworth .
Definition: Hdiv.cpp:137
MoFEM::Hdiv_Demkowicz_Face_MBTET_ON_FACE
MoFEMErrorCode Hdiv_Demkowicz_Face_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb)
Definition: Hdiv.cpp:617
MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET
MoFEMErrorCode Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f_e[4][3], double *diff_phi_f_e[4][3], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base functions, Edge-based face functions by Ainsworth .
Definition: Hdiv.cpp:13
MoFEM::Hdiv_Demkowicz_Interior_MBTET
MoFEMErrorCode Hdiv_Demkowicz_Interior_MBTET(int p, double *N, double *diffN, int p_face[], double *phi_f[4], double *diff_phi_f[4], double *phi_v, double *diff_phi_v, int gdim)
Definition: Hdiv.cpp:762
MoFEM::Hcurl_Ainsworth_VolumeFunctions_MBTET
MoFEMErrorCode Hcurl_Ainsworth_VolumeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v, double *diff_phi_v, int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
H-curl volume base functions.
Definition: Hcurl.cpp:1403
MoFEM::Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET
MoFEMErrorCode Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_f[], double *diff_phi_v_f[], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: Hdiv.cpp:368
MoFEM::Hcurl_Ainsworth_EdgeBaseFunctions_MBTET
MoFEMErrorCode Hcurl_Ainsworth_EdgeBaseFunctions_MBTET(int *sense, int *p, double *N, double *diffN, double *edgeN[], double *diff_edgeN[], int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Edge based H-curl base functions on tetrahedral.
Definition: Hcurl.cpp:16
MoFEM::Hcurl_Demkowicz_FaceBaseFunctions_MBTET
MoFEMErrorCode Hcurl_Demkowicz_FaceBaseFunctions_MBTET(int *faces_nodes, int *p, double *n, double *diff_n, double *phi[], double *diff_phi[], int nb_integration_pts)
Face base interior function.
Definition: Hcurl.cpp:2395
MoFEM::Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET
MoFEMErrorCode Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_e[6], double *diff_phi_v_e[6], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base function, Edge-based interior (volume) functions by Ainsworth .
Definition: Hdiv.cpp:280
MoFEM::Hcurl_Demkowicz_VolumeBaseFunctions_MBTET
MoFEMErrorCode Hcurl_Demkowicz_VolumeBaseFunctions_MBTET(int p, double *n, double *diff_n, double *phi, double *diff_phi, int nb_integration_pts)
Volume base interior function.
Definition: Hcurl.cpp:2468
MoFEM::Hcurl_Demkowicz_EdgeBaseFunctions_MBTET
MoFEMErrorCode Hcurl_Demkowicz_EdgeBaseFunctions_MBTET(int *sense, int *p, double *n, double *diff_n, double *phi[], double *diff_phi[], int nb_integration_pts)
Edge based H-curl base functions on tetrahedral.
Definition: Hcurl.cpp:2068
MoFEM::Hcurl_Ainsworth_FaceFunctions_MBTET
MoFEMErrorCode Hcurl_Ainsworth_FaceFunctions_MBTET(int *face_nodes, int *p, double *N, double *diffN, double *phi_f[4], double *diff_phi_f[4], int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face H-curl functions.
Definition: Hcurl.cpp:1052
MoFEM::Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET
MoFEMErrorCode Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v, double *diff_phi_v, int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Interior bubble functions by Ainsworth .
Definition: Hdiv.cpp:484
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
g
constexpr double g
Definition: shallow_wave.cpp:63
MoFEM::BernsteinBezier::generateIndicesTriTet
static MoFEMErrorCode generateIndicesTriTet(const int N[], int *alpha[])
Definition: BernsteinBezier.cpp:174
MoFEM::BernsteinBezier::baseFunctionsTet
static MoFEMErrorCode baseFunctionsTet(const int N, const int gdim, const int n_alpha, const int *alpha, const double *lambda, const double *grad_lambda, double *base, double *grad_base)
Definition: BernsteinBezier.cpp:443
MoFEM::BernsteinBezier::generateIndicesEdgeTet
static MoFEMErrorCode generateIndicesEdgeTet(const int N[], int *alpha[])
Definition: BernsteinBezier.cpp:138
MoFEM::BernsteinBezier::generateIndicesVertexTet
static MoFEMErrorCode generateIndicesVertexTet(const int N, int *alpha)
Definition: BernsteinBezier.cpp:128
MoFEM::BernsteinBezier::generateIndicesTetTet
static MoFEMErrorCode generateIndicesTetTet(const int N, int *alpha)
Definition: BernsteinBezier.cpp:203
MoFEM::EntPolynomialBaseCtx
Class used to pass element data to calculate base functions on tet,triangle,edge.
Definition: EntPolynomialBaseCtx.hpp:22
MoFEM::EntPolynomialBaseCtx::copyNodeBase
const FieldApproximationBase copyNodeBase
Definition: EntPolynomialBaseCtx.hpp:40
MoFEM::EntPolynomialBaseCtx::basePolynomialsType0
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Definition: EntPolynomialBaseCtx.hpp:27
MoFEM::EntPolynomialBaseCtx::bAse
const FieldApproximationBase bAse
Definition: EntPolynomialBaseCtx.hpp:38
MoFEM::EntitiesFieldData
data structure for finite element entity
Definition: EntitiesFieldData.hpp:40
MoFEM::EntitiesFieldData::facesNodes
MatrixInt facesNodes
nodes on finite element faces
Definition: EntitiesFieldData.hpp:46
MoFEM::EntitiesFieldData::spacesOnEntities
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
Definition: EntitiesFieldData.hpp:50
MoFEM::EntitiesFieldData::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition: EntitiesFieldData.hpp:56
MoFEM::TetPolynomialBase
Calculate base functions on tetrahedral.
Definition: TetPolynomialBase.hpp:19
MoFEM::TetPolynomialBase::diffN_face_edge
ublas::matrix< MatrixDouble > diffN_face_edge
Definition: TetPolynomialBase.hpp:83
MoFEM::TetPolynomialBase::cTx
EntPolynomialBaseCtx * cTx
Definition: TetPolynomialBase.hpp:31
MoFEM::TetPolynomialBase::getValueH1AinsworthBase
MoFEMErrorCode getValueH1AinsworthBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:37
MoFEM::TetPolynomialBase::diffN_volume_face
ublas::vector< MatrixDouble > diffN_volume_face
Definition: TetPolynomialBase.hpp:86
MoFEM::TetPolynomialBase::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition: TetPolynomialBase.cpp:11
MoFEM::TetPolynomialBase::getValueL2BernsteinBezierBase
MoFEMErrorCode getValueL2BernsteinBezierBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:479
MoFEM::TetPolynomialBase::getValueHcurlDemkowiczBase
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:1197
MoFEM::TetPolynomialBase::getValueHcurl
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
Get base functions for Hcurl space.
Definition: TetPolynomialBase.cpp:1329
MoFEM::TetPolynomialBase::getValueH1BernsteinBezierBase
MoFEMErrorCode getValueH1BernsteinBezierBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:149
MoFEM::TetPolynomialBase::getValueH1
MoFEMErrorCode getValueH1(MatrixDouble &pts)
Get base functions for H1 space.
Definition: TetPolynomialBase.cpp:19
MoFEM::TetPolynomialBase::diffN_face_bubble
ublas::vector< MatrixDouble > diffN_face_bubble
Definition: TetPolynomialBase.hpp:84
MoFEM::TetPolynomialBase::N_volume_edge
ublas::vector< MatrixDouble > N_volume_edge
Definition: TetPolynomialBase.hpp:79
MoFEM::TetPolynomialBase::N_volume_face
ublas::vector< MatrixDouble > N_volume_face
Definition: TetPolynomialBase.hpp:80
MoFEM::TetPolynomialBase::getValueHdiv
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
Get base functions for Hdiv space.
Definition: TetPolynomialBase.cpp:1062
MoFEM::TetPolynomialBase::getValueHdivAinsworthBase
MoFEMErrorCode getValueHdivAinsworthBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:655
MoFEM::TetPolynomialBase::getValueL2AinsworthBase
MoFEMErrorCode getValueL2AinsworthBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:449
MoFEM::TetPolynomialBase::N_face_bubble
ublas::vector< MatrixDouble > N_face_bubble
Definition: TetPolynomialBase.hpp:78
MoFEM::TetPolynomialBase::N_face_edge
ublas::matrix< MatrixDouble > N_face_edge
Definition: TetPolynomialBase.hpp:77
MoFEM::TetPolynomialBase::getValueL2
MoFEMErrorCode getValueL2(MatrixDouble &pts)
Get base functions for L2 space.
Definition: TetPolynomialBase.cpp:431
MoFEM::TetPolynomialBase::getValueHdivDemkowiczBase
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:995
MoFEM::TetPolynomialBase::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: TetPolynomialBase.cpp:1348
MoFEM::TetPolynomialBase::diffN_volume_edge
ublas::vector< MatrixDouble > diffN_volume_edge
Definition: TetPolynomialBase.hpp:85
MoFEM::TetPolynomialBase::getValueHcurlAinsworthBase
MoFEMErrorCode getValueHcurlAinsworthBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:1079
MoFEM::Tools::shapeFunMBTET
static MoFEMErrorCode shapeFunMBTET(double *shape, const double *ksi, const double *eta, const double *zeta, const double nb)
Calculate shape functions on tetrahedron.
Definition: Tools.hpp:673
MoFEM::Tools::diffShapeFunMBTET
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition: Tools.hpp:271
MoFEM::UnknownInterface
base class for all interface classes
Definition: UnknownInterface.hpp:34
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
Generated by Doxygen 1.9.5 and hosted at University of Glasgow