v0.14.0
HexPolynomialBase.cpp
Go to the documentation of this file.
1 /** \file HexPolynomialBase.cpp
2 \brief Implementation of hierarchical bases on tetrahedral
3 
4 A l2, h1, h-div and h-curl spaces are implemented.
5 
6 */
7 
8 
9 
10 using namespace MoFEM;
11 
12 MoFEMErrorCode
13 HexPolynomialBase::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 
20 MoFEMErrorCode 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 
34 MoFEMErrorCode 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 
202 MoFEMErrorCode 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 
216 MoFEMErrorCode 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 
256 MoFEMErrorCode 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 
270 MoFEMErrorCode 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 
436 MoFEMErrorCode 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 
450 MoFEMErrorCode
451 HexPolynomialBase::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 
710 MoFEMErrorCode
711 HexPolynomialBase::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 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::HexPolynomialBase
Calculate base functions on tetrahedral.
Definition: HexPolynomialBase.hpp:19
H1
@ H1
continuous field
Definition: definitions.h:85
NBVOLUMEHEX_DEMKOWICZ_HCURL
#define NBVOLUMEHEX_DEMKOWICZ_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:124
NBEDGE_H1
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
Definition: h1_hdiv_hcurl_l2.h:55
MoFEM::HexPolynomialBase::getValueH1DemkowiczBase
MoFEMErrorCode getValueH1DemkowiczBase(MatrixDouble &pts)
Definition: HexPolynomialBase.cpp:34
MoFEM::HexPolynomialBase::getValueHdiv
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
Get base functions for Hdiv space.
Definition: HexPolynomialBase.cpp:256
MoFEM::EntPolynomialBaseCtx
Class used to pass element data to calculate base functions on tet,triangle,edge.
Definition: EntPolynomialBaseCtx.hpp:22
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::EntPolynomialBaseCtx::dAta
EntitiesFieldData & dAta
Definition: EntPolynomialBaseCtx.hpp:36
MoFEM::HexPolynomialBase::cTx
EntPolynomialBaseCtx * cTx
Definition: HexPolynomialBase.hpp:30
MoFEM::HexPolynomialBase::getValueL2DemkowiczBase
MoFEMErrorCode getValueL2DemkowiczBase(MatrixDouble &pts)
Definition: HexPolynomialBase.cpp:216
NBFACEQUAD_DEMKOWICZ_HCURL
#define NBFACEQUAD_DEMKOWICZ_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:118
order
constexpr int order
Definition: dg_projection.cpp:18
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::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_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::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::EntPolynomialBaseCtx::copyNodeBase
const FieldApproximationBase copyNodeBase
Definition: EntPolynomialBaseCtx.hpp:40
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::HexPolynomialBase::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: HexPolynomialBase.cpp:711
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::HexPolynomialBase::getValueHcurl
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
Get base functions for Hcurl space.
Definition: HexPolynomialBase.cpp:436
NBVOLUMEHEX_DEMKOWICZ_HDIV
#define NBVOLUMEHEX_DEMKOWICZ_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:149
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::EntPolynomialBaseCtx::sPace
const FieldSpace sPace
Definition: EntPolynomialBaseCtx.hpp:37
NBEDGE_DEMKOWICZ_HCURL
#define NBEDGE_DEMKOWICZ_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:108
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)
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
NBFACEQUAD_H1
#define NBFACEQUAD_H1(P)
Number of base functions on quad for H1 space.
Definition: h1_hdiv_hcurl_l2.h:65
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::EntPolynomialBaseCtx::bAse
const FieldApproximationBase bAse
Definition: EntPolynomialBaseCtx.hpp:38
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
convert.n
n
Definition: convert.py:82
MoFEM::HexPolynomialBase::volFamily
MatrixDouble volFamily
Definition: HexPolynomialBase.hpp:85
MoFEM::UnknownInterface
base class for all interface classes
Definition: UnknownInterface.hpp:34
MoFEM::HexPolynomialBase::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition: HexPolynomialBase.cpp:13
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
NBFACEQUAD_DEMKOWICZ_HDIV
#define NBFACEQUAD_DEMKOWICZ_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:145
MoFEM::HexPolynomialBase::getValueH1
MoFEMErrorCode getValueH1(MatrixDouble &pts)
Definition: HexPolynomialBase.cpp:20
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::HexPolynomialBase::getValueHdivDemkowiczBase
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)
Definition: HexPolynomialBase.cpp:270
MoFEM::HexPolynomialBase::diffFaceFamily
std::array< MatrixDouble, 6 > diffFaceFamily
Definition: HexPolynomialBase.hpp:83
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
NBVOLUMEHEX_DEMKOWICZ_FAMILY_HCURL
#define NBVOLUMEHEX_DEMKOWICZ_FAMILY_HCURL(P, Q, R)
Definition: h1_hdiv_hcurl_l2.h:121
NBVOLUMEHEX_H1
#define NBVOLUMEHEX_H1(P)
Number of base functions on hex for H1 space.
Definition: h1_hdiv_hcurl_l2.h:93
MoFEM::HexPolynomialBase::getValueHcurlDemkowiczBase
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)
Definition: HexPolynomialBase.cpp:451
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
MoFEM::HexPolynomialBase::getValueL2
MoFEMErrorCode getValueL2(MatrixDouble &pts)
Get base functions for L2 space.
Definition: HexPolynomialBase.cpp:202
NBVOLUMEHEX_DEMKOWICZ_FAMILY_HDIV
#define NBVOLUMEHEX_DEMKOWICZ_FAMILY_HDIV(P, Q, R)
Definition: h1_hdiv_hcurl_l2.h:147
MoFEM::HexPolynomialBase::faceFamily
std::array< MatrixDouble, 6 > faceFamily
Definition: HexPolynomialBase.hpp:82
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
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
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
NBVOLUMEHEX_L2
#define NBVOLUMEHEX_L2(P)
Number of base functions on hexahedron for L2 space.
Definition: h1_hdiv_hcurl_l2.h:37
MoFEM::HexPolynomialBase::diffVolFamily
MatrixDouble diffVolFamily
Definition: HexPolynomialBase.hpp:86
Generated by Doxygen 1.8.17 and hosted at University of Glasgow