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