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 using namespace MoFEM;
9 
10 MoFEMErrorCode
11 HexPolynomialBase::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 
18 MoFEMErrorCode 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 
32 MoFEMErrorCode 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 
200 MoFEMErrorCode 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 
214 MoFEMErrorCode 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 
254 MoFEMErrorCode 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 
268 MoFEMErrorCode 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  SETERRQ1(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 
434 MoFEMErrorCode 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 
448 MoFEMErrorCode
449 HexPolynomialBase::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  SETERRQ1(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  SETERRQ1(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 
708 MoFEMErrorCode
709 HexPolynomialBase::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 }
UBlasMatrix< double >
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference 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)
Number 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:32
MoFEM::HexPolynomialBase::getValueHdiv
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
Get base functions for Hdiv space.
Definition: HexPolynomialBase.cpp:254
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::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:214
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:41
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::HexPolynomialBase::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: HexPolynomialBase.cpp:709
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:434
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:39
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::EntPolynomialBaseCtx::spaceContinuity
const FieldContinuity spaceContinuity
Definition: EntPolynomialBaseCtx.hpp:38
MoFEM::HexPolynomialBase::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition: HexPolynomialBase.cpp:11
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:18
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::HexPolynomialBase::getValueHdivDemkowiczBase
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)
Definition: HexPolynomialBase.cpp:268
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:449
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MoFEM::HexPolynomialBase::getValueL2
MoFEMErrorCode getValueL2(MatrixDouble &pts)
Get base functions for L2 space.
Definition: HexPolynomialBase.cpp:200
NBVOLUMEHEX_DEMKOWICZ_FAMILY_HDIV
#define NBVOLUMEHEX_DEMKOWICZ_FAMILY_HDIV(P, Q, R)
Definition: h1_hdiv_hcurl_l2.h:147
CONTINUOUS
@ CONTINUOUS
Regular field.
Definition: definitions.h:100
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:429
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:359
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