v0.14.0
Loading...
Searching...
No Matches
HexPolynomialBase.cpp
Go to the documentation of this file.
1/** \file HexPolynomialBase.cpp
2\brief Implementation of hierarchical bases on tetrahedral
3
4A l2, h1, h-div and h-curl spaces are implemented.
5
6*/
7
8
9
10using namespace MoFEM;
11
12MoFEMErrorCode
13HexPolynomialBase::query_interface(boost::typeindex::type_index type_index,
14 UnknownInterface **iface) const {
15 MoFEMFunctionBegin;
16 *iface = const_cast<HexPolynomialBase *>(this);
17 MoFEMFunctionReturn(0);
18}
19
20MoFEMErrorCode HexPolynomialBase::getValueH1(MatrixDouble &pts) {
21 MoFEMFunctionBegin;
22
23 switch (cTx->bAse) {
24 case DEMKOWICZ_JACOBI_BASE:
25 CHKERR getValueH1DemkowiczBase(pts);
26 break;
27 default:
28 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
29 }
30
31 MoFEMFunctionReturn(0);
32}
33
34MoFEMErrorCode HexPolynomialBase::getValueH1DemkowiczBase(MatrixDouble &pts) {
35 MoFEMFunctionBegin;
36
37 auto &data = cTx->dAta;
38 const auto base = cTx->bAse;
39 const auto copy_base = cTx->copyNodeBase;
40 int nb_gauss_pts = pts.size2();
41
42 auto &copy_base_fun = data.dataOnEntities[MBVERTEX][0].getN(copy_base);
43 auto &copy_diff_base_fun =
44 data.dataOnEntities[MBVERTEX][0].getDiffN(copy_base);
45 if (copy_base_fun.size1() != nb_gauss_pts)
46 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
47 "Inconsistent number of integration pts");
48
49 auto add_base_on_verts = [&] {
50 MoFEMFunctionBeginHot;
51 data.dataOnEntities[MBVERTEX][0].getN(base).resize(
52 nb_gauss_pts, copy_base_fun.size2(), false);
53 data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(
54 nb_gauss_pts, copy_diff_base_fun.size2(), false);
55 noalias(data.dataOnEntities[MBVERTEX][0].getN(base)) = copy_base_fun;
56 noalias(data.dataOnEntities[MBVERTEX][0].getDiffN(base)) =
57 copy_diff_base_fun;
58 MoFEMFunctionReturnHot(0);
59 };
60
61 // Edges
62 auto add_base_on_edges = [&] {
63 MoFEMFunctionBeginHot;
64 std::array<int, 12> sense;
65 std::array<int, 12> order;
66 if (data.spacesOnEntities[MBEDGE].test(H1)) {
67 if (data.dataOnEntities[MBEDGE].size() != 12)
68 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
69 "Expected 12 data on entities");
70
71 std::array<double *, 12> h1_edge_n;
72 std::array<double *, 12> diff_h1_egde_n;
73 bool nb_dofs_sum = false;
74 for (int ee = 0; ee != 12; ++ee) {
75 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0)
76 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
77 "Sense of entity not set");
78
79 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
80 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
81
82 int nb_dofs = NBEDGE_H1(data.dataOnEntities[MBEDGE][ee].getOrder());
83 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
84 false);
85 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(
86 nb_gauss_pts, 3 * nb_dofs, false);
87 h1_edge_n[ee] =
88 &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
89 diff_h1_egde_n[ee] =
90 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
91
92 nb_dofs_sum |= (nb_dofs > 0);
93 }
94 if (nb_dofs_sum) {
95 CHKERR DemkowiczHexAndQuad::H1_EdgeShapeFunctions_ONHEX(
96 sense.data(), order.data(), &*copy_base_fun.data().begin(),
97 &*copy_diff_base_fun.data().begin(), h1_edge_n.data(),
98 diff_h1_egde_n.data(), nb_gauss_pts);
99 }
100 } else {
101 for (int ee = 0; ee != 12; ++ee) {
102 data.dataOnEntities[MBEDGE][ee].getN(base).resize(0, 0, false);
103 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(0, 0, false);
104 }
105 }
106 MoFEMFunctionReturnHot(0);
107 };
108
109 // Face
110 auto add_base_on_quads = [&]() {
111 MoFEMFunctionBeginHot;
112 std::array<int, 6> order;
113 if (data.spacesOnEntities[MBQUAD].test(H1)) {
114 // faces
115 if (data.dataOnEntities[MBQUAD].size() != 6)
116 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
117 "Expected six faces on hex");
118
119 std::array<double *, 6> h1_face_n;
120 std::array<double *, 6> diff_h1_face_n;
121 bool nb_dofs_sum = false;
122 for (int ff = 0; ff != 6; ++ff) {
123
124 order[ff] = data.dataOnEntities[MBQUAD][ff].getOrder();
125 const int nb_dofs = NBFACEQUAD_H1(order[ff]);
126
127 data.dataOnEntities[MBQUAD][ff].getN(base).resize(nb_gauss_pts, nb_dofs,
128 false);
129 data.dataOnEntities[MBQUAD][ff].getDiffN(base).resize(
130 nb_gauss_pts, 3 * nb_dofs, false);
131
132 h1_face_n[ff] =
133 &*data.dataOnEntities[MBQUAD][ff].getN(base).data().begin();
134 diff_h1_face_n[ff] =
135 &*data.dataOnEntities[MBQUAD][ff].getDiffN(base).data().begin();
136
137 nb_dofs_sum |= (nb_dofs > 0);
138 }
139 if (data.facesNodes.size1() != 6)
140 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
141 "Expected six face nodes");
142 if (data.facesNodes.size2() != 4)
143 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
144 "Expected four nodes on face");
145
146 if (nb_dofs_sum) {
147 CHKERR DemkowiczHexAndQuad::H1_FaceShapeFunctions_ONHEX(
148 &*data.facesNodes.data().begin(),
149 &*data.facesNodesOrder.data().begin(), order.data(),
150 &*copy_base_fun.data().begin(), &*copy_diff_base_fun.data().begin(),
151 h1_face_n.data(), diff_h1_face_n.data(), nb_gauss_pts);
152 }
153
154 } else {
155 for (int ff = 0; ff != 6; ++ff) {
156 data.dataOnEntities[MBQUAD][ff].getN(base).resize(0, false);
157 data.dataOnEntities[MBQUAD][ff].getDiffN(base).resize(0, 0, false);
158 }
159 }
160
161 MoFEMFunctionReturnHot(0);
162 };
163
164 // Face
165 auto add_base_on_volume = [&]() {
166 MoFEMFunctionBeginHot;
167
168 if (data.spacesOnEntities[MBHEX].test(H1)) {
169 // volume
170 int order = data.dataOnEntities[MBHEX][0].getOrder();
171 int nb_vol_dofs = NBVOLUMEHEX_H1(order);
172 data.dataOnEntities[MBHEX][0].getN(base).resize(nb_gauss_pts, nb_vol_dofs,
173 false);
174 data.dataOnEntities[MBHEX][0].getDiffN(base).resize(
175 nb_gauss_pts, 3 * nb_vol_dofs, false);
176
177 if (nb_vol_dofs) {
178 const std::array<int, 3> p{order, order, order};
179 CHKERR DemkowiczHexAndQuad::H1_InteriorShapeFunctions_ONHEX(
180 p.data(), &*copy_base_fun.data().begin(),
181 &*copy_diff_base_fun.data().begin(),
182 &*data.dataOnEntities[MBHEX][0].getN(base).data().begin(),
183 &*data.dataOnEntities[MBHEX][0].getDiffN(base).data().begin(),
184 nb_gauss_pts);
185 }
186 } else {
187 data.dataOnEntities[MBHEX][0].getN(base).resize(0, 0, false);
188 data.dataOnEntities[MBHEX][0].getDiffN(base).resize(0, 0, false);
189 }
190
191 MoFEMFunctionReturnHot(0);
192 };
193
194 CHKERR add_base_on_verts();
195 CHKERR add_base_on_edges();
196 CHKERR add_base_on_quads();
197 CHKERR add_base_on_volume();
198
199 MoFEMFunctionReturn(0);
200}
201
202MoFEMErrorCode HexPolynomialBase::getValueL2(MatrixDouble &pts) {
203 MoFEMFunctionBegin;
204
205 switch (cTx->bAse) {
206 case DEMKOWICZ_JACOBI_BASE:
207 CHKERR getValueL2DemkowiczBase(pts);
208 break;
209 default:
210 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
211 }
212
213 MoFEMFunctionReturn(0);
214}
215
216MoFEMErrorCode HexPolynomialBase::getValueL2DemkowiczBase(MatrixDouble &pts) {
217 MoFEMFunctionBegin;
218
219 auto &data = cTx->dAta;
220 const auto base = cTx->bAse;
221 const auto copy_base = cTx->copyNodeBase;
222
223 int nb_gauss_pts = pts.size2();
224
225 auto &copy_base_fun = data.dataOnEntities[MBVERTEX][0].getN(copy_base);
226 auto &copy_diff_base_fun =
227 data.dataOnEntities[MBVERTEX][0].getDiffN(copy_base);
228
229 if (nb_gauss_pts != copy_base_fun.size1())
230 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
231 "Wrong number of gauss pts");
232
233 if (nb_gauss_pts != copy_diff_base_fun.size1())
234 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
235 "Wrong number of gauss pts");
236
237 auto &base_fun = data.dataOnEntities[MBHEX][0].getN(base);
238 auto &diff_base_fun = data.dataOnEntities[MBHEX][0].getDiffN(base);
239 const int vol_order = data.dataOnEntities[MBHEX][0].getOrder();
240
241 const int nb_dofs = NBVOLUMEHEX_L2(vol_order);
242 base_fun.resize(nb_gauss_pts, nb_dofs, false);
243 diff_base_fun.resize(nb_gauss_pts, 3 * nb_dofs, false);
244
245 if (nb_dofs) {
246 const std::array<int, 3> p{vol_order, vol_order, vol_order};
247 CHKERR DemkowiczHexAndQuad::L2_InteriorShapeFunctions_ONHEX(
248 p.data(), &*copy_base_fun.data().begin(),
249 &*copy_diff_base_fun.data().begin(), &*base_fun.data().begin(),
250 &*diff_base_fun.data().begin(), nb_gauss_pts);
251 }
252
253 MoFEMFunctionReturn(0);
254}
255
256MoFEMErrorCode HexPolynomialBase::getValueHdiv(MatrixDouble &pts) {
257 MoFEMFunctionBegin;
258
259 switch (cTx->bAse) {
260 case DEMKOWICZ_JACOBI_BASE:
261 CHKERR getValueHdivDemkowiczBase(pts);
262 break;
263 default:
264 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
265 }
266
267 MoFEMFunctionReturn(0);
268}
269
270MoFEMErrorCode HexPolynomialBase::getValueHdivDemkowiczBase(MatrixDouble &pts) {
271 MoFEMFunctionBegin;
272
273 auto &data = cTx->dAta;
274 const auto base = cTx->bAse;
275 const auto copy_base = cTx->copyNodeBase;
276 const int nb_gauss_pts = pts.size2();
277
278 auto &copy_base_fun = data.dataOnEntities[MBVERTEX][0].getN(copy_base);
279 auto &copy_diff_base_fun =
280 data.dataOnEntities[MBVERTEX][0].getDiffN(copy_base);
281
282 if (nb_gauss_pts != copy_base_fun.size1())
283 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
284 "Wrong number of gauss pts");
285
286 if (nb_gauss_pts != copy_diff_base_fun.size1())
287 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
288 "Wrong number of gauss pts");
289
290 // Quad
291 if (data.spacesOnEntities[MBQUAD].test(HDIV)) {
292
293 if (data.dataOnEntities[MBQUAD].size() != 6)
294 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
295 "Expected six data structures on Hex");
296
297 std::array<int, 6> order;
298 std::array<double *, 6> hdiv_face_n;
299 std::array<double *, 6> diff_hdiv_face_n;
300
301 bool sum_nb_dofs = false;
302 for (int ff = 0; ff != 6; ff++) {
303 if (data.dataOnEntities[MBQUAD][ff].getSense() == 0)
304 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
305 "Sense pn quad <%d> not set", ff);
306
307 order[ff] = data.dataOnEntities[MBQUAD][ff].getOrder();
308 if (data.facesNodes.size1() != 6)
309 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
310 "Expected six faces");
311 if (data.facesNodes.size2() != 4)
312 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
313 "Expected four nodes on face");
314
315 const int nb_dofs = NBFACEQUAD_DEMKOWICZ_HDIV(order[ff]);
316 auto &face_n = data.dataOnEntities[MBQUAD][ff].getN(base);
317 auto &diff_face_n = data.dataOnEntities[MBQUAD][ff].getDiffN(base);
318 face_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
319 diff_face_n.resize(nb_gauss_pts, 9 * nb_dofs, false);
320
321 hdiv_face_n[ff] = &*face_n.data().begin();
322 diff_hdiv_face_n[ff] = &*diff_face_n.data().begin();
323
324 sum_nb_dofs |= (nb_dofs > 0);
325 }
326
327 if (sum_nb_dofs) {
328 DemkowiczHexAndQuad::Hdiv_FaceShapeFunctions_ONHEX(
329 &*data.facesNodes.data().begin(),
330 &*data.facesNodesOrder.data().begin(), order.data(),
331 &*copy_base_fun.data().begin(), &*copy_diff_base_fun.data().begin(),
332 hdiv_face_n.data(), diff_hdiv_face_n.data(), nb_gauss_pts);
333
334 } else {
335 for (int ff = 0; ff != 6; ff++) {
336 data.dataOnEntities[MBQUAD][ff].getN(base).resize(nb_gauss_pts, 0,
337 false);
338 data.dataOnEntities[MBQUAD][ff].getDiffN(base).resize(nb_gauss_pts, 0,
339 false);
340 }
341 }
342
343 } else {
344 for (int ff = 0; ff != 6; ff++) {
345 data.dataOnEntities[MBQUAD][ff].getN(base).resize(nb_gauss_pts, 0, false);
346 data.dataOnEntities[MBQUAD][ff].getDiffN(base).resize(nb_gauss_pts, 0,
347 false);
348 }
349 }
350
351 // Hex
352 if (data.spacesOnEntities[MBHEX].test(HDIV)) {
353
354 const int order = data.dataOnEntities[MBHEX][0].getOrder();
355 const int nb_dofs_family =
356 NBVOLUMEHEX_DEMKOWICZ_FAMILY_HDIV(order, order, order);
357
358 volFamily.resize(3, 3 * nb_dofs_family * nb_gauss_pts);
359 diffVolFamily.resize(3, 9 * nb_dofs_family * nb_gauss_pts);
360 if (nb_dofs_family) {
361
362 std::array<double *, 3> family_ptr = {&volFamily(0, 0), &volFamily(1, 0),
363 &volFamily(2, 0)};
364 std::array<double *, 3> diff_family_ptr = {
365 &diffVolFamily(0, 0), &diffVolFamily(1, 0), &diffVolFamily(2, 0)};
366
367 std::array<int, 3> p{order, order, order};
368 DemkowiczHexAndQuad::Hdiv_InteriorShapeFunctions_ONHEX(
369 p.data(), &*copy_base_fun.data().begin(),
370 &*copy_diff_base_fun.data().begin(), family_ptr.data(),
371 diff_family_ptr.data(), nb_gauss_pts);
372
373 const int nb_vol_dofs = NBVOLUMEHEX_DEMKOWICZ_HDIV(order);
374 auto &face_n = data.dataOnEntities[MBHEX][0].getN(base);
375 auto &diff_face_n = data.dataOnEntities[MBHEX][0].getDiffN(base);
376 face_n.resize(nb_gauss_pts, 3 * nb_vol_dofs, false);
377 diff_face_n.resize(nb_gauss_pts, 9 * nb_vol_dofs, false);
378
379 auto ptr_f0 = &(volFamily(0, 0));
380 auto ptr_f1 = &(volFamily(1, 0));
381 auto ptr_f2 = &(volFamily(2, 0));
382 double *ptr = &face_n(0, 0);
383 for (int n = 0; n != volFamily.size2() / 3; ++n) {
384 for (int j = 0; j != 3; ++j) {
385 *ptr = *ptr_f0;
386 ++ptr;
387 ++ptr_f0;
388 }
389 for (int j = 0; j != 3; ++j) {
390 *ptr = *ptr_f1;
391 ++ptr;
392 ++ptr_f1;
393 }
394 for (int j = 0; j != 3; ++j) {
395 *ptr = *ptr_f2;
396 ++ptr;
397 ++ptr_f2;
398 }
399 }
400
401 auto diff_ptr_f0 = &(diffVolFamily(0, 0));
402 auto diff_ptr_f1 = &(diffVolFamily(1, 0));
403 auto diff_ptr_f2 = &(diffVolFamily(2, 0));
404 double *diff_ptr = &diff_face_n(0, 0);
405 for (int n = 0; n != diffVolFamily.size2() / 9; ++n) {
406 for (int j = 0; j != 9; ++j) {
407 *diff_ptr = *diff_ptr_f0;
408 ++diff_ptr;
409 ++diff_ptr_f0;
410 }
411 for (int j = 0; j != 9; ++j) {
412 *diff_ptr = *diff_ptr_f1;
413 ++diff_ptr;
414 ++diff_ptr_f1;
415 }
416 for (int j = 0; j != 9; ++j) {
417 *diff_ptr = *diff_ptr_f2;
418 ++diff_ptr;
419 ++diff_ptr_f2;
420 }
421 }
422 } else {
423 data.dataOnEntities[MBHEX][0].getN(base).resize(nb_gauss_pts, 0, false);
424 data.dataOnEntities[MBHEX][0].getDiffN(base).resize(nb_gauss_pts, 0,
425 false);
426 }
427
428 } else {
429 data.dataOnEntities[MBHEX][0].getN(base).resize(nb_gauss_pts, 0, false);
430 data.dataOnEntities[MBHEX][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
431 }
432
433 MoFEMFunctionReturn(0);
434}
435
436MoFEMErrorCode HexPolynomialBase::getValueHcurl(MatrixDouble &pts) {
437 MoFEMFunctionBegin;
438
439 switch (cTx->bAse) {
440 case DEMKOWICZ_JACOBI_BASE:
441 CHKERR getValueHcurlDemkowiczBase(pts);
442 break;
443 default:
444 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
445 }
446
447 MoFEMFunctionReturn(0);
448}
449
450MoFEMErrorCode
451HexPolynomialBase::getValueHcurlDemkowiczBase(MatrixDouble &pts) {
452 MoFEMFunctionBegin;
453
454 auto &data = cTx->dAta;
455 const auto base = cTx->bAse;
456 const auto copy_base = cTx->copyNodeBase;
457 const int nb_gauss_pts = pts.size2();
458
459 auto &copy_base_fun = data.dataOnEntities[MBVERTEX][0].getN(copy_base);
460 auto &copy_diff_base_fun =
461 data.dataOnEntities[MBVERTEX][0].getDiffN(copy_base);
462
463 if (nb_gauss_pts != copy_base_fun.size1())
464 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
465 "Wrong number of gauss pts");
466
467 if (nb_gauss_pts != copy_diff_base_fun.size1())
468 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
469 "Wrong number of gauss pts");
470
471 // edges
472 if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
473 std::array<int, 12> sense;
474 std::array<int, 12> order;
475 if (data.dataOnEntities[MBEDGE].size() != 12)
476 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
477 "Expected 12 edges data structures on Hex");
478
479 std::array<double *, 12> hcurl_edge_n;
480 std::array<double *, 12> diff_hcurl_edge_n;
481 bool sum_nb_dofs = false;
482
483 for (int ee = 0; ee != 12; ++ee) {
484 if (data.dataOnEntities[MBEDGE][ee].getSense() == 0)
485 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
486 "Sense on edge <%d> on Hex not set", ee);
487
488 sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
489 order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
490 const int nb_dofs = NBEDGE_DEMKOWICZ_HCURL(order[ee]);
491 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
492 3 * nb_dofs, false);
493 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
494 9 * nb_dofs, false);
495 hcurl_edge_n[ee] =
496 &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
497 diff_hcurl_edge_n[ee] =
498 &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
499
500 sum_nb_dofs |= (nb_dofs > 0);
501 }
502
503 if (sum_nb_dofs) {
504 CHKERR DemkowiczHexAndQuad::Hcurl_EdgeShapeFunctions_ONHEX(
505 sense.data(), order.data(), &*copy_base_fun.data().begin(),
506 &*copy_diff_base_fun.data().begin(), hcurl_edge_n.data(),
507 diff_hcurl_edge_n.data(), nb_gauss_pts);
508 }
509
510 } else {
511 for (int ee = 0; ee != 12; ee++) {
512 data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
513 data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
514 false);
515 }
516 }
517
518 // Quad
519 if (data.spacesOnEntities[MBQUAD].test(HCURL)) {
520
521 if (data.dataOnEntities[MBQUAD].size() != 6)
522 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
523 "Expected six data structures on Hex");
524
525 std::array<int, 6> order;
526 double *face_family_ptr[6][2];
527 double *diff_face_family_ptr[6][2];
528
529 bool sum_nb_dofs = false;
530 for (int ff = 0; ff != 6; ff++) {
531 if (data.dataOnEntities[MBQUAD][ff].getSense() == 0)
532 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
533 "Sense pn quad <%d> not set", ff);
534
535 order[ff] = data.dataOnEntities[MBQUAD][ff].getOrder();
536 if (data.facesNodes.size1() != 6)
537 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
538 "Expected six faces");
539 if (data.facesNodes.size2() != 4)
540 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
541 "Expected four nodes on face");
542
543 const int nb_family_dofs =
544 NBFACEQUAD_DEMKOWICZ_FAMILY_HCURL(order[ff], order[ff]);
545 faceFamily[ff].resize(2, 3 * nb_family_dofs * nb_gauss_pts, false);
546 diffFaceFamily[ff].resize(2, 9 * nb_family_dofs * nb_gauss_pts, false);
547
548 if (nb_family_dofs) {
549 face_family_ptr[ff][0] = &(faceFamily[ff](0, 0));
550 face_family_ptr[ff][1] = &(faceFamily[ff](1, 0));
551 diff_face_family_ptr[ff][0] = &(diffFaceFamily[ff](0, 0));
552 diff_face_family_ptr[ff][1] = &(diffFaceFamily[ff](1, 0));
553 }
554
555 sum_nb_dofs |= (nb_family_dofs > 0);
556 }
557
558 if (sum_nb_dofs) {
559 DemkowiczHexAndQuad::Hcurl_FaceShapeFunctions_ONHEX(
560 &*data.facesNodes.data().begin(),
561 &*data.facesNodesOrder.data().begin(), order.data(),
562 &*copy_base_fun.data().begin(), &*copy_diff_base_fun.data().begin(),
563 face_family_ptr, diff_face_family_ptr, nb_gauss_pts);
564
565 for (int ff = 0; ff != 6; ++ff) {
566
567 // put family back
568
569 int nb_dofs = NBFACEQUAD_DEMKOWICZ_HCURL(order[ff]);
570 if (nb_dofs) {
571 auto &face_n = data.dataOnEntities[MBQUAD][ff].getN(base);
572 auto &diff_face_n = data.dataOnEntities[MBQUAD][ff].getDiffN(base);
573 face_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
574 diff_face_n.resize(nb_gauss_pts, 9 * nb_dofs, false);
575
576 auto ptr_f0 = &(faceFamily[ff](0, 0));
577 auto ptr_f1 = &(faceFamily[ff](1, 0));
578 double *ptr = &face_n(0, 0);
579 for (int n = 0; n != faceFamily[ff].size2() / 3; ++n) {
580 for (int j = 0; j != 3; ++j) {
581 *ptr = *ptr_f0;
582 ++ptr;
583 ++ptr_f0;
584 }
585 for (int j = 0; j != 3; ++j) {
586 *ptr = *ptr_f1;
587 ++ptr;
588 ++ptr_f1;
589 }
590 }
591
592 auto diff_ptr_f0 = &(diffFaceFamily[ff](0, 0));
593 auto diff_ptr_f1 = &(diffFaceFamily[ff](1, 0));
594 double *diff_ptr = &diff_face_n(0, 0);
595 for (int n = 0; n != diffFaceFamily[ff].size2() / 9; ++n) {
596 for (int j = 0; j != 9; ++j) {
597 *diff_ptr = *diff_ptr_f0;
598 ++diff_ptr;
599 ++diff_ptr_f0;
600 }
601 for (int j = 0; j != 9; ++j) {
602 *diff_ptr = *diff_ptr_f1;
603 ++diff_ptr;
604 ++diff_ptr_f1;
605 }
606 }
607 }
608 }
609 } else {
610 for (int ff = 0; ff != 6; ff++) {
611 data.dataOnEntities[MBQUAD][ff].getN(base).resize(nb_gauss_pts, 0,
612 false);
613 data.dataOnEntities[MBQUAD][ff].getDiffN(base).resize(nb_gauss_pts, 0,
614 false);
615 }
616 }
617
618 } else {
619 for (int ff = 0; ff != 6; ff++) {
620 data.dataOnEntities[MBQUAD][ff].getN(base).resize(nb_gauss_pts, 0, false);
621 data.dataOnEntities[MBQUAD][ff].getDiffN(base).resize(nb_gauss_pts, 0,
622 false);
623 }
624 }
625
626 // Hex
627 if (data.spacesOnEntities[MBHEX].test(HCURL)) {
628
629 const int order = data.dataOnEntities[MBHEX][0].getOrder();
630 const int nb_dofs = NBVOLUMEHEX_DEMKOWICZ_FAMILY_HCURL(order, order, order);
631
632 volFamily.resize(3, 3 * nb_dofs * nb_gauss_pts);
633 diffVolFamily.resize(3, 9 * nb_dofs * nb_gauss_pts);
634 if (nb_dofs) {
635
636 std::array<double *, 3> family_ptr = {&volFamily(0, 0), &volFamily(1, 0),
637 &volFamily(2, 0)};
638 std::array<double *, 3> diff_family_ptr = {
639 &diffVolFamily(0, 0), &diffVolFamily(1, 0), &diffVolFamily(2, 0)};
640
641 std::array<int, 3> p{order, order, order};
642 DemkowiczHexAndQuad::Hcurl_InteriorShapeFunctions_ONHEX(
643 p.data(), &*copy_base_fun.data().begin(),
644 &*copy_diff_base_fun.data().begin(), family_ptr.data(),
645 diff_family_ptr.data(), nb_gauss_pts);
646
647 const int nb_vol_dofs = NBVOLUMEHEX_DEMKOWICZ_HCURL(order);
648 auto &face_n = data.dataOnEntities[MBHEX][0].getN(base);
649 auto &diff_face_n = data.dataOnEntities[MBHEX][0].getDiffN(base);
650 face_n.resize(nb_gauss_pts, 3 * nb_vol_dofs, false);
651 diff_face_n.resize(nb_gauss_pts, 9 * nb_vol_dofs, false);
652
653 auto ptr_f0 = &(volFamily(0, 0));
654 auto ptr_f1 = &(volFamily(1, 0));
655 auto ptr_f2 = &(volFamily(2, 0));
656 double *ptr = &face_n(0, 0);
657 for (int n = 0; n != volFamily.size2() / 3; ++n) {
658 for (int j = 0; j != 3; ++j) {
659 *ptr = *ptr_f0;
660 ++ptr;
661 ++ptr_f0;
662 }
663 for (int j = 0; j != 3; ++j) {
664 *ptr = *ptr_f1;
665 ++ptr;
666 ++ptr_f1;
667 }
668 for (int j = 0; j != 3; ++j) {
669 *ptr = *ptr_f2;
670 ++ptr;
671 ++ptr_f2;
672 }
673 }
674
675 auto diff_ptr_f0 = &(diffVolFamily(0, 0));
676 auto diff_ptr_f1 = &(diffVolFamily(1, 0));
677 auto diff_ptr_f2 = &(diffVolFamily(2, 0));
678 double *diff_ptr = &diff_face_n(0, 0);
679 for (int n = 0; n != diffVolFamily.size2() / 9; ++n) {
680 for (int j = 0; j != 9; ++j) {
681 *diff_ptr = *diff_ptr_f0;
682 ++diff_ptr;
683 ++diff_ptr_f0;
684 }
685 for (int j = 0; j != 9; ++j) {
686 *diff_ptr = *diff_ptr_f1;
687 ++diff_ptr;
688 ++diff_ptr_f1;
689 }
690 for (int j = 0; j != 9; ++j) {
691 *diff_ptr = *diff_ptr_f2;
692 ++diff_ptr;
693 ++diff_ptr_f2;
694 }
695 }
696 } else {
697 data.dataOnEntities[MBHEX][0].getN(base).resize(nb_gauss_pts, 0, false);
698 data.dataOnEntities[MBHEX][0].getDiffN(base).resize(nb_gauss_pts, 0,
699 false);
700 }
701
702 } else {
703 data.dataOnEntities[MBHEX][0].getN(base).resize(nb_gauss_pts, 0, false);
704 data.dataOnEntities[MBHEX][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
705 }
706
707 MoFEMFunctionReturn(0);
708}
709
710MoFEMErrorCode
711HexPolynomialBase::getValue(MatrixDouble &pts,
712 boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
713 MoFEMFunctionBegin;
714
715 cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
716
717 int nb_gauss_pts = pts.size2();
718 if (!nb_gauss_pts)
719 MoFEMFunctionReturnHot(0);
720
721 if (pts.size1() < 3)
722 SETERRQ(
723 PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
724 "Wrong dimension of pts, should be at least 3 rows with coordinates");
725
726 switch (cTx->sPace) {
727 case H1:
728 CHKERR getValueH1(pts);
729 break;
730 case HDIV:
731 CHKERR getValueHdiv(pts);
732 break;
733 case HCURL:
734 CHKERR getValueHcurl(pts);
735 break;
736 case L2:
737 CHKERR getValueL2(pts);
738 break;
739 default:
740 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space");
741 }
742
743 MoFEMFunctionReturn(0);
744}
p
static Index< 'p', 3 > p
Definition: BasicFeTools.hpp:90
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
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
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
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
order
constexpr int order
Definition: dg_projection.cpp:18
n
FTensor::Index< 'n', SPACE_DIM > n
Definition: fluid_structure_eigenproblem.cpp:66
NBFACEQUAD_DEMKOWICZ_HDIV
#define NBFACEQUAD_DEMKOWICZ_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:145
NBFACEQUAD_DEMKOWICZ_FAMILY_HCURL
#define NBFACEQUAD_DEMKOWICZ_FAMILY_HCURL(P, Q)
Number of base functions on quad for Hcurl space.
Definition: h1_hdiv_hcurl_l2.h:116
NBEDGE_DEMKOWICZ_HCURL
#define NBEDGE_DEMKOWICZ_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:108
NBVOLUMEHEX_DEMKOWICZ_HCURL
#define NBVOLUMEHEX_DEMKOWICZ_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:124
NBFACEQUAD_H1
#define NBFACEQUAD_H1(P)
Number of base functions on quad for H1 space.
Definition: h1_hdiv_hcurl_l2.h:65
NBVOLUMEHEX_H1
#define NBVOLUMEHEX_H1(P)
Number of base functions on hex for H1 space.
Definition: h1_hdiv_hcurl_l2.h:93
NBVOLUMEHEX_DEMKOWICZ_HDIV
#define NBVOLUMEHEX_DEMKOWICZ_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:149
NBFACEQUAD_DEMKOWICZ_HCURL
#define NBFACEQUAD_DEMKOWICZ_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:118
NBEDGE_H1
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
Definition: h1_hdiv_hcurl_l2.h:55
NBVOLUMEHEX_L2
#define NBVOLUMEHEX_L2(P)
Number of base functions on hexahedron for L2 space.
Definition: h1_hdiv_hcurl_l2.h:37
NBVOLUMEHEX_DEMKOWICZ_FAMILY_HCURL
#define NBVOLUMEHEX_DEMKOWICZ_FAMILY_HCURL(P, Q, R)
Definition: h1_hdiv_hcurl_l2.h:121
NBVOLUMEHEX_DEMKOWICZ_FAMILY_HDIV
#define NBVOLUMEHEX_DEMKOWICZ_FAMILY_HDIV(P, Q, R)
Definition: h1_hdiv_hcurl_l2.h:147
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::DemkowiczHexAndQuad::Hcurl_EdgeShapeFunctions_ONHEX
MoFEMErrorCode Hcurl_EdgeShapeFunctions_ONHEX(int *sense, int *p, double *N, double *N_diff, double *edgeN[12], double *diff_edgeN[12], int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:1019
MoFEM::DemkowiczHexAndQuad::H1_EdgeShapeFunctions_ONHEX
MoFEMErrorCode H1_EdgeShapeFunctions_ONHEX(int *sense, int *p, double *N, double *N_diff, double *edgeN[12], double *diff_edgeN[12], int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:707
MoFEM::DemkowiczHexAndQuad::L2_InteriorShapeFunctions_ONHEX
MoFEMErrorCode L2_InteriorShapeFunctions_ONHEX(const int *p, double *N, double *N_diff, double *volN, double *diff_volN, int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:966
MoFEM::DemkowiczHexAndQuad::H1_FaceShapeFunctions_ONHEX
MoFEMErrorCode H1_FaceShapeFunctions_ONHEX(int *face_nodes, int *face_nodes_order, int *p, double *N, double *N_diff, double *faceN[6], double *diff_faceN[6], int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:795
MoFEM::DemkowiczHexAndQuad::Hdiv_FaceShapeFunctions_ONHEX
MoFEMErrorCode Hdiv_FaceShapeFunctions_ONHEX(int *face_nodes, int *face_nodes_order, int *p, double *N, double *diffN, double *faceN[6], double *div_faceN[6], int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:1404
MoFEM::DemkowiczHexAndQuad::Hcurl_FaceShapeFunctions_ONHEX
MoFEMErrorCode Hcurl_FaceShapeFunctions_ONHEX(int *face_nodes, int *face_nodes_order, int *p, double *N, double *N_diff, double *faceN[6][2], double *diff_faceN[6][2], int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:1139
MoFEM::DemkowiczHexAndQuad::H1_InteriorShapeFunctions_ONHEX
MoFEMErrorCode H1_InteriorShapeFunctions_ONHEX(const int *p, double *N, double *N_diff, double *faceN, double *diff_faceN, int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:905
MoFEM::DemkowiczHexAndQuad::Hdiv_InteriorShapeFunctions_ONHEX
MoFEMErrorCode Hdiv_InteriorShapeFunctions_ONHEX(int *p, double *N, double *N_diff, double *bubleN[3], double *div_bubleN[3], int nb_integration_pts)
Definition: EdgeQuadHexPolynomials.cpp:1544
MoFEM::DemkowiczHexAndQuad::Hcurl_InteriorShapeFunctions_ONHEX
MoFEMErrorCode Hcurl_InteriorShapeFunctions_ONHEX(int *p, double *N, double *N_diff, double *volN[3], double *diff_volN[3], int nb_integration_pts)
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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::bAse
const FieldApproximationBase bAse
Definition: EntPolynomialBaseCtx.hpp:38
MoFEM::HexPolynomialBase
Calculate base functions on tetrahedral.
Definition: HexPolynomialBase.hpp:19
MoFEM::HexPolynomialBase::getValueH1DemkowiczBase
MoFEMErrorCode getValueH1DemkowiczBase(MatrixDouble &pts)
Definition: HexPolynomialBase.cpp:34
MoFEM::HexPolynomialBase::getValueL2
MoFEMErrorCode getValueL2(MatrixDouble &pts)
Get base functions for L2 space.
Definition: HexPolynomialBase.cpp:202
MoFEM::HexPolynomialBase::getValueH1
MoFEMErrorCode getValueH1(MatrixDouble &pts)
Definition: HexPolynomialBase.cpp:20
MoFEM::HexPolynomialBase::faceFamily
std::array< MatrixDouble, 6 > faceFamily
Definition: HexPolynomialBase.hpp:82
MoFEM::HexPolynomialBase::getValueHdivDemkowiczBase
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)
Definition: HexPolynomialBase.cpp:270
MoFEM::HexPolynomialBase::getValueHcurl
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
Get base functions for Hcurl space.
Definition: HexPolynomialBase.cpp:436
MoFEM::HexPolynomialBase::getValueL2DemkowiczBase
MoFEMErrorCode getValueL2DemkowiczBase(MatrixDouble &pts)
Definition: HexPolynomialBase.cpp:216
MoFEM::HexPolynomialBase::getValueHdiv
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
Get base functions for Hdiv space.
Definition: HexPolynomialBase.cpp:256
MoFEM::HexPolynomialBase::getValueHcurlDemkowiczBase
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)
Definition: HexPolynomialBase.cpp:451
MoFEM::HexPolynomialBase::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: HexPolynomialBase.cpp:711
MoFEM::HexPolynomialBase::diffFaceFamily
std::array< MatrixDouble, 6 > diffFaceFamily
Definition: HexPolynomialBase.hpp:83
MoFEM::HexPolynomialBase::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition: HexPolynomialBase.cpp:13
MoFEM::HexPolynomialBase::cTx
EntPolynomialBaseCtx * cTx
Definition: HexPolynomialBase.hpp:30
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