v0.14.0
TetPolynomialBase.cpp
Go to the documentation of this file.
1 /** \file TetPolynomialBase.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 struct TetBaseCache {
11 
12  struct BaseCacheItem {
13  int order;
14  int nb_gauss_pts;
15  mutable MatrixDouble N;
16  mutable MatrixDouble diffN;
17  };
18 
19  struct HDivBaseCacheItem {
20 
21  int order;
22  int nb_gauss_pts;
23 
24  // Number of permeations for tetrahedron
25  // That is P(3, 4) = 24
26 
27  int n0;
28  int n1;
29  int n2;
30 
31  mutable MatrixDouble N;
32  mutable MatrixDouble diffN;
33  };
34 
35  using BaseCacheMI = boost::multi_index_container<
36  BaseCacheItem,
37  boost::multi_index::indexed_by<
38 
39  boost::multi_index::hashed_unique<
40 
41  composite_key<
42 
43  BaseCacheItem,
44  member<BaseCacheItem, int, &BaseCacheItem::order>,
45  member<BaseCacheItem, int, &BaseCacheItem::nb_gauss_pts>>>>
46 
47  >;
48 
49  using HDivBaseFaceCacheMI = boost::multi_index_container<
50  HDivBaseCacheItem,
51  boost::multi_index::indexed_by<
52 
53  boost::multi_index::hashed_unique<
54 
55  composite_key<
56 
57  HDivBaseCacheItem,
58  member<HDivBaseCacheItem, int, &HDivBaseCacheItem::order>,
59  member<HDivBaseCacheItem, int,
60  &HDivBaseCacheItem::nb_gauss_pts>,
61  member<HDivBaseCacheItem, int, &HDivBaseCacheItem::n0>,
62  member<HDivBaseCacheItem, int, &HDivBaseCacheItem::n1>,
63  member<HDivBaseCacheItem, int, &HDivBaseCacheItem::n2>>>>
64 
65  >;
66 
67  static std::map<const void *, BaseCacheMI> hdivBaseInteriorDemkowicz;
68  static std::map<const void *, BaseCacheMI> hdivBrokenBaseInteriorDemkowicz;
69  static std::map<const void *, HDivBaseFaceCacheMI> hDivBaseFaceDemkowicz;
70 
71  static std::map<const void *, BaseCacheMI> hdivBrokenBaseInteriorAinsworth;
72 };
73 
74 std::map<const void *, TetBaseCache::BaseCacheMI>
75  TetBaseCache::hdivBaseInteriorDemkowicz;
76 std::map<const void *, TetBaseCache::BaseCacheMI>
77  TetBaseCache::hdivBrokenBaseInteriorDemkowicz;
78 std::map<const void *, TetBaseCache::HDivBaseFaceCacheMI>
79  TetBaseCache::hDivBaseFaceDemkowicz;
80 
81 std::map<const void *, TetBaseCache::BaseCacheMI>
82  TetBaseCache::hdivBrokenBaseInteriorAinsworth;
83 
84 MoFEMErrorCode
85 TetPolynomialBase::query_interface(boost::typeindex::type_index type_index,
86  UnknownInterface **iface) const {
87 
88  MoFEMFunctionBeginHot;
89  *iface = const_cast<TetPolynomialBase *>(this);
90  MoFEMFunctionReturnHot(0);
91 }
92 
93 TetPolynomialBase::TetPolynomialBase(const void *ptr) : vPtr(ptr) {}
94 
95 TetPolynomialBase::~TetPolynomialBase() {
96  if (vPtr) {
97  if (TetBaseCache::hdivBaseInteriorDemkowicz.find(vPtr) !=
98  TetBaseCache::hdivBaseInteriorDemkowicz.end())
99  TetBaseCache::hdivBaseInteriorDemkowicz.erase(vPtr);
100  if (TetBaseCache::hDivBaseFaceDemkowicz.find(vPtr) !=
101  TetBaseCache::hDivBaseFaceDemkowicz.end())
102  TetBaseCache::hDivBaseFaceDemkowicz.erase(vPtr);
103  }
104 }
105 
106 MoFEMErrorCode TetPolynomialBase::getValueH1(MatrixDouble &pts) {
107  MoFEMFunctionBegin;
108 
109  switch (cTx->bAse) {
110  case AINSWORTH_LEGENDRE_BASE:
111  case AINSWORTH_LOBATTO_BASE:
112  CHKERR getValueH1AinsworthBase(pts);
113  break;
114  case AINSWORTH_BERNSTEIN_BEZIER_BASE:
115  CHKERR getValueH1BernsteinBezierBase(pts);
116  break;
117  default:
118  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
119  }
120 
121  MoFEMFunctionReturn(0);
122 }
123 
124 MoFEMErrorCode TetPolynomialBase::getValueH1AinsworthBase(MatrixDouble &pts) {
125  MoFEMFunctionBegin;
126 
127  EntitiesFieldData &data = cTx->dAta;
128  const FieldApproximationBase base = cTx->bAse;
129  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
130  double *diffL, const int dim) =
131  cTx->basePolynomialsType0;
132 
133  int nb_gauss_pts = pts.size2();
134 
135  int sense[6], order[6];
136  if (data.spacesOnEntities[MBEDGE].test(H1)) {
137  // edges
138  if (data.dataOnEntities[MBEDGE].size() != 6) {
139  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
140  }
141  double *h1_edge_n[6], *diff_h1_egde_n[6];
142  for (int ee = 0; ee != 6; ++ee) {
143  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
144  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
145  "data inconsistency");
146  }
147  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
148  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
149  int nb_dofs = NBEDGE_H1(data.dataOnEntities[MBEDGE][ee].getOrder());
150  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
151  false);
152  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
153  3 * nb_dofs, false);
154  h1_edge_n[ee] =
155  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
156  diff_h1_egde_n[ee] =
157  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
158  }
159  CHKERR H1_EdgeShapeFunctions_MBTET(
160  sense, order,
161  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
162  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
163  h1_edge_n, diff_h1_egde_n, nb_gauss_pts, base_polynomials);
164  } else {
165  for (int ee = 0; ee != 6; ++ee) {
166  data.dataOnEntities[MBEDGE][ee].getN(base).resize(0, 0, false);
167  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(0, 0, false);
168  }
169  }
170 
171  if (data.spacesOnEntities[MBTRI].test(H1)) {
172  // faces
173  if (data.dataOnEntities[MBTRI].size() != 4) {
174  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
175  }
176  double *h1_face_n[4], *diff_h1_face_n[4];
177  for (int ff = 0; ff != 4; ++ff) {
178  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
179  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
180  "data inconsistency");
181  }
182  int nb_dofs = NBFACETRI_H1(data.dataOnEntities[MBTRI][ff].getOrder());
183  order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
184  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, nb_dofs,
185  false);
186  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
187  3 * nb_dofs, false);
188  h1_face_n[ff] =
189  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
190  diff_h1_face_n[ff] =
191  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
192  }
193  if (data.facesNodes.size1() != 4) {
194  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
195  }
196  if (data.facesNodes.size2() != 3) {
197  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
198  }
199  CHKERR H1_FaceShapeFunctions_MBTET(
200  &*data.facesNodes.data().begin(), order,
201  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
202  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
203  h1_face_n, diff_h1_face_n, nb_gauss_pts, base_polynomials);
204 
205  } else {
206  for (int ff = 0; ff != 4; ++ff) {
207  data.dataOnEntities[MBTRI][ff].getN(base).resize(0, false);
208  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(0, 0, false);
209  }
210  }
211 
212  if (data.spacesOnEntities[MBTET].test(H1)) {
213  // volume
214  int order = data.dataOnEntities[MBTET][0].getOrder();
215  int nb_vol_dofs = NBVOLUMETET_H1(order);
216  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, nb_vol_dofs,
217  false);
218  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
219  3 * nb_vol_dofs, false);
220  CHKERR H1_VolumeShapeFunctions_MBTET(
221  data.dataOnEntities[MBTET][0].getOrder(),
222  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
223  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
224  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
225  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
226  nb_gauss_pts, base_polynomials);
227  } else {
228  data.dataOnEntities[MBTET][0].getN(base).resize(0, 0, false);
229  data.dataOnEntities[MBTET][0].getDiffN(base).resize(0, 0, false);
230  }
231 
232  MoFEMFunctionReturn(0);
233 }
234 
235 MoFEMErrorCode
236 TetPolynomialBase::getValueH1BernsteinBezierBase(MatrixDouble &pts) {
237  MoFEMFunctionBegin;
238 
239  EntitiesFieldData &data = cTx->dAta;
240  const std::string field_name = cTx->fieldName;
241  const int nb_gauss_pts = pts.size2();
242 
243  if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
244  (unsigned int)nb_gauss_pts)
245  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
246  "Base functions or nodes has wrong number of integration points "
247  "for base %s",
248  ApproximationBaseNames[NOBASE]);
249  auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
250 
251  auto get_alpha = [field_name](auto &data) -> MatrixInt & {
252  auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
253  if (!ptr)
254  ptr.reset(new MatrixInt());
255  return *ptr;
256  };
257 
258  auto get_base = [field_name](auto &data) -> MatrixDouble & {
259  auto &ptr = data.getBBNSharedPtr(field_name);
260  if (!ptr)
261  ptr.reset(new MatrixDouble());
262  return *ptr;
263  };
264 
265  auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
266  auto &ptr = data.getBBDiffNSharedPtr(field_name);
267  if (!ptr)
268  ptr.reset(new MatrixDouble());
269  return *ptr;
270  };
271 
272  auto get_alpha_by_name_ptr =
273  [](auto &data,
274  const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
275  return data.getBBAlphaIndicesSharedPtr(field_name);
276  };
277 
278  auto get_base_by_name_ptr =
279  [](auto &data,
280  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
281  return data.getBBNSharedPtr(field_name);
282  };
283 
284  auto get_diff_base_by_name_ptr =
285  [](auto &data,
286  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
287  return data.getBBDiffNSharedPtr(field_name);
288  };
289 
290  auto get_alpha_by_order_ptr =
291  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
292  return data.getBBAlphaIndicesByOrderSharedPtr(o);
293  };
294 
295  auto get_base_by_order_ptr =
296  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
297  return data.getBBNByOrderSharedPtr(o);
298  };
299 
300  auto get_diff_base_by_order_ptr =
301  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
302  return data.getBBDiffNByOrderSharedPtr(o);
303  };
304 
305  auto &vert_ent_data = data.dataOnEntities[MBVERTEX][0];
306  auto &vertex_alpha = get_alpha(vert_ent_data);
307  vertex_alpha.resize(4, 4, false);
308  vertex_alpha.clear();
309  for (int n = 0; n != 4; ++n)
310  vertex_alpha(n, n) = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n];
311 
312  auto &vert_get_n = get_base(vert_ent_data);
313  auto &vert_get_diff_n = get_diff_base(vert_ent_data);
314  vert_get_n.resize(nb_gauss_pts, 4, false);
315  vert_get_diff_n.resize(nb_gauss_pts, 12, false);
316  CHKERR BernsteinBezier::baseFunctionsTet(
317  1, lambda.size1(), vertex_alpha.size1(), &vertex_alpha(0, 0),
318  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &vert_get_n(0, 0),
319  &vert_get_diff_n(0, 0));
320  for (int n = 0; n != 4; ++n) {
321  const double f = boost::math::factorial<double>(
322  data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n]);
323  for (int g = 0; g != nb_gauss_pts; ++g) {
324  vert_get_n(g, n) *= f;
325  for (int d = 0; d != 3; ++d)
326  vert_get_diff_n(g, 3 * n + d) *= f;
327  }
328  }
329 
330  // edges
331  if (data.spacesOnEntities[MBEDGE].test(H1)) {
332  if (data.dataOnEntities[MBEDGE].size() != 6)
333  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
334  "Wrong size of ent data");
335 
336  constexpr int edges_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
337  {0, 3}, {1, 3}, {2, 3}};
338  for (int ee = 0; ee != 6; ++ee) {
339  auto &ent_data = data.dataOnEntities[MBEDGE][ee];
340  const int sense = ent_data.getSense();
341  if (sense == 0)
342  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
343  "Sense of the edge unknown");
344  const int order = ent_data.getOrder();
345  const int nb_dofs = NBEDGE_H1(order);
346 
347  if (nb_dofs) {
348  if (get_alpha_by_order_ptr(ent_data, order)) {
349  get_alpha_by_name_ptr(ent_data, field_name) =
350  get_alpha_by_order_ptr(ent_data, order);
351  get_base_by_name_ptr(ent_data, field_name) =
352  get_base_by_order_ptr(ent_data, order);
353  get_diff_base_by_name_ptr(ent_data, field_name) =
354  get_diff_base_by_order_ptr(ent_data, order);
355  } else {
356  auto &get_n = get_base(ent_data);
357  auto &get_diff_n = get_diff_base(ent_data);
358  get_n.resize(nb_gauss_pts, nb_dofs, false);
359  get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
360 
361  auto &edge_alpha = get_alpha(data.dataOnEntities[MBEDGE][ee]);
362  edge_alpha.resize(nb_dofs, 4, false);
363  CHKERR BernsteinBezier::generateIndicesEdgeTet(ee, order,
364  &edge_alpha(0, 0));
365  if (sense == -1) {
366  for (int i = 0; i != edge_alpha.size1(); ++i) {
367  int a = edge_alpha(i, edges_nodes[ee][0]);
368  edge_alpha(i, edges_nodes[ee][0]) =
369  edge_alpha(i, edges_nodes[ee][1]);
370  edge_alpha(i, edges_nodes[ee][1]) = a;
371  }
372  }
373  CHKERR BernsteinBezier::baseFunctionsTet(
374  order, lambda.size1(), edge_alpha.size1(), &edge_alpha(0, 0),
375  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
376  &get_diff_n(0, 0));
377 
378  get_alpha_by_order_ptr(ent_data, order) =
379  get_alpha_by_name_ptr(ent_data, field_name);
380  get_base_by_order_ptr(ent_data, order) =
381  get_base_by_name_ptr(ent_data, field_name);
382  get_diff_base_by_order_ptr(ent_data, order) =
383  get_diff_base_by_name_ptr(ent_data, field_name);
384  }
385  }
386  }
387  } else {
388  for (int ee = 0; ee != 6; ++ee) {
389  auto &ent_data = data.dataOnEntities[MBEDGE][ee];
390  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
391  auto &get_n = get_base(ent_data);
392  auto &get_diff_n = get_diff_base(ent_data);
393  get_n.resize(nb_gauss_pts, 0, false);
394  get_diff_n.resize(nb_gauss_pts, 0, false);
395  }
396  }
397 
398  // face
399  if (data.spacesOnEntities[MBTRI].test(H1)) {
400  if (data.dataOnEntities[MBTRI].size() != 4)
401  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
402  "Wrong size of ent data");
403  if (data.facesNodes.size1() != 4)
404  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
405  if (data.facesNodes.size2() != 3)
406  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
407 
408  for (int ff = 0; ff != 4; ++ff) {
409  auto &ent_data = data.dataOnEntities[MBTRI][ff];
410  const int order = ent_data.getOrder();
411  const int nb_dofs = NBFACETRI_H1(order);
412 
413  if (nb_dofs) {
414  if (get_alpha_by_order_ptr(ent_data, order)) {
415  get_alpha_by_name_ptr(ent_data, field_name) =
416  get_alpha_by_order_ptr(ent_data, order);
417  get_base_by_name_ptr(ent_data, field_name) =
418  get_base_by_order_ptr(ent_data, order);
419  get_diff_base_by_name_ptr(ent_data, field_name) =
420  get_diff_base_by_order_ptr(ent_data, order);
421  } else {
422 
423  auto &get_n = get_base(ent_data);
424  auto &get_diff_n = get_diff_base(ent_data);
425  get_n.resize(nb_gauss_pts, nb_dofs, false);
426  get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
427 
428  auto &face_alpha = get_alpha(ent_data);
429  face_alpha.resize(nb_dofs, 4, false);
430 
431  CHKERR BernsteinBezier::generateIndicesTriTet(ff, order,
432  &face_alpha(0, 0));
433  senseFaceAlpha.resize(face_alpha.size1(), face_alpha.size2(), false);
434  senseFaceAlpha.clear();
435  constexpr int tri_nodes[4][3] = {
436  {0, 1, 3}, {1, 2, 3}, {0, 2, 3}, {0, 1, 2}};
437  for (int d = 0; d != nb_dofs; ++d)
438  for (int n = 0; n != 3; ++n)
439  senseFaceAlpha(d, data.facesNodes(ff, n)) =
440  face_alpha(d, tri_nodes[ff][n]);
441  face_alpha.swap(senseFaceAlpha);
442  CHKERR BernsteinBezier::baseFunctionsTet(
443  order, lambda.size1(), face_alpha.size1(), &face_alpha(0, 0),
444  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
445  &get_diff_n(0, 0));
446 
447  get_alpha_by_order_ptr(ent_data, order) =
448  get_alpha_by_name_ptr(ent_data, field_name);
449  get_base_by_order_ptr(ent_data, order) =
450  get_base_by_name_ptr(ent_data, field_name);
451  get_diff_base_by_order_ptr(ent_data, order) =
452  get_diff_base_by_name_ptr(ent_data, field_name);
453  }
454  }
455  }
456  } else {
457  for (int ff = 0; ff != 4; ++ff) {
458  auto &ent_data = data.dataOnEntities[MBTRI][ff];
459  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
460  auto &get_n = get_base(ent_data);
461  auto &get_diff_n = get_diff_base(ent_data);
462  get_n.resize(nb_gauss_pts, 0, false);
463  get_diff_n.resize(nb_gauss_pts, 0, false);
464  }
465  }
466 
467  if (data.spacesOnEntities[MBTET].test(H1)) {
468  if (data.dataOnEntities[MBTET].size() != 1)
469  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
470  "Wrong size ent of ent data");
471 
472  auto &ent_data = data.dataOnEntities[MBTET][0];
473  const int order = ent_data.getOrder();
474  const int nb_dofs = NBVOLUMETET_H1(order);
475  if (get_alpha_by_order_ptr(ent_data, order)) {
476  get_alpha_by_name_ptr(ent_data, field_name) =
477  get_alpha_by_order_ptr(ent_data, order);
478  get_base_by_name_ptr(ent_data, field_name) =
479  get_base_by_order_ptr(ent_data, order);
480  get_diff_base_by_name_ptr(ent_data, field_name) =
481  get_diff_base_by_order_ptr(ent_data, order);
482  } else {
483 
484  auto &get_n = get_base(ent_data);
485  auto &get_diff_n = get_diff_base(ent_data);
486  get_n.resize(nb_gauss_pts, nb_dofs, false);
487  get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
488  if (nb_dofs) {
489  auto &tet_alpha = get_alpha(ent_data);
490  tet_alpha.resize(nb_dofs, 4, false);
491 
492  CHKERR BernsteinBezier::generateIndicesTetTet(order, &tet_alpha(0, 0));
493  CHKERR BernsteinBezier::baseFunctionsTet(
494  order, lambda.size1(), tet_alpha.size1(), &tet_alpha(0, 0),
495  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
496  &get_diff_n(0, 0));
497 
498  get_alpha_by_order_ptr(ent_data, order) =
499  get_alpha_by_name_ptr(ent_data, field_name);
500  get_base_by_order_ptr(ent_data, order) =
501  get_base_by_name_ptr(ent_data, field_name);
502  get_diff_base_by_order_ptr(ent_data, order) =
503  get_diff_base_by_name_ptr(ent_data, field_name);
504  }
505  }
506  } else {
507  auto &ent_data = data.dataOnEntities[MBTET][0];
508  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
509  auto &get_n = get_base(ent_data);
510  auto &get_diff_n = get_diff_base(ent_data);
511  get_n.resize(nb_gauss_pts, 0, false);
512  get_diff_n.resize(nb_gauss_pts, 0, false);
513  }
514 
515  MoFEMFunctionReturn(0);
516 }
517 
518 MoFEMErrorCode TetPolynomialBase::getValueL2(MatrixDouble &pts) {
519  MoFEMFunctionBegin;
520 
521  switch (cTx->bAse) {
522  case AINSWORTH_LEGENDRE_BASE:
523  case AINSWORTH_LOBATTO_BASE:
524  case DEMKOWICZ_JACOBI_BASE:
525  CHKERR getValueL2AinsworthBase(pts);
526  break;
527  case AINSWORTH_BERNSTEIN_BEZIER_BASE:
528  CHKERR getValueL2BernsteinBezierBase(pts);
529  break;
530  default:
531  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
532  }
533 
534  MoFEMFunctionReturn(0);
535 }
536 
537 MoFEMErrorCode TetPolynomialBase::getValueL2AinsworthBase(MatrixDouble &pts) {
538  MoFEMFunctionBegin;
539 
540  EntitiesFieldData &data = cTx->dAta;
541  const FieldApproximationBase base = cTx->bAse;
542  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
543  double *diffL, const int dim) =
544  cTx->basePolynomialsType0;
545 
546  int nb_gauss_pts = pts.size2();
547 
548  data.dataOnEntities[MBTET][0].getN(base).resize(
549  nb_gauss_pts, NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getOrder()),
550  false);
551  data.dataOnEntities[MBTET][0].getDiffN(base).resize(
552  nb_gauss_pts,
553  3 * NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getOrder()), false);
554 
555  CHKERR L2_Ainsworth_ShapeFunctions_MBTET(
556  data.dataOnEntities[MBTET][0].getOrder(),
557  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
558  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
559  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
560  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
561  nb_gauss_pts, base_polynomials);
562 
563  MoFEMFunctionReturn(0);
564 }
565 
566 MoFEMErrorCode
567 TetPolynomialBase::getValueL2BernsteinBezierBase(MatrixDouble &pts) {
568  MoFEMFunctionBegin;
569 
570  EntitiesFieldData &data = cTx->dAta;
571  const std::string field_name = cTx->fieldName;
572  const int nb_gauss_pts = pts.size2();
573 
574  if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
575  (unsigned int)nb_gauss_pts)
576  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
577  "Base functions or nodes has wrong number of integration points "
578  "for base %s",
579  ApproximationBaseNames[NOBASE]);
580  auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
581 
582  auto get_alpha = [field_name](auto &data) -> MatrixInt & {
583  auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
584  if (!ptr)
585  ptr.reset(new MatrixInt());
586  return *ptr;
587  };
588 
589  auto get_base = [field_name](auto &data) -> MatrixDouble & {
590  auto &ptr = data.getBBNSharedPtr(field_name);
591  if (!ptr)
592  ptr.reset(new MatrixDouble());
593  return *ptr;
594  };
595 
596  auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
597  auto &ptr = data.getBBDiffNSharedPtr(field_name);
598  if (!ptr)
599  ptr.reset(new MatrixDouble());
600  return *ptr;
601  };
602 
603  auto get_alpha_by_name_ptr =
604  [](auto &data,
605  const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
606  return data.getBBAlphaIndicesSharedPtr(field_name);
607  };
608 
609  auto get_base_by_name_ptr =
610  [](auto &data,
611  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
612  return data.getBBNSharedPtr(field_name);
613  };
614 
615  auto get_diff_base_by_name_ptr =
616  [](auto &data,
617  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
618  return data.getBBDiffNSharedPtr(field_name);
619  };
620 
621  auto get_alpha_by_order_ptr =
622  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
623  return data.getBBAlphaIndicesByOrderSharedPtr(o);
624  };
625 
626  auto get_base_by_order_ptr =
627  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
628  return data.getBBNByOrderSharedPtr(o);
629  };
630 
631  auto get_diff_base_by_order_ptr =
632  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
633  return data.getBBDiffNByOrderSharedPtr(o);
634  };
635 
636  if (data.spacesOnEntities[MBTET].test(L2)) {
637  if (data.dataOnEntities[MBTET].size() != 1)
638  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
639  "Wrong size ent of ent data");
640 
641  auto &ent_data = data.dataOnEntities[MBTET][0];
642  const int order = ent_data.getOrder();
643  const int nb_dofs = NBVOLUMETET_L2(order);
644 
645  if (get_alpha_by_order_ptr(ent_data, order)) {
646  get_alpha_by_name_ptr(ent_data, field_name) =
647  get_alpha_by_order_ptr(ent_data, order);
648  get_base_by_name_ptr(ent_data, field_name) =
649  get_base_by_order_ptr(ent_data, order);
650  get_diff_base_by_name_ptr(ent_data, field_name) =
651  get_diff_base_by_order_ptr(ent_data, order);
652  } else {
653 
654  auto &get_n = get_base(ent_data);
655  auto &get_diff_n = get_diff_base(ent_data);
656  get_n.resize(nb_gauss_pts, nb_dofs, false);
657  get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
658 
659  if (nb_dofs) {
660 
661  if (order == 0) {
662 
663  if (nb_dofs != 1)
664  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
665  "Inconsistent number of DOFs");
666 
667  auto &tri_alpha = get_alpha(ent_data);
668  tri_alpha.clear();
669  get_n(0, 0) = 1;
670  get_diff_n.clear();
671 
672  } else {
673 
674  if (nb_dofs != 4 + 6 * NBEDGE_H1(order) + 4 * NBFACETRI_H1(order) +
675  NBVOLUMETET_H1(order) ||
676  nb_dofs != NBVOLUMETET_L2(order))
677  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
678  "Inconsistent number of DOFs");
679 
680  auto &tet_alpha = get_alpha(ent_data);
681  tet_alpha.resize(nb_dofs, 4, false);
682 
683  CHKERR BernsteinBezier::generateIndicesVertexTet(order,
684  &tet_alpha(0, 0));
685  if (order > 1) {
686  std::array<int, 6> edge_n{order, order, order, order, order, order};
687  std::array<int *, 6> tet_edge_ptr{
688  &tet_alpha(4, 0),
689  &tet_alpha(4 + 1 * NBEDGE_H1(order), 0),
690  &tet_alpha(4 + 2 * NBEDGE_H1(order), 0),
691  &tet_alpha(4 + 3 * NBEDGE_H1(order), 0),
692  &tet_alpha(4 + 4 * NBEDGE_H1(order), 0),
693  &tet_alpha(4 + 5 * NBEDGE_H1(order), 0)};
694  CHKERR BernsteinBezier::generateIndicesEdgeTet(edge_n.data(),
695  tet_edge_ptr.data());
696  if (order > 2) {
697  std::array<int, 6> face_n{order, order, order, order};
698  std::array<int *, 6> tet_face_ptr{
699  &tet_alpha(4 + 6 * NBEDGE_H1(order), 0),
700  &tet_alpha(4 + 6 * NBEDGE_H1(order) + 1 * NBFACETRI_H1(order),
701  0),
702  &tet_alpha(4 + 6 * NBEDGE_H1(order) + 2 * NBFACETRI_H1(order),
703  0),
704  &tet_alpha(4 + 6 * NBEDGE_H1(order) + 3 * NBFACETRI_H1(order),
705  0),
706  };
707  CHKERR BernsteinBezier::generateIndicesTriTet(
708  face_n.data(), tet_face_ptr.data());
709  if (order > 3)
710  CHKERR BernsteinBezier::generateIndicesTetTet(
711  order,
712  &tet_alpha(
713  4 + 6 * NBEDGE_H1(order) + 4 * NBFACETRI_H1(order), 0));
714  }
715  }
716 
717  CHKERR BernsteinBezier::baseFunctionsTet(
718  order, lambda.size1(), tet_alpha.size1(), &tet_alpha(0, 0),
719  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
720  &get_diff_n(0, 0));
721 
722  get_alpha_by_order_ptr(ent_data, order) =
723  get_alpha_by_name_ptr(ent_data, field_name);
724  get_base_by_order_ptr(ent_data, order) =
725  get_base_by_name_ptr(ent_data, field_name);
726  get_diff_base_by_order_ptr(ent_data, order) =
727  get_diff_base_by_name_ptr(ent_data, field_name);
728  }
729  }
730  }
731  } else {
732  auto &ent_data = data.dataOnEntities[MBTET][0];
733  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
734  auto &get_n = get_base(ent_data);
735  auto &get_diff_n = get_diff_base(ent_data);
736  get_n.resize(nb_gauss_pts, 0, false);
737  get_diff_n.resize(nb_gauss_pts, 0, false);
738  }
739 
740  MoFEMFunctionReturn(0);
741 }
742 
743 MoFEMErrorCode TetPolynomialBase::getValueHdivAinsworthBaseImpl(
744 
745  MatrixDouble &pts,
746 
747  MatrixDouble &shape_functions, MatrixDouble &diff_shape_functions,
748 
749  int volume_order, std::array<int, 4> &faces_order,
750  std::array<int, 3 * 4> &faces_nodes
751 
752 ) {
753  MoFEMFunctionBegin;
754 
755  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
756  double *diffL, const int dim) =
757  cTx->basePolynomialsType0;
758 
759  int nb_gauss_pts = pts.size2();
760 
761  // face shape functions
762 
763  double *phi_f_e[4][3];
764  double *phi_f[4];
765  double *diff_phi_f_e[4][3];
766  double *diff_phi_f[4];
767 
768  N_face_edge.resize(4, 3, false);
769  N_face_bubble.resize(4, false);
770  diffN_face_edge.resize(4, 3, false);
771  diffN_face_bubble.resize(4, false);
772 
773  for (int ff = 0; ff != 4; ++ff) {
774  // three edges on face
775  for (int ee = 0; ee < 3; ee++) {
776  N_face_edge(ff, ee).resize(
777  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(faces_order[ff]),
778  false);
779  diffN_face_edge(ff, ee).resize(
780  nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_EDGE_HDIV(faces_order[ff]),
781  false);
782  phi_f_e[ff][ee] = &*N_face_edge(ff, ee).data().begin();
783  diff_phi_f_e[ff][ee] = &*diffN_face_edge(ff, ee).data().begin();
784  }
785  N_face_bubble[ff].resize(nb_gauss_pts,
786  3 * NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]),
787  false);
788  diffN_face_bubble[ff].resize(
789  nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]),
790  false);
791  phi_f[ff] = &*(N_face_bubble[ff].data().begin());
792  diff_phi_f[ff] = &*(diffN_face_bubble[ff].data().begin());
793  }
794 
795  CHKERR Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET(
796  faces_nodes.data(), faces_order.data(), &*shape_functions.data().begin(),
797  &*diff_shape_functions.data().begin(), phi_f_e, diff_phi_f_e,
798  nb_gauss_pts, base_polynomials);
799 
800  CHKERR Hdiv_Ainsworth_FaceBubbleShapeFunctions(
801  faces_nodes.data(), faces_order.data(), &*shape_functions.data().begin(),
802  &*diff_shape_functions.data().begin(), phi_f, diff_phi_f, nb_gauss_pts,
803  base_polynomials);
804 
805  // volume shape functions
806 
807  double *phi_v_e[6];
808  double *phi_v_f[4];
809  double *phi_v;
810  double *diff_phi_v_e[6];
811  double *diff_phi_v_f[4];
812  double *diff_phi_v;
813 
814  N_volume_edge.resize(6, false);
815  diffN_volume_edge.resize(6, false);
816  for (int ee = 0; ee != 6; ++ee) {
817  N_volume_edge[ee].resize(
818  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_EDGE_HDIV(volume_order), false);
819  diffN_volume_edge[ee].resize(
820  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_EDGE_HDIV(volume_order), false);
821  phi_v_e[ee] = &*(N_volume_edge[ee].data().begin());
822  diff_phi_v_e[ee] = &*(diffN_volume_edge[ee].data().begin());
823  }
824  if (NBVOLUMETET_AINSWORTH_EDGE_HDIV(volume_order))
825  CHKERR Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET(
826  volume_order, &*shape_functions.data().begin(),
827  &*diff_shape_functions.data().begin(), phi_v_e, diff_phi_v_e,
828  nb_gauss_pts, base_polynomials);
829 
830  N_volume_face.resize(4, false);
831  diffN_volume_face.resize(4, false);
832  for (int ff = 0; ff != 4; ++ff) {
833  N_volume_face[ff].resize(
834  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order), false);
835  diffN_volume_face[ff].resize(
836  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order), false);
837  phi_v_f[ff] = &*(N_volume_face[ff].data().begin());
838  diff_phi_v_f[ff] = &*(diffN_volume_face[ff].data().begin());
839  }
840  if (NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order))
841  CHKERR Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET(
842  volume_order, &*shape_functions.data().begin(),
843  &*diff_shape_functions.data().begin(), phi_v_f, diff_phi_v_f,
844  nb_gauss_pts, base_polynomials);
845 
846  N_volume_bubble.resize(
847  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order), false);
848  diffN_volume_bubble.resize(
849  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order), false);
850  phi_v = &*(N_volume_bubble.data().begin());
851  diff_phi_v = &*(diffN_volume_bubble.data().begin());
852  if (NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order))
853  CHKERR Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET(
854  volume_order, &*shape_functions.data().begin(),
855  &*diff_shape_functions.data().begin(), phi_v, diff_phi_v, nb_gauss_pts,
856  base_polynomials);
857 
858  MoFEMFunctionReturn(0);
859 }
860 
861 MoFEMErrorCode TetPolynomialBase::getValueHdivAinsworthBase(MatrixDouble &pts) {
862  MoFEMFunctionBegin;
863 
864  std::array<int, 4> faces_order;
865  std::array<int, 4 * 3> faces_nodes;
866 
867  FieldApproximationBase base = cTx->bAse;
868  EntitiesFieldData &data = cTx->dAta;
869 
870  int volume_order = data.dataOnEntities[MBTET][0].getOrder();
871  std::copy(data.facesNodes.data().begin(), data.facesNodes.data().end(),
872  faces_nodes.begin());
873  for (int ff = 0; ff != 4; ff++) {
874  faces_order[ff] = cTx->dAta.dataOnEntities[MBTRI][ff].getOrder();
875  }
876 
877  CHKERR getValueHdivAinsworthBaseImpl(
878  pts, data.dataOnEntities[MBVERTEX][0].getN(base),
879  data.dataOnEntities[MBVERTEX][0].getDiffN(base), volume_order,
880  faces_order, faces_nodes);
881 
882  // Set shape functions into data structure Shape functions hast to be put
883  // in arrays in order which guarantee hierarchical series of degrees of
884  // freedom, i.e. in other words dofs form sub-entities has to be group
885  // by order.
886 
887  FTensor::Index<'i', 3> i;
888  FTensor::Index<'j', 3> j;
889 
890  int nb_gauss_pts = pts.size2();
891 
892  // faces
893  if (data.dataOnEntities[MBTRI].size() != 4) {
894  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
895  }
896 
897  // face-face
898  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_f_f[] = {
899  getFTensor1FromPtr<3>(&*(N_face_bubble[0].data().begin())),
900  getFTensor1FromPtr<3>(&*(N_face_bubble[1].data().begin())),
901  getFTensor1FromPtr<3>(&*(N_face_bubble[2].data().begin())),
902  getFTensor1FromPtr<3>(&*(N_face_bubble[3].data().begin()))};
903  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_f_f[] = {
904  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[0].data().begin())),
905  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[1].data().begin())),
906  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[2].data().begin())),
907  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[3].data().begin()))};
908  // face-edge
909  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_f_e[] = {
910  getFTensor1FromPtr<3>(&*(N_face_edge(0, 0).data().begin())),
911  getFTensor1FromPtr<3>(&*(N_face_edge(0, 1).data().begin())),
912  getFTensor1FromPtr<3>(&*(N_face_edge(0, 2).data().begin())),
913  getFTensor1FromPtr<3>(&*(N_face_edge(1, 0).data().begin())),
914  getFTensor1FromPtr<3>(&*(N_face_edge(1, 1).data().begin())),
915  getFTensor1FromPtr<3>(&*(N_face_edge(1, 2).data().begin())),
916  getFTensor1FromPtr<3>(&*(N_face_edge(2, 0).data().begin())),
917  getFTensor1FromPtr<3>(&*(N_face_edge(2, 1).data().begin())),
918  getFTensor1FromPtr<3>(&*(N_face_edge(2, 2).data().begin())),
919  getFTensor1FromPtr<3>(&*(N_face_edge(3, 0).data().begin())),
920  getFTensor1FromPtr<3>(&*(N_face_edge(3, 1).data().begin())),
921  getFTensor1FromPtr<3>(&*(N_face_edge(3, 2).data().begin()))};
922  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_f_e[] = {
923  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(0, 0).data().begin())),
924  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(0, 1).data().begin())),
925  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(0, 2).data().begin())),
926  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(1, 0).data().begin())),
927  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(1, 1).data().begin())),
928  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(1, 2).data().begin())),
929  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(2, 0).data().begin())),
930  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(2, 1).data().begin())),
931  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(2, 2).data().begin())),
932  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(3, 0).data().begin())),
933  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(3, 1).data().begin())),
934  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(3, 2).data().begin()))};
935 
936  for (int ff = 0; ff != 4; ff++) {
937  data.dataOnEntities[MBTRI][ff].getN(base).resize(
938  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_HDIV(faces_order[ff]), false);
939  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
940  nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_HDIV(faces_order[ff]), false);
941  if (NBFACETRI_AINSWORTH_HDIV(faces_order[ff]) == 0)
942  continue;
943  double *base_ptr =
944  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
945  double *diff_base_ptr =
946  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
947  auto t_base = getFTensor1FromPtr<3>(base_ptr);
948  auto t_diff_base = getFTensor2HVecFromPtr<3, 3>(diff_base_ptr);
949 
950  for (int gg = 0; gg != nb_gauss_pts; gg++) {
951  for (int oo = 0; oo != faces_order[ff]; oo++) {
952 
953  // face-edge
954  for (int dd = NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
955  dd != NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
956  for (int ee = 0; ee != 3; ++ee) {
957  t_base(i) = t_base_f_e[ff * 3 + ee](i);
958  ++t_base;
959  ++t_base_f_e[ff * 3 + ee];
960  }
961  for (int ee = 0; ee != 3; ++ee) {
962  t_diff_base(i, j) = t_diff_base_f_e[ff * 3 + ee](i, j);
963  ++t_diff_base;
964  ++t_diff_base_f_e[ff * 3 + ee];
965  }
966  }
967 
968  // face-face
969  for (int dd = NBFACETRI_AINSWORTH_FACE_HDIV(oo);
970  dd != NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
971  t_base(i) = t_base_f_f[ff](i);
972  ++t_base;
973  ++t_base_f_f[ff];
974  t_diff_base(i, j) = t_diff_base_f_f[ff](i, j);
975  ++t_diff_base;
976  ++t_diff_base_f_f[ff];
977  }
978  }
979  }
980  }
981 
982  // volume
983  data.dataOnEntities[MBTET][0].getN(base).resize(
984  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_HDIV(volume_order), false);
985  data.dataOnEntities[MBTET][0].getDiffN(base).resize(
986  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_HDIV(volume_order), false);
987  if (NBVOLUMETET_AINSWORTH_HDIV(volume_order) > 0) {
988  double *base_ptr =
989  &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
990  double *diff_base_ptr =
991  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
992  auto t_base = getFTensor1FromPtr<3>(base_ptr);
993  auto t_diff_base = getFTensor2HVecFromPtr<3, 3>(diff_base_ptr);
994 
995  // volume-edge
996  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_v_e[] = {
997  getFTensor1FromPtr<3>(&*N_volume_edge[0].data().begin()),
998  getFTensor1FromPtr<3>(&*N_volume_edge[1].data().begin()),
999  getFTensor1FromPtr<3>(&*N_volume_edge[2].data().begin()),
1000  getFTensor1FromPtr<3>(&*N_volume_edge[3].data().begin()),
1001  getFTensor1FromPtr<3>(&*N_volume_edge[4].data().begin()),
1002  getFTensor1FromPtr<3>(&*N_volume_edge[5].data().begin())};
1003  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_e[] = {
1004  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[0].data().begin()),
1005  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[1].data().begin()),
1006  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[2].data().begin()),
1007  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[3].data().begin()),
1008  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[4].data().begin()),
1009  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[5].data().begin())};
1010 
1011  // volume-faces
1012  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_v_f[] = {
1013  getFTensor1FromPtr<3>(&*N_volume_face[0].data().begin()),
1014  getFTensor1FromPtr<3>(&*N_volume_face[1].data().begin()),
1015  getFTensor1FromPtr<3>(&*N_volume_face[2].data().begin()),
1016  getFTensor1FromPtr<3>(&*N_volume_face[3].data().begin())};
1017  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_f[] = {
1018  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[0].data().begin()),
1019  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[1].data().begin()),
1020  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[2].data().begin()),
1021  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[3].data().begin())};
1022 
1023  // volume-bubble
1024  base_ptr = &*(N_volume_bubble.data().begin());
1025  diff_base_ptr = &*(diffN_volume_bubble.data().begin());
1026  auto t_base_v = getFTensor1FromPtr<3>(base_ptr);
1027  auto t_diff_base_v = getFTensor2HVecFromPtr<3, 3>(diff_base_ptr);
1028 
1029  for (int gg = 0; gg != nb_gauss_pts; gg++) {
1030  for (int oo = 0; oo < volume_order; oo++) {
1031 
1032  // volume-edge
1033  for (int dd = NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo);
1034  dd < NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
1035  for (int ee = 0; ee < 6; ee++) {
1036  t_base(i) = t_base_v_e[ee](i);
1037  ++t_base;
1038  ++t_base_v_e[ee];
1039  t_diff_base(i, j) = t_diff_base_v_e[ee](i, j);
1040  ++t_diff_base;
1041  ++t_diff_base_v_e[ee];
1042  }
1043  }
1044 
1045  // volume-face
1046  for (int dd = NBVOLUMETET_AINSWORTH_FACE_HDIV(oo);
1047  dd < NBVOLUMETET_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
1048  for (int ff = 0; ff < 4; ff++) {
1049  t_base(i) = t_base_v_f[ff](i);
1050  ++t_base;
1051  ++t_base_v_f[ff];
1052  t_diff_base(i, j) = t_diff_base_v_f[ff](i, j);
1053  ++t_diff_base;
1054  ++t_diff_base_v_f[ff];
1055  }
1056  }
1057 
1058  // volume-bubble
1059  for (int dd = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo);
1060  dd < NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo + 1); dd++) {
1061  t_base(i) = t_base_v(i);
1062  ++t_base;
1063  ++t_base_v;
1064  t_diff_base(i, j) = t_diff_base_v(i, j);
1065  ++t_diff_base;
1066  ++t_diff_base_v;
1067  }
1068  }
1069  }
1070  }
1071 
1072  MoFEMFunctionReturn(0);
1073 }
1074 
1075 MoFEMErrorCode
1076 TetPolynomialBase::getValueHdivAinsworthBrokenBase(MatrixDouble &pts) {
1077  MoFEMFunctionBegin;
1078 
1079  // Set shape functions into data structure Shape functions has to be put
1080  // in arrays in order which guarantee hierarchical series of degrees of
1081  // freedom, i.e. in other words dofs form sub-entities has to be group
1082  // by order.
1083 
1084  FieldApproximationBase base = cTx->bAse;
1085  EntitiesFieldData &data = cTx->dAta;
1086 
1087  FTensor::Index<'i', 3> i;
1088  FTensor::Index<'j', 3> j;
1089 
1090  int volume_order = data.dataOnEntities[MBTET][0].getOrder();
1091  int nb_gauss_pts = pts.size2();
1092  int nb_dofs_face = NBFACETRI_AINSWORTH_HDIV(volume_order);
1093  int nb_dofs_volume = NBVOLUMETET_AINSWORTH_HDIV(volume_order);
1094  int nb_dofs = 4 * nb_dofs_face + nb_dofs_volume;
1095  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
1096  false);
1097  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 9 * nb_dofs,
1098  false);
1099  if (nb_dofs == 0)
1100  MoFEMFunctionReturnHot(0);
1101 
1102 
1103  auto get_interior_cache = [this]() -> TetBaseCache::BaseCacheMI * {
1104  if (vPtr) {
1105  auto it = TetBaseCache::hdivBrokenBaseInteriorAinsworth.find(vPtr);
1106  if (it != TetBaseCache::hdivBrokenBaseInteriorAinsworth.end()) {
1107  return &it->second;
1108  }
1109  }
1110  return nullptr;
1111  };
1112 
1113  auto interior_cache_ptr = get_interior_cache();
1114 
1115  if (interior_cache_ptr) {
1116  auto it =
1117  interior_cache_ptr->find(boost::make_tuple(volume_order, nb_gauss_pts));
1118  if (it != interior_cache_ptr->end()) {
1119  noalias(data.dataOnEntities[MBTET][0].getN(base)) = it->N;
1120  noalias(data.dataOnEntities[MBTET][0].getDiffN(base)) = it->diffN;
1121  MoFEMFunctionBeginHot(0);
1122  }
1123  }
1124 
1125  std::array<int, 4 * 3> faces_nodes{
1126  0, 1, 3, // face 0
1127 
1128  1, 2, 3, // face 1
1129 
1130  0, 2, 3, // face 2
1131 
1132  0, 1, 2 // face 3
1133 
1134  };
1135  std::array<int, 4> faces_order{volume_order, volume_order, volume_order,
1136  volume_order};
1137  CHKERR getValueHdivAinsworthBaseImpl(
1138  pts, data.dataOnEntities[MBVERTEX][0].getN(base),
1139  data.dataOnEntities[MBVERTEX][0].getDiffN(base), volume_order,
1140  faces_order, faces_nodes);
1141 
1142  auto *base_ptr = &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
1143  auto *diff_base_ptr =
1144  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
1145  auto t_base = getFTensor1FromPtr<3>(base_ptr);
1146  auto t_diff_base = getFTensor2HVecFromPtr<3, 3>(diff_base_ptr);
1147 
1148  // face-edge
1149  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_f_e[] = {
1150  getFTensor1FromPtr<3>(&*(N_face_edge(0, 0).data().begin())),
1151  getFTensor1FromPtr<3>(&*(N_face_edge(0, 1).data().begin())),
1152  getFTensor1FromPtr<3>(&*(N_face_edge(0, 2).data().begin())),
1153  getFTensor1FromPtr<3>(&*(N_face_edge(1, 0).data().begin())),
1154  getFTensor1FromPtr<3>(&*(N_face_edge(1, 1).data().begin())),
1155  getFTensor1FromPtr<3>(&*(N_face_edge(1, 2).data().begin())),
1156  getFTensor1FromPtr<3>(&*(N_face_edge(2, 0).data().begin())),
1157  getFTensor1FromPtr<3>(&*(N_face_edge(2, 1).data().begin())),
1158  getFTensor1FromPtr<3>(&*(N_face_edge(2, 2).data().begin())),
1159  getFTensor1FromPtr<3>(&*(N_face_edge(3, 0).data().begin())),
1160  getFTensor1FromPtr<3>(&*(N_face_edge(3, 1).data().begin())),
1161  getFTensor1FromPtr<3>(&*(N_face_edge(3, 2).data().begin()))};
1162  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_f_e[] = {
1163  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(0, 0).data().begin())),
1164  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(0, 1).data().begin())),
1165  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(0, 2).data().begin())),
1166  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(1, 0).data().begin())),
1167  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(1, 1).data().begin())),
1168  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(1, 2).data().begin())),
1169  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(2, 0).data().begin())),
1170  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(2, 1).data().begin())),
1171  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(2, 2).data().begin())),
1172  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(3, 0).data().begin())),
1173  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(3, 1).data().begin())),
1174  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_edge(3, 2).data().begin()))};
1175 
1176  // face-face
1177  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_f_f[] = {
1178  getFTensor1FromPtr<3>(&*(N_face_bubble[0].data().begin())),
1179  getFTensor1FromPtr<3>(&*(N_face_bubble[1].data().begin())),
1180  getFTensor1FromPtr<3>(&*(N_face_bubble[2].data().begin())),
1181  getFTensor1FromPtr<3>(&*(N_face_bubble[3].data().begin()))};
1182  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_f_f[] = {
1183  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[0].data().begin())),
1184  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[1].data().begin())),
1185  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[2].data().begin())),
1186  getFTensor2HVecFromPtr<3, 3>(&*(diffN_face_bubble[3].data().begin()))};
1187 
1188  // volume-edge
1189  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_v_e[] = {
1190  getFTensor1FromPtr<3>(&*N_volume_edge[0].data().begin()),
1191  getFTensor1FromPtr<3>(&*N_volume_edge[1].data().begin()),
1192  getFTensor1FromPtr<3>(&*N_volume_edge[2].data().begin()),
1193  getFTensor1FromPtr<3>(&*N_volume_edge[3].data().begin()),
1194  getFTensor1FromPtr<3>(&*N_volume_edge[4].data().begin()),
1195  getFTensor1FromPtr<3>(&*N_volume_edge[5].data().begin())};
1196  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_e[] = {
1197  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[0].data().begin()),
1198  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[1].data().begin()),
1199  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[2].data().begin()),
1200  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[3].data().begin()),
1201  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[4].data().begin()),
1202  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_edge[5].data().begin())};
1203 
1204  // volume-faces
1205  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_v_f[] = {
1206  getFTensor1FromPtr<3>(&*N_volume_face[0].data().begin()),
1207  getFTensor1FromPtr<3>(&*N_volume_face[1].data().begin()),
1208  getFTensor1FromPtr<3>(&*N_volume_face[2].data().begin()),
1209  getFTensor1FromPtr<3>(&*N_volume_face[3].data().begin())};
1210  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_f[] = {
1211  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[0].data().begin()),
1212  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[1].data().begin()),
1213  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[2].data().begin()),
1214  getFTensor2HVecFromPtr<3, 3>(&*diffN_volume_face[3].data().begin())};
1215 
1216  // volume-bubble
1217  auto *base_vol_ptr = &*(N_volume_bubble.data().begin());
1218  auto *diff_base_vol_ptr = &*(diffN_volume_bubble.data().begin());
1219  auto t_base_v = getFTensor1FromPtr<3>(base_vol_ptr);
1220  auto t_diff_base_v = getFTensor2HVecFromPtr<3, 3>(diff_base_vol_ptr);
1221 
1222  // int count_dofs = 0;
1223 
1224  for (int gg = 0; gg != nb_gauss_pts; gg++) {
1225  for (int oo = 0; oo < volume_order; oo++) {
1226 
1227  // faces-edge (((P) > 0) ? (P) : 0)
1228  for (int dd = NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
1229  dd != NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
1230  for (auto ff = 0; ff != 4; ++ff) {
1231  for (int ee = 0; ee != 3; ++ee) {
1232  t_base(i) = t_base_f_e[ff * 3 + ee](i);
1233  ++t_base;
1234  ++t_base_f_e[ff * 3 + ee];
1235  // ++count_dofs;
1236  }
1237  for (int ee = 0; ee != 3; ++ee) {
1238  t_diff_base(i, j) = t_diff_base_f_e[ff * 3 + ee](i, j);
1239  ++t_diff_base;
1240  ++t_diff_base_f_e[ff * 3 + ee];
1241  }
1242  }
1243  }
1244 
1245  // face-face (P - 1) * (P - 2) / 2
1246  for (int dd = NBFACETRI_AINSWORTH_FACE_HDIV(oo);
1247  dd != NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
1248  for (auto ff = 0; ff != 4; ++ff) {
1249  t_base(i) = t_base_f_f[ff](i);
1250  ++t_base;
1251  ++t_base_f_f[ff];
1252  t_diff_base(i, j) = t_diff_base_f_f[ff](i, j);
1253  ++t_diff_base;
1254  ++t_diff_base_f_f[ff];
1255  // ++count_dofs;
1256  }
1257  }
1258 
1259  // volume-edge (P - 1)
1260  for (int dd = NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo);
1261  dd != NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
1262  for (int ee = 0; ee < 6; ee++) {
1263  t_base(i) = t_base_v_e[ee](i);
1264  ++t_base;
1265  ++t_base_v_e[ee];
1266  t_diff_base(i, j) = t_diff_base_v_e[ee](i, j);
1267  ++t_diff_base;
1268  ++t_diff_base_v_e[ee];
1269  // ++count_dofs;
1270  }
1271  }
1272  // volume-face (P - 1) * (P - 2)
1273  for (int dd = NBVOLUMETET_AINSWORTH_FACE_HDIV(oo);
1274  dd < NBVOLUMETET_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
1275  for (int ff = 0; ff < 4; ff++) {
1276  t_base(i) = t_base_v_f[ff](i);
1277  ++t_base;
1278  ++t_base_v_f[ff];
1279  t_diff_base(i, j) = t_diff_base_v_f[ff](i, j);
1280  ++t_diff_base;
1281  ++t_diff_base_v_f[ff];
1282  // ++count_dofs;
1283  }
1284  }
1285  // volume-bubble (P - 3) * (P - 2) * (P - 1) / 2
1286  for (int dd = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo);
1287  dd < NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo + 1); dd++) {
1288  t_base(i) = t_base_v(i);
1289  ++t_base;
1290  ++t_base_v;
1291  t_diff_base(i, j) = t_diff_base_v(i, j);
1292  ++t_diff_base;
1293  ++t_diff_base_v;
1294  // ++count_dofs;
1295  }
1296  }
1297  }
1298 
1299 // #ifdef NDEBUG
1300 // if (nb_dofs != count_dofs / nb_gauss_pts)
1301 // SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1302 // "Number of dofs %d is different than expected %d", count_dofs,
1303 // nb_dofs);
1304 // #endif // NDEBUG
1305 
1306  if (interior_cache_ptr) {
1307  auto p = interior_cache_ptr->emplace(
1308  TetBaseCache::BaseCacheItem{volume_order, nb_gauss_pts});
1309  p.first->N = data.dataOnEntities[MBTET][0].getN(base);
1310  p.first->diffN = data.dataOnEntities[MBTET][0].getDiffN(base);
1311  }
1312 
1313  MoFEMFunctionReturn(0);
1314 }
1315 
1316 MoFEMErrorCode TetPolynomialBase::getValueHdivDemkowiczBase(MatrixDouble &pts) {
1317  MoFEMFunctionBegin;
1318 
1319  EntitiesFieldData &data = cTx->dAta;
1320  const FieldApproximationBase base = cTx->bAse;
1321  if (base != DEMKOWICZ_JACOBI_BASE) {
1322  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1323  "This should be used only with DEMKOWICZ_JACOBI_BASE "
1324  "but base is %s",
1325  ApproximationBaseNames[base]);
1326  }
1327  int nb_gauss_pts = pts.size2();
1328 
1329  int volume_order = data.dataOnEntities[MBTET][0].getOrder();
1330 
1331  int p_f[4];
1332  double *phi_f[4];
1333  double *diff_phi_f[4];
1334 
1335  auto get_face_cache_ptr = [this]() -> TetBaseCache::HDivBaseFaceCacheMI * {
1336  if (vPtr) {
1337  auto it = TetBaseCache::hDivBaseFaceDemkowicz.find(vPtr);
1338  if (it != TetBaseCache::hDivBaseFaceDemkowicz.end()) {
1339  return &it->second;
1340  }
1341  }
1342  return nullptr;
1343  };
1344 
1345  auto face_cache_ptr = get_face_cache_ptr();
1346 
1347  // Calculate base function on tet faces
1348  for (int ff = 0; ff != 4; ff++) {
1349  int face_order = data.dataOnEntities[MBTRI][ff].getOrder();
1350  int order = volume_order > face_order ? volume_order : face_order;
1351  data.dataOnEntities[MBTRI][ff].getN(base).resize(
1352  nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
1353  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
1354  nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
1355  if (NBFACETRI_DEMKOWICZ_HDIV(order) == 0)
1356  continue;
1357 
1358  if (face_cache_ptr) {
1359  auto it = face_cache_ptr->find(boost::make_tuple(
1360 
1361  face_order, nb_gauss_pts,
1362 
1363  data.facesNodes(ff, 0), data.facesNodes(ff, 1), data.facesNodes(ff, 2)
1364 
1365  ));
1366  if (it != face_cache_ptr->end()) {
1367  noalias(data.dataOnEntities[MBTRI][ff].getN(base)) = it->N;
1368  noalias(data.dataOnEntities[MBTRI][ff].getDiffN(base)) = it->diffN;
1369  continue;
1370  }
1371  }
1372 
1373  p_f[ff] = order;
1374  phi_f[ff] = &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1375  diff_phi_f[ff] =
1376  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1377 
1378  CHKERR Hdiv_Demkowicz_Face_MBTET_ON_FACE(
1379  &data.facesNodes(ff, 0), order,
1380  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1381  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1382  phi_f[ff], diff_phi_f[ff], nb_gauss_pts, 4);
1383  if (face_cache_ptr) {
1384  auto p = face_cache_ptr->emplace(TetBaseCache::HDivBaseCacheItem{
1385  face_order, nb_gauss_pts, data.facesNodes(ff, 0),
1386  data.facesNodes(ff, 1), data.facesNodes(ff, 2)});
1387  p.first->N = data.dataOnEntities[MBTRI][ff].getN(base);
1388  p.first->diffN = data.dataOnEntities[MBTRI][ff].getDiffN(base);
1389  }
1390  }
1391 
1392  auto get_interior_cache = [this]() -> TetBaseCache::BaseCacheMI * {
1393  if (vPtr) {
1394  auto it = TetBaseCache::hdivBaseInteriorDemkowicz.find(vPtr);
1395  if (it != TetBaseCache::hdivBaseInteriorDemkowicz.end()) {
1396  return &it->second;
1397  }
1398  }
1399  return nullptr;
1400  };
1401 
1402  auto interior_cache_ptr = get_interior_cache();
1403 
1404  // Calculate base functions in tet interior
1405  if (NBVOLUMETET_DEMKOWICZ_HDIV(volume_order) > 0) {
1406  data.dataOnEntities[MBTET][0].getN(base).resize(
1407  nb_gauss_pts, 3 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
1408  data.dataOnEntities[MBTET][0].getDiffN(base).resize(
1409  nb_gauss_pts, 9 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
1410 
1411  for (int v = 0; v != 1; ++v) {
1412  if (interior_cache_ptr) {
1413  auto it = interior_cache_ptr->find(
1414  boost::make_tuple(volume_order, nb_gauss_pts));
1415  if (it != interior_cache_ptr->end()) {
1416  noalias(data.dataOnEntities[MBTET][0].getN(base)) = it->N;
1417  noalias(data.dataOnEntities[MBTET][0].getDiffN(base)) = it->diffN;
1418  continue;
1419  }
1420  }
1421 
1422  double *phi_v = &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
1423  double *diff_phi_v =
1424  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
1425 
1426  CHKERR Hdiv_Demkowicz_Interior_MBTET(
1427  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
1428  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), p_f, phi_f,
1429  diff_phi_f, phi_v, diff_phi_v, nb_gauss_pts);
1430  if (interior_cache_ptr) {
1431  auto p = interior_cache_ptr->emplace(
1432  TetBaseCache::BaseCacheItem{volume_order, nb_gauss_pts});
1433  p.first->N = data.dataOnEntities[MBTET][0].getN(base);
1434  p.first->diffN = data.dataOnEntities[MBTET][0].getDiffN(base);
1435  }
1436  }
1437  }
1438 
1439  // Set size of face base correctly
1440  for (int ff = 0; ff != 4; ff++) {
1441  int face_order = data.dataOnEntities[MBTRI][ff].getOrder();
1442  data.dataOnEntities[MBTRI][ff].getN(base).resize(
1443  nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
1444  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
1445  nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
1446  }
1447 
1448  MoFEMFunctionReturn(0);
1449 }
1450 
1451 MoFEMErrorCode
1452 TetPolynomialBase::getValueHdivDemkowiczBrokenBase(MatrixDouble &pts) {
1453  MoFEMFunctionBegin;
1454 
1455  EntitiesFieldData &data = cTx->dAta;
1456  const FieldApproximationBase base = cTx->bAse;
1457  if (base != DEMKOWICZ_JACOBI_BASE) {
1458  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1459  "This should be used only with DEMKOWICZ_JACOBI_BASE "
1460  "but base is %s",
1461  ApproximationBaseNames[base]);
1462  }
1463  int nb_gauss_pts = pts.size2();
1464 
1465  int volume_order = data.dataOnEntities[MBTET][0].getOrder();
1466  int nb_dofs_face = NBFACETRI_DEMKOWICZ_HDIV(volume_order);
1467  int nb_dofs_volume = NBVOLUMETET_DEMKOWICZ_HDIV(volume_order);
1468  int nb_dofs = 4 * nb_dofs_face + nb_dofs_volume;
1469  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
1470  false);
1471  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 9 * nb_dofs,
1472  false);
1473  if (nb_dofs == 0)
1474  MoFEMFunctionReturnHot(0);
1475 
1476  auto get_interior_cache = [this]() -> TetBaseCache::BaseCacheMI * {
1477  if (vPtr) {
1478  auto it = TetBaseCache::hdivBrokenBaseInteriorDemkowicz.find(vPtr);
1479  if (it != TetBaseCache::hdivBrokenBaseInteriorDemkowicz.end()) {
1480  return &it->second;
1481  }
1482  }
1483  return nullptr;
1484  };
1485 
1486  auto interior_cache_ptr = get_interior_cache();
1487 
1488  if (interior_cache_ptr) {
1489  auto it =
1490  interior_cache_ptr->find(boost::make_tuple(volume_order, nb_gauss_pts));
1491  if (it != interior_cache_ptr->end()) {
1492  noalias(data.dataOnEntities[MBTET][0].getN(base)) = it->N;
1493  noalias(data.dataOnEntities[MBTET][0].getDiffN(base)) = it->diffN;
1494  MoFEMFunctionBeginHot(0);
1495  }
1496  }
1497 
1498  std::array<MatrixDouble, 4> face_base_fun{
1499  MatrixDouble(nb_gauss_pts, 3 * nb_dofs_face),
1500  MatrixDouble(nb_gauss_pts, 3 * nb_dofs_face),
1501  MatrixDouble(nb_gauss_pts, 3 * nb_dofs_face),
1502  MatrixDouble(nb_gauss_pts, 3 * nb_dofs_face)};
1503  std::array<MatrixDouble, 4> face_diff_base{
1504  MatrixDouble(nb_gauss_pts, 9 * nb_dofs_face),
1505  MatrixDouble(nb_gauss_pts, 9 * nb_dofs_face),
1506  MatrixDouble(nb_gauss_pts, 9 * nb_dofs_face),
1507  MatrixDouble(nb_gauss_pts, 9 * nb_dofs_face)};
1508 
1509  int face_node[4][3] = {{0, 1, 3}, {1, 2, 3}, {0, 2, 3}, {0, 1, 2}};
1510 
1511  std::array<int, 4> p_f{volume_order, volume_order, volume_order,
1512  volume_order};
1513  std::array<double *, 4> phi_f{
1514  &*face_base_fun[0].data().begin(), &*face_base_fun[1].data().begin(),
1515  &*face_base_fun[2].data().begin(), &*face_base_fun[3].data().begin()};
1516  std::array<double *, 4> diff_phi_f{
1517  &*face_diff_base[0].data().begin(), &*face_diff_base[1].data().begin(),
1518  &*face_diff_base[2].data().begin(), &*face_diff_base[3].data().begin()};
1519 
1520  // Calculate base function on tet faces
1521  for (int ff = 0; ff != 4; ff++) {
1522  CHKERR Hdiv_Demkowicz_Face_MBTET_ON_FACE(
1523  face_node[ff], p_f[ff],
1524  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1525  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1526  phi_f[ff], diff_phi_f[ff], nb_gauss_pts, 4);
1527  }
1528 
1529  MatrixDouble vol_bases(nb_gauss_pts, 3 * nb_dofs_volume);
1530  MatrixDouble vol_diff_bases(nb_gauss_pts, 9 * nb_dofs_volume);
1531  auto *phi_v = &*vol_bases.data().begin();
1532  auto *diff_phi_v = &*vol_diff_bases.data().begin();
1533  CHKERR Hdiv_Demkowicz_Interior_MBTET(
1534  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
1535  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), p_f.data(),
1536  phi_f.data(), diff_phi_f.data(), phi_v, diff_phi_v, nb_gauss_pts);
1537 
1538  // faces
1539  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_v_f[] = {
1540  getFTensor1FromPtr<3>(phi_f[0]), getFTensor1FromPtr<3>(phi_f[1]),
1541  getFTensor1FromPtr<3>(phi_f[2]), getFTensor1FromPtr<3>(phi_f[3])};
1542  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v_f[] = {
1543  getFTensor2HVecFromPtr<3, 3>(diff_phi_f[0]),
1544  getFTensor2HVecFromPtr<3, 3>(diff_phi_f[1]),
1545  getFTensor2HVecFromPtr<3, 3>(diff_phi_f[2]),
1546  getFTensor2HVecFromPtr<3, 3>(diff_phi_f[3])};
1547 
1548  // volumes
1549  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_base_v =
1550  getFTensor1FromPtr<3>(&*vol_bases.data().begin());
1551  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_diff_base_v =
1552  getFTensor2HVecFromPtr<3, 3>(&*vol_diff_bases.data().begin());
1553 
1554  auto t_base = getFTensor1FromPtr<3>(
1555  &*data.dataOnEntities[MBTET][0].getN(base).data().begin());
1556  auto t_diff_base = getFTensor2HVecFromPtr<3, 3>(
1557  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin());
1558 
1559  i_FTIndex<3> i;
1560  j_FTIndex<3> j;
1561 
1562  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
1563  for (int oo = 0; oo < volume_order; oo++) {
1564  // face
1565  for (auto dd = NBFACETRI_DEMKOWICZ_HDIV(oo);
1566  dd != NBFACETRI_DEMKOWICZ_HDIV(oo + 1); ++dd) {
1567  for (auto ff = 0; ff != 4; ++ff) {
1568  t_base(i) = t_base_v_f[ff](i);
1569  ++t_base;
1570  ++t_base_v_f[ff];
1571  t_diff_base(i, j) = t_diff_base_v_f[ff](i, j);
1572  ++t_diff_base;
1573  ++t_diff_base_v_f[ff];
1574  }
1575  }
1576  // volume
1577  for (auto dd = NBVOLUMETET_DEMKOWICZ_HDIV(oo);
1578  dd != NBVOLUMETET_DEMKOWICZ_HDIV(oo + 1); ++dd) {
1579  t_base(i) = t_base_v(i);
1580  ++t_base;
1581  ++t_base_v;
1582  t_diff_base(i, j) = t_diff_base_v(i, j);
1583  ++t_diff_base;
1584  ++t_diff_base_v;
1585  }
1586  }
1587  }
1588 
1589  if (interior_cache_ptr) {
1590  auto p = interior_cache_ptr->emplace(
1591  TetBaseCache::BaseCacheItem{volume_order, nb_gauss_pts});
1592  p.first->N = data.dataOnEntities[MBTET][0].getN(base);
1593  p.first->diffN = data.dataOnEntities[MBTET][0].getDiffN(base);
1594  }
1595 
1596  MoFEMFunctionReturn(0);
1597 }
1598 
1599 MoFEMErrorCode TetPolynomialBase::getValueHdiv(MatrixDouble &pts) {
1600  MoFEMFunctionBeginHot;
1601 
1602  switch (cTx->spaceContinuity) {
1603  case CONTINUOUS:
1604  switch (cTx->bAse) {
1605  case AINSWORTH_LEGENDRE_BASE:
1606  case AINSWORTH_LOBATTO_BASE:
1607  return getValueHdivAinsworthBase(pts);
1608  case DEMKOWICZ_JACOBI_BASE:
1609  return getValueHdivDemkowiczBase(pts);
1610  default:
1611  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1612  }
1613  break;
1614  case DISCONTINUOUS:
1615  switch (cTx->bAse) {
1616  case AINSWORTH_LEGENDRE_BASE:
1617  case AINSWORTH_LOBATTO_BASE:
1618  return getValueHdivAinsworthBrokenBase(pts);
1619  case DEMKOWICZ_JACOBI_BASE:
1620  return getValueHdivDemkowiczBrokenBase(pts);
1621  default:
1622  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1623  }
1624  break;
1625 
1626  default:
1627  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown continuity");
1628  }
1629 
1630  MoFEMFunctionReturnHot(0);
1631 }
1632 
1633 MoFEMErrorCode
1634 TetPolynomialBase::getValueHcurlAinsworthBase(MatrixDouble &pts) {
1635  MoFEMFunctionBegin;
1636 
1637  EntitiesFieldData &data = cTx->dAta;
1638  const FieldApproximationBase base = cTx->bAse;
1639  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
1640  double *diffL, const int dim) =
1641  cTx->basePolynomialsType0;
1642 
1643  int nb_gauss_pts = pts.size2();
1644 
1645  // edges
1646  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
1647  int sense[6], order[6];
1648  if (data.dataOnEntities[MBEDGE].size() != 6) {
1649  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1650  }
1651  double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
1652  for (int ee = 0; ee != 6; ee++) {
1653  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
1654  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1655  "data inconsistency");
1656  }
1657  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
1658  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
1659  int nb_dofs =
1660  NBEDGE_AINSWORTH_HCURL(data.dataOnEntities[MBEDGE][ee].getOrder());
1661  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
1662  3 * nb_dofs, false);
1663  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
1664  9 * nb_dofs, false);
1665  hcurl_edge_n[ee] =
1666  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
1667  diff_hcurl_edge_n[ee] =
1668  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
1669  }
1670  CHKERR Hcurl_Ainsworth_EdgeBaseFunctions_MBTET(
1671  sense, order,
1672  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1673  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1674  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
1675  } else {
1676  for (int ee = 0; ee != 6; ee++) {
1677  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
1678  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
1679  false);
1680  }
1681  }
1682 
1683  // triangles
1684  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
1685  int order[4];
1686  // faces
1687  if (data.dataOnEntities[MBTRI].size() != 4) {
1688  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1689  }
1690  double *hcurl_base_n[4], *diff_hcurl_base_n[4];
1691  for (int ff = 0; ff != 4; ff++) {
1692  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
1693  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1694  "data inconsistency");
1695  }
1696  order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
1697  int nb_dofs = NBFACETRI_AINSWORTH_HCURL(order[ff]);
1698  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1699  3 * nb_dofs, false);
1700  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1701  9 * nb_dofs, false);
1702  hcurl_base_n[ff] =
1703  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1704  diff_hcurl_base_n[ff] =
1705  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1706  }
1707  if (data.facesNodes.size1() != 4) {
1708  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1709  }
1710  if (data.facesNodes.size2() != 3) {
1711  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1712  }
1713  CHKERR Hcurl_Ainsworth_FaceFunctions_MBTET(
1714  &*data.facesNodes.data().begin(), order,
1715  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1716  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1717  hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts, base_polynomials);
1718  } else {
1719  for (int ff = 0; ff != 4; ff++) {
1720  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
1721  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
1722  false);
1723  }
1724  }
1725 
1726  if (data.spacesOnEntities[MBTET].test(HCURL)) {
1727 
1728  // volume
1729  int order = data.dataOnEntities[MBTET][0].getOrder();
1730  int nb_vol_dofs = NBVOLUMETET_AINSWORTH_HCURL(order);
1731  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
1732  3 * nb_vol_dofs, false);
1733  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
1734  9 * nb_vol_dofs, false);
1735  CHKERR Hcurl_Ainsworth_VolumeFunctions_MBTET(
1736  data.dataOnEntities[MBTET][0].getOrder(),
1737  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1738  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1739  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
1740  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
1741  nb_gauss_pts, base_polynomials);
1742 
1743  } else {
1744  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
1745  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
1746  }
1747 
1748  MoFEMFunctionReturn(0);
1749 }
1750 
1751 MoFEMErrorCode
1752 TetPolynomialBase::getValueHcurlDemkowiczBase(MatrixDouble &pts) {
1753  MoFEMFunctionBegin;
1754 
1755  EntitiesFieldData &data = cTx->dAta;
1756  const FieldApproximationBase base = cTx->bAse;
1757  if (base != DEMKOWICZ_JACOBI_BASE) {
1758  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1759  "This should be used only with DEMKOWICZ_JACOBI_BASE "
1760  "but base is %s",
1761  ApproximationBaseNames[base]);
1762  }
1763 
1764  int nb_gauss_pts = pts.size2();
1765 
1766  // edges
1767  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
1768  int sense[6], order[6];
1769  if (data.dataOnEntities[MBEDGE].size() != 6) {
1770  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1771  "wrong size of data structure, expected space for six edges "
1772  "but is %d",
1773  data.dataOnEntities[MBEDGE].size());
1774  }
1775  double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
1776  for (int ee = 0; ee != 6; ee++) {
1777  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
1778  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1779  "orintation of edges is not set");
1780  }
1781  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
1782  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
1783  int nb_dofs =
1784  NBEDGE_DEMKOWICZ_HCURL(data.dataOnEntities[MBEDGE][ee].getOrder());
1785  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
1786  3 * nb_dofs, false);
1787  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
1788  9 * nb_dofs, false);
1789  hcurl_edge_n[ee] =
1790  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
1791  diff_hcurl_edge_n[ee] =
1792  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
1793  }
1794  CHKERR Hcurl_Demkowicz_EdgeBaseFunctions_MBTET(
1795  sense, order,
1796  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1797  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1798  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts);
1799  } else {
1800  // No DOFs on edges, resize base function matrices, indicating that no
1801  // dofs on them.
1802  for (int ee = 0; ee != 6; ee++) {
1803  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
1804  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
1805  false);
1806  }
1807  }
1808 
1809  // triangles
1810  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
1811  int order[4];
1812  // faces
1813  if (data.dataOnEntities[MBTRI].size() != 4) {
1814  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1815  "data structure for storing face h-curl base have wrong size "
1816  "should be four but is %d",
1817  data.dataOnEntities[MBTRI].size());
1818  }
1819  double *hcurl_base_n[4], *diff_hcurl_base_n[4];
1820  for (int ff = 0; ff != 4; ff++) {
1821  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
1822  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1823  "orintation of face is not set");
1824  }
1825  order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
1826  int nb_dofs = NBFACETRI_DEMKOWICZ_HCURL(order[ff]);
1827  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1828  3 * nb_dofs, false);
1829  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1830  9 * nb_dofs, false);
1831  hcurl_base_n[ff] =
1832  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1833  diff_hcurl_base_n[ff] =
1834  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1835  }
1836  if (data.facesNodes.size1() != 4) {
1837  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1838  "data inconsistency, should be four faces");
1839  }
1840  if (data.facesNodes.size2() != 3) {
1841  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1842  "data inconsistency, should be three nodes on face");
1843  }
1844  CHKERR Hcurl_Demkowicz_FaceBaseFunctions_MBTET(
1845  &*data.facesNodes.data().begin(), order,
1846  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1847  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1848  hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts);
1849  } else {
1850  // No DOFs on faces, resize base function matrices, indicating that no
1851  // dofs on them.
1852  for (int ff = 0; ff != 4; ff++) {
1853  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
1854  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
1855  false);
1856  }
1857  }
1858 
1859  if (data.spacesOnEntities[MBTET].test(HCURL)) {
1860  // volume
1861  int order = data.dataOnEntities[MBTET][0].getOrder();
1862  int nb_vol_dofs = NBVOLUMETET_DEMKOWICZ_HCURL(order);
1863  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
1864  3 * nb_vol_dofs, false);
1865  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
1866  9 * nb_vol_dofs, false);
1867  CHKERR Hcurl_Demkowicz_VolumeBaseFunctions_MBTET(
1868  data.dataOnEntities[MBTET][0].getOrder(),
1869  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1870  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1871  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
1872  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
1873  nb_gauss_pts);
1874  } else {
1875  // No DOFs on faces, resize base function matrices, indicating that no
1876  // dofs on them.
1877  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
1878  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
1879  }
1880 
1881  MoFEMFunctionReturn(0);
1882 }
1883 
1884 MoFEMErrorCode TetPolynomialBase::getValueHcurl(MatrixDouble &pts) {
1885  MoFEMFunctionBegin;
1886 
1887  switch (cTx->bAse) {
1888  case AINSWORTH_LEGENDRE_BASE:
1889  case AINSWORTH_LOBATTO_BASE:
1890  CHKERR getValueHcurlAinsworthBase(pts);
1891  break;
1892  case DEMKOWICZ_JACOBI_BASE:
1893  CHKERR getValueHcurlDemkowiczBase(pts);
1894  break;
1895  default:
1896  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1897  }
1898 
1899  MoFEMFunctionReturn(0);
1900 }
1901 
1902 MoFEMErrorCode
1903 TetPolynomialBase::getValue(MatrixDouble &pts,
1904  boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
1905  MoFEMFunctionBegin;
1906 
1907  cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
1908 
1909  int nb_gauss_pts = pts.size2();
1910  if (!nb_gauss_pts)
1911  MoFEMFunctionReturnHot(0);
1912 
1913  if (pts.size1() < 3)
1914  SETERRQ(
1915  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1916  "Wrong dimension of pts, should be at least 3 rows with coordinates");
1917 
1918  if (cTx->bAse != AINSWORTH_BERNSTEIN_BEZIER_BASE) {
1919  const FieldApproximationBase base = cTx->bAse;
1920  EntitiesFieldData &data = cTx->dAta;
1921  if (cTx->copyNodeBase == LASTBASE) {
1922  data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4,
1923  false);
1924  CHKERR Tools::shapeFunMBTET(
1925  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1926  &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
1927  } else {
1928  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
1929  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
1930  }
1931  if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
1932  (unsigned int)nb_gauss_pts) {
1933  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1934  "Base functions or nodes has wrong number of integration points "
1935  "for base %s",
1936  ApproximationBaseNames[base]);
1937  }
1938  data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, false);
1939  std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
1940  data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
1941  }
1942 
1943  switch (cTx->spaceContinuity) {
1944  case CONTINUOUS:
1945 
1946  switch (cTx->sPace) {
1947  case H1:
1948  CHKERR getValueH1(pts);
1949  break;
1950  case HDIV:
1951  CHKERR getValueHdiv(pts);
1952  break;
1953  case HCURL:
1954  CHKERR getValueHcurl(pts);
1955  break;
1956  case L2:
1957  CHKERR getValueL2(pts);
1958  break;
1959  default:
1960  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space %s",
1961  FieldSpaceNames[cTx->sPace]);
1962  }
1963  break;
1964 
1965  case DISCONTINUOUS:
1966 
1967  switch (cTx->sPace) {
1968  case HDIV:
1969  CHKERR getValueHdiv(pts);
1970  break;
1971  default:
1972  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space %s",
1973  FieldSpaceNames[cTx->sPace]);
1974  }
1975  break;
1976 
1977  default:
1978  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown continuity");
1979  }
1980 
1981  MoFEMFunctionReturn(0);
1982 }
1983 
1984 template <typename T>
1985 auto tetCacheSwitch(const void *ptr, T &cache, std::string cache_name) {
1986  auto it = cache.find(ptr);
1987  if (it != cache.end()) {
1988  MOFEM_LOG_CHANNEL("WORLD");
1989  MOFEM_TAG_AND_LOG("WORLD", Sev::noisy, "TetPolynomialBase")
1990  << "Cache off " << cache_name << ": " << it->second.size();
1991  cache.erase(it);
1992  return false;
1993  } else {
1994  MOFEM_LOG_CHANNEL("WORLD");
1995  MOFEM_TAG_AND_LOG("WORLD", Sev::noisy, "TetPolynomialBase")
1996  << "Cache on " << cache_name;
1997  cache[ptr];
1998  return true;
1999  }
2000 }
2001 
2002 bool TetPolynomialBase::swichCacheHDivBaseFaceDemkowicz(const void *ptr) {
2003  return tetCacheSwitch(ptr, TetBaseCache::hDivBaseFaceDemkowicz,
2004  "hDivBaseFaceDemkowicz");
2005 }
2006 
2007 bool TetPolynomialBase::swichCacheHdivBaseInteriorDemkowicz(const void *ptr) {
2008  return tetCacheSwitch(ptr, TetBaseCache::hdivBaseInteriorDemkowicz,
2009  "hdivBaseInteriorDemkowicz");
2010 }
2011 
2012 bool TetPolynomialBase::swichCacheHdivBrokenBaseInteriorDemkowicz(
2013  const void *ptr) {
2014  return tetCacheSwitch(ptr, TetBaseCache::hdivBrokenBaseInteriorDemkowicz,
2015  "hdivBrokenBaseInteriorDemkowicz");
2016 }
2017 
2018 void TetPolynomialBase::swichCacheHDivBaseDemkowiczOn(std::vector<void *> v) {
2019  for (auto fe_ptr : v) {
2020  if (!TetPolynomialBase::swichCacheHDivBaseFaceDemkowicz(fe_ptr)) {
2021  TetPolynomialBase::swichCacheHDivBaseFaceDemkowicz(fe_ptr);
2022  }
2023  if (!TetPolynomialBase::swichCacheHdivBaseInteriorDemkowicz(fe_ptr)) {
2024  TetPolynomialBase::swichCacheHdivBaseInteriorDemkowicz(fe_ptr);
2025  }
2026  if (!TetPolynomialBase::swichCacheHdivBrokenBaseInteriorDemkowicz(fe_ptr)) {
2027  TetPolynomialBase::swichCacheHdivBrokenBaseInteriorDemkowicz(fe_ptr);
2028  }
2029  }
2030 }
2031 
2032 void TetPolynomialBase::swichCacheHDivBaseDemkowiczOff(std::vector<void *> v) {
2033  for (auto fe_ptr : v) {
2034  if (TetPolynomialBase::swichCacheHDivBaseFaceDemkowicz(fe_ptr)) {
2035  TetPolynomialBase::swichCacheHDivBaseFaceDemkowicz(fe_ptr);
2036  }
2037  if (TetPolynomialBase::swichCacheHdivBaseInteriorDemkowicz(fe_ptr)) {
2038  TetPolynomialBase::swichCacheHdivBaseInteriorDemkowicz(fe_ptr);
2039  }
2040  if (TetPolynomialBase::swichCacheHdivBrokenBaseInteriorDemkowicz(fe_ptr)) {
2041  TetPolynomialBase::swichCacheHdivBrokenBaseInteriorDemkowicz(fe_ptr);
2042  }
2043  }
2044 }
2045 
2046 bool TetPolynomialBase::swichCacheHdivBrokenBaseInteriorAinsworth(
2047  const void *ptr) {
2048  return tetCacheSwitch(ptr, TetBaseCache::hdivBrokenBaseInteriorAinsworth,
2049  "hdivBrokenBaseInteriorAinsworth");
2050 }
2051 
2052 void TetPolynomialBase::swichCacheHDivBaseAinsworthOn(std::vector<void *> v) {
2053  for (auto fe_ptr : v) {
2054  if (!TetPolynomialBase::swichCacheHdivBrokenBaseInteriorAinsworth(fe_ptr)) {
2055  TetPolynomialBase::swichCacheHdivBrokenBaseInteriorAinsworth(fe_ptr);
2056  }
2057  }
2058 }
2059 
2060 void TetPolynomialBase::swichCacheHDivBaseAinsworthOff(std::vector<void *> v) {
2061  for (auto fe_ptr : v) {
2062  if (TetPolynomialBase::swichCacheHdivBrokenBaseInteriorAinsworth(fe_ptr)) {
2063  TetPolynomialBase::swichCacheHdivBrokenBaseInteriorAinsworth(fe_ptr);
2064  }
2065  }
2066 }
2067 
2068 void TetPolynomialBase::swichCacheHDivBaseOn(std::vector<void *> v) {
2069  swichCacheHDivBaseDemkowiczOn(v);
2070  swichCacheHDivBaseAinsworthOn(v);
2071 }
2072 
2073 void TetPolynomialBase::swichCacheHDivBaseOff(std::vector<void *> v) {
2074  swichCacheHDivBaseDemkowiczOff(v);
2075  swichCacheHDivBaseAinsworthOff(v);
2076 }
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::TetPolynomialBase::swichCacheHDivBaseAinsworthOn
static void swichCacheHDivBaseAinsworthOn(std::vector< void * > v)
Definition: TetPolynomialBase.cpp:2052
MoFEM::Hcurl_Ainsworth_FaceFunctions_MBTET
MoFEMErrorCode Hcurl_Ainsworth_FaceFunctions_MBTET(int *face_nodes, int *p, double *N, double *diffN, double *phi_f[4], double *diff_phi_f[4], int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face H-curl functions.
Definition: Hcurl.cpp:1052
g
constexpr double g
Definition: shallow_wave.cpp:63
MoFEM::TetPolynomialBase::getValueH1
MoFEMErrorCode getValueH1(MatrixDouble &pts)
Get base functions for H1 space.
Definition: TetPolynomialBase.cpp:106
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET
MoFEMErrorCode Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_f[], double *diff_phi_v_f[], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: Hdiv.cpp:368
NBEDGE_H1
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
Definition: h1_hdiv_hcurl_l2.h:55
MoFEM::Hdiv_Demkowicz_Face_MBTET_ON_FACE
MoFEMErrorCode Hdiv_Demkowicz_Face_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb)
Definition: Hdiv.cpp:617
TetBaseCache::BaseCacheItem::nb_gauss_pts
int nb_gauss_pts
Definition: TetPolynomialBase.cpp:14
TetBaseCache::HDivBaseCacheItem::diffN
MatrixDouble diffN
Definition: TetPolynomialBase.cpp:32
LASTBASE
@ LASTBASE
Definition: definitions.h:69
TetBaseCache::HDivBaseCacheItem::n2
int n2
Definition: TetPolynomialBase.cpp:29
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
L2_Ainsworth_ShapeFunctions_MBTET
PetscErrorCode L2_Ainsworth_ShapeFunctions_MBTET(int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Get base functions on tetrahedron for L2 space.
Definition: l2.c:74
MoFEM::TetPolynomialBase::getValueL2
MoFEMErrorCode getValueL2(MatrixDouble &pts)
Get base functions for L2 space.
Definition: TetPolynomialBase.cpp:518
NBVOLUMETET_DEMKOWICZ_HCURL
#define NBVOLUMETET_DEMKOWICZ_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:110
TetBaseCache::HDivBaseCacheItem::N
MatrixDouble N
Definition: TetPolynomialBase.cpp:31
MoFEM::EntPolynomialBaseCtx
Class used to pass element data to calculate base functions on tet,triangle,edge.
Definition: EntPolynomialBaseCtx.hpp:22
MoFEM::TetPolynomialBase::N_face_bubble
ublas::vector< MatrixDouble > N_face_bubble
Definition: TetPolynomialBase.hpp:119
TetBaseCache::hDivBaseFaceDemkowicz
static std::map< const void *, HDivBaseFaceCacheMI > hDivBaseFaceDemkowicz
Definition: TetPolynomialBase.cpp:69
MoFEM::TetPolynomialBase::getValueHcurlAinsworthBase
MoFEMErrorCode getValueHcurlAinsworthBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:1634
MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET
MoFEMErrorCode Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f_e[4][3], double *diff_phi_f_e[4][3], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base functions, Edge-based face functions by Ainsworth .
Definition: Hdiv.cpp:13
NOBASE
@ NOBASE
Definition: definitions.h:59
TetBaseCache::BaseCacheMI
boost::multi_index_container< BaseCacheItem, boost::multi_index::indexed_by< boost::multi_index::hashed_unique< composite_key< BaseCacheItem, member< BaseCacheItem, int, &BaseCacheItem::order >, member< BaseCacheItem, int, &BaseCacheItem::nb_gauss_pts > >> > > BaseCacheMI
Definition: TetPolynomialBase.cpp:47
ApproximationBaseNames
const static char *const ApproximationBaseNames[]
Definition: definitions.h:72
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
NBFACETRI_AINSWORTH_EDGE_HDIV
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:130
MoFEM::Hcurl_Demkowicz_EdgeBaseFunctions_MBTET
MoFEMErrorCode Hcurl_Demkowicz_EdgeBaseFunctions_MBTET(int *sense, int *p, double *n, double *diff_n, double *phi[], double *diff_phi[], int nb_integration_pts)
Edge based H-curl base functions on tetrahedral.
Definition: Hcurl.cpp:2068
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
TetBaseCache
Definition: TetPolynomialBase.cpp:10
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
TetBaseCache::HDivBaseCacheItem::n0
int n0
Definition: TetPolynomialBase.cpp:27
MoFEM::TetPolynomialBase::senseFaceAlpha
MatrixInt senseFaceAlpha
Definition: TetPolynomialBase.hpp:116
MoFEM::EntPolynomialBaseCtx::dAta
EntitiesFieldData & dAta
Definition: EntPolynomialBaseCtx.hpp:36
NBVOLUMETET_AINSWORTH_FACE_HDIV
#define NBVOLUMETET_AINSWORTH_FACE_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:134
MoFEM::Tools::diffShapeFunMBTET
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition: Tools.hpp:271
MoFEM::TetPolynomialBase::N_face_edge
ublas::matrix< MatrixDouble > N_face_edge
Definition: TetPolynomialBase.hpp:118
NBVOLUMETET_H1
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
Definition: h1_hdiv_hcurl_l2.h:75
MoFEM::Hdiv_Demkowicz_Interior_MBTET
MoFEMErrorCode Hdiv_Demkowicz_Interior_MBTET(int p, double *N, double *diffN, int p_face[], double *phi_f[4], double *diff_phi_f[4], double *phi_v, double *diff_phi_v, int gdim)
Definition: Hdiv.cpp:762
MoFEM::TetPolynomialBase::getValueH1AinsworthBase
MoFEMErrorCode getValueH1AinsworthBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:124
MoFEM::TetPolynomialBase::diffN_volume_edge
ublas::vector< MatrixDouble > diffN_volume_edge
Definition: TetPolynomialBase.hpp:126
NBVOLUMETET_L2
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
Definition: h1_hdiv_hcurl_l2.h:27
MoFEM::TetPolynomialBase::getValueL2AinsworthBase
MoFEMErrorCode getValueL2AinsworthBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:537
MoFEM::TetPolynomialBase::TetPolynomialBase
TetPolynomialBase(const void *ptr=nullptr)
Definition: TetPolynomialBase.cpp:93
MoFEM::EntitiesFieldData::spacesOnEntities
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
Definition: EntitiesFieldData.hpp:49
MoFEM::EntitiesFieldData::facesNodes
MatrixInt facesNodes
nodes on finite element faces
Definition: EntitiesFieldData.hpp:45
TetBaseCache::BaseCacheItem::N
MatrixDouble N
Definition: TetPolynomialBase.cpp:15
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
H1_EdgeShapeFunctions_MBTET
PetscErrorCode H1_EdgeShapeFunctions_MBTET(int *sense, int *p, double *N, double *diffN, double *edgeN[], double *diff_edgeN[], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: h1.c:274
order
constexpr int order
Definition: dg_projection.cpp:18
TetBaseCache::hdivBaseInteriorDemkowicz
static std::map< const void *, BaseCacheMI > hdivBaseInteriorDemkowicz
Definition: TetPolynomialBase.cpp:67
MoFEM::EntPolynomialBaseCtx::copyNodeBase
const FieldApproximationBase copyNodeBase
Definition: EntPolynomialBaseCtx.hpp:41
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::TetPolynomialBase::getValueHcurl
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
Get base functions for Hcurl space.
Definition: TetPolynomialBase.cpp:1884
MoFEM::BernsteinBezier::generateIndicesEdgeTet
static MoFEMErrorCode generateIndicesEdgeTet(const int N[], int *alpha[])
Definition: BernsteinBezier.cpp:138
TetBaseCache::HDivBaseCacheItem
Definition: TetPolynomialBase.cpp:19
MoFEM::BernsteinBezier::generateIndicesTriTet
static MoFEMErrorCode generateIndicesTriTet(const int N[], int *alpha[])
Definition: BernsteinBezier.cpp:174
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::TetPolynomialBase::swichCacheHdivBaseInteriorDemkowicz
static bool swichCacheHdivBaseInteriorDemkowicz(const void *ptr)
Definition: TetPolynomialBase.cpp:2007
NBVOLUMETET_AINSWORTH_VOLUME_HDIV
#define NBVOLUMETET_AINSWORTH_VOLUME_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:135
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::TetPolynomialBase::diffN_volume_bubble
MatrixDouble diffN_volume_bubble
Definition: TetPolynomialBase.hpp:128
MoFEM::getFTensor2HVecFromPtr< 3, 3 >
FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > getFTensor2HVecFromPtr< 3, 3 >(double *ptr)
Definition: Templates.hpp:946
MoFEM::TetPolynomialBase::N_volume_edge
ublas::vector< MatrixDouble > N_volume_edge
Definition: TetPolynomialBase.hpp:120
NBFACETRI_DEMKOWICZ_HDIV
#define NBFACETRI_DEMKOWICZ_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:139
MoFEM::Types::MatrixInt
UBlasMatrix< int > MatrixInt
Definition: Types.hpp:76
MoFEM::TetPolynomialBase::vPtr
const void * vPtr
Definition: TetPolynomialBase.hpp:44
TetBaseCache::HDivBaseFaceCacheMI
boost::multi_index_container< HDivBaseCacheItem, boost::multi_index::indexed_by< boost::multi_index::hashed_unique< composite_key< HDivBaseCacheItem, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::order >, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::nb_gauss_pts >, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::n0 >, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::n1 >, member< HDivBaseCacheItem, int, &HDivBaseCacheItem::n2 > >> > > HDivBaseFaceCacheMI
Definition: TetPolynomialBase.cpp:65
MoFEM::EntPolynomialBaseCtx::sPace
const FieldSpace sPace
Definition: EntPolynomialBaseCtx.hpp:37
MoFEM::BernsteinBezier::generateIndicesVertexTet
static MoFEMErrorCode generateIndicesVertexTet(const int N, int *alpha)
Definition: BernsteinBezier.cpp:128
NBVOLUMETET_DEMKOWICZ_HDIV
#define NBVOLUMETET_DEMKOWICZ_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:140
MoFEM::Hcurl_Demkowicz_VolumeBaseFunctions_MBTET
MoFEMErrorCode Hcurl_Demkowicz_VolumeBaseFunctions_MBTET(int p, double *n, double *diff_n, double *phi, double *diff_phi, int nb_integration_pts)
Volume base interior function.
Definition: Hcurl.cpp:2468
MoFEM::TetPolynomialBase::diffN_face_edge
ublas::matrix< MatrixDouble > diffN_face_edge
Definition: TetPolynomialBase.hpp:124
MoFEM::TetPolynomialBase::getValueHdivAinsworthBaseImpl
MoFEMErrorCode getValueHdivAinsworthBaseImpl(MatrixDouble &pts, MatrixDouble &shape_functions, MatrixDouble &diff_shape_functions, int volume_order, std::array< int, 4 > &faces_order, std::array< int, 3 *4 > &faces_nodes)
Definition: TetPolynomialBase.cpp:743
NBEDGE_DEMKOWICZ_HCURL
#define NBEDGE_DEMKOWICZ_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:108
TetBaseCache::BaseCacheItem
Definition: TetPolynomialBase.cpp:12
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
MoFEM::TetPolynomialBase::N_volume_face
ublas::vector< MatrixDouble > N_volume_face
Definition: TetPolynomialBase.hpp:121
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
MoFEM::TetPolynomialBase::swichCacheHDivBaseOff
static void swichCacheHDivBaseOff(std::vector< void * > v)
Definition: TetPolynomialBase.cpp:2073
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::TetPolynomialBase::getValueHdiv
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
Get base functions for Hdiv space.
Definition: TetPolynomialBase.cpp:1599
MoFEM::TetPolynomialBase::getValueHdivAinsworthBase
MoFEMErrorCode getValueHdivAinsworthBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:861
MoFEM::EntPolynomialBaseCtx::bAse
const FieldApproximationBase bAse
Definition: EntPolynomialBaseCtx.hpp:39
TetBaseCache::BaseCacheItem::diffN
MatrixDouble diffN
Definition: TetPolynomialBase.cpp:16
MoFEM::TetPolynomialBase::diffN_volume_face
ublas::vector< MatrixDouble > diffN_volume_face
Definition: TetPolynomialBase.hpp:127
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
MoFEM::TetPolynomialBase::getValueHcurlDemkowiczBase
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:1752
FTensor::Index< 'i', 3 >
NBFACETRI_DEMKOWICZ_HCURL
#define NBFACETRI_DEMKOWICZ_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:109
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
MoFEM::TetPolynomialBase::getValueL2BernsteinBezierBase
MoFEMErrorCode getValueL2BernsteinBezierBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:567
convert.n
n
Definition: convert.py:82
MoFEM::TetPolynomialBase::swichCacheHDivBaseDemkowiczOff
static void swichCacheHDivBaseDemkowiczOff(std::vector< void * > v)
Definition: TetPolynomialBase.cpp:2032
MoFEM::UnknownInterface
base class for all interface classes
Definition: UnknownInterface.hpp:34
MoFEM::EntPolynomialBaseCtx::spaceContinuity
const FieldContinuity spaceContinuity
Definition: EntPolynomialBaseCtx.hpp:38
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
MoFEM::Tools::shapeFunMBTET
static MoFEMErrorCode shapeFunMBTET(double *shape, const double *ksi, const double *eta, const double *zeta, const double nb)
Calculate shape functions on tetrahedron.
Definition: Tools.hpp:738
MoFEM::TetPolynomialBase::getValueHdivDemkowiczBrokenBase
MoFEMErrorCode getValueHdivDemkowiczBrokenBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:1452
NBEDGE_AINSWORTH_HCURL
#define NBEDGE_AINSWORTH_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:97
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
NBFACETRI_H1
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
Definition: h1_hdiv_hcurl_l2.h:60
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
MoFEM::BernsteinBezier::generateIndicesTetTet
static MoFEMErrorCode generateIndicesTetTet(const int N, int *alpha)
Definition: BernsteinBezier.cpp:203
FieldSpaceNames
const static char *const FieldSpaceNames[]
Definition: definitions.h:92
MoFEM::TetPolynomialBase::swichCacheHDivBaseDemkowiczOn
static void swichCacheHDivBaseDemkowiczOn(std::vector< void * > v)
Definition: TetPolynomialBase.cpp:2018
MoFEM::TetPolynomialBase::cTx
EntPolynomialBaseCtx * cTx
Definition: TetPolynomialBase.hpp:45
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::TetPolynomialBase::swichCacheHdivBrokenBaseInteriorAinsworth
static bool swichCacheHdivBrokenBaseInteriorAinsworth(const void *ptr)
Definition: TetPolynomialBase.cpp:2046
MoFEM::Hcurl_Ainsworth_EdgeBaseFunctions_MBTET
MoFEMErrorCode Hcurl_Ainsworth_EdgeBaseFunctions_MBTET(int *sense, int *p, double *N, double *diffN, double *edgeN[], double *diff_edgeN[], int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Edge based H-curl base functions on tetrahedral.
Definition: Hcurl.cpp:16
MoFEM::TetPolynomialBase::swichCacheHDivBaseOn
static void swichCacheHDivBaseOn(std::vector< void * > v)
Definition: TetPolynomialBase.cpp:2068
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
MoFEM::Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET
MoFEMErrorCode Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v, double *diff_phi_v, int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Interior bubble functions by Ainsworth .
Definition: Hdiv.cpp:484
NBVOLUMETET_AINSWORTH_HCURL
#define NBVOLUMETET_AINSWORTH_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:105
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
DISCONTINUOUS
@ DISCONTINUOUS
Broken continuity (No effect on L2 space)
Definition: definitions.h:101
TetBaseCache::hdivBrokenBaseInteriorAinsworth
static std::map< const void *, BaseCacheMI > hdivBrokenBaseInteriorAinsworth
Definition: TetPolynomialBase.cpp:71
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::TetPolynomialBase::swichCacheHDivBaseAinsworthOff
static void swichCacheHDivBaseAinsworthOff(std::vector< void * > v)
Definition: TetPolynomialBase.cpp:2060
lambda
static double lambda
Definition: incompressible_elasticity.cpp:199
MoFEM::TetPolynomialBase::diffN_face_bubble
ublas::vector< MatrixDouble > diffN_face_bubble
Definition: TetPolynomialBase.hpp:125
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MoFEM::TetPolynomialBase::getValueHdivDemkowiczBase
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:1316
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::Hdiv_Ainsworth_FaceBubbleShapeFunctions
MoFEMErrorCode Hdiv_Ainsworth_FaceBubbleShapeFunctions(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f[], double *diff_phi_f[], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face bubble functions by Ainsworth .
Definition: Hdiv.cpp:137
NBVOLUMETET_AINSWORTH_HDIV
#define NBVOLUMETET_AINSWORTH_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:137
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
MoFEM::TetPolynomialBase::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: TetPolynomialBase.cpp:1903
MoFEM::EntPolynomialBaseCtx::fieldName
const std::string fieldName
Definition: EntPolynomialBaseCtx.hpp:40
MoFEM::TetPolynomialBase::swichCacheHdivBrokenBaseInteriorDemkowicz
static bool swichCacheHdivBrokenBaseInteriorDemkowicz(const void *ptr)
Definition: TetPolynomialBase.cpp:2012
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
TetBaseCache::HDivBaseCacheItem::order
int order
Definition: TetPolynomialBase.cpp:21
MoFEM::Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET
MoFEMErrorCode Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_e[6], double *diff_phi_v_e[6], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base function, Edge-based interior (volume) functions by Ainsworth .
Definition: Hdiv.cpp:280
MoFEM::Hcurl_Demkowicz_FaceBaseFunctions_MBTET
MoFEMErrorCode Hcurl_Demkowicz_FaceBaseFunctions_MBTET(int *faces_nodes, int *p, double *n, double *diff_n, double *phi[], double *diff_phi[], int nb_integration_pts)
Face base interior function.
Definition: Hcurl.cpp:2395
MoFEM::TetPolynomialBase::N_volume_bubble
MatrixDouble N_volume_bubble
Definition: TetPolynomialBase.hpp:122
MoFEM::EntitiesFieldData::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition: EntitiesFieldData.hpp:57
NBFACETRI_AINSWORTH_FACE_HDIV
#define NBFACETRI_AINSWORTH_FACE_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:131
tetCacheSwitch
auto tetCacheSwitch(const void *ptr, T &cache, std::string cache_name)
Definition: TetPolynomialBase.cpp:1985
H1_FaceShapeFunctions_MBTET
PetscErrorCode H1_FaceShapeFunctions_MBTET(int *faces_nodes, int *p, double *N, double *diffN, double *faceN[], double *diff_faceN[], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: h1.c:373
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
NBFACETRI_AINSWORTH_HCURL
#define NBFACETRI_AINSWORTH_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:100
TetBaseCache::HDivBaseCacheItem::n1
int n1
Definition: TetPolynomialBase.cpp:28
MoFEM::BernsteinBezier::baseFunctionsTet
static MoFEMErrorCode baseFunctionsTet(const int N, const int gdim, const int n_alpha, const int *alpha, const double *lambda, const double *grad_lambda, double *base, double *grad_base)
Definition: BernsteinBezier.cpp:443
H1_VolumeShapeFunctions_MBTET
PetscErrorCode H1_VolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *volumeN, double *diff_volumeN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: h1.c:475
MoFEM::EntitiesFieldData
data structure for finite element entity
Definition: EntitiesFieldData.hpp:40
MoFEM::TetPolynomialBase::getValueH1BernsteinBezierBase
MoFEMErrorCode getValueH1BernsteinBezierBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:236
NBFACETRI_AINSWORTH_HDIV
#define NBFACETRI_AINSWORTH_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:132
TetBaseCache::BaseCacheItem::order
int order
Definition: TetPolynomialBase.cpp:13
CONTINUOUS
@ CONTINUOUS
Regular field.
Definition: definitions.h:100
TetBaseCache::hdivBrokenBaseInteriorDemkowicz
static std::map< const void *, BaseCacheMI > hdivBrokenBaseInteriorDemkowicz
Definition: TetPolynomialBase.cpp:68
NBVOLUMETET_AINSWORTH_EDGE_HDIV
#define NBVOLUMETET_AINSWORTH_EDGE_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:133
convert.int
int
Definition: convert.py:64
MoFEM::Hcurl_Ainsworth_VolumeFunctions_MBTET
MoFEMErrorCode Hcurl_Ainsworth_VolumeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v, double *diff_phi_v, int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
H-curl volume base functions.
Definition: Hcurl.cpp:1403
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::TetPolynomialBase::~TetPolynomialBase
virtual ~TetPolynomialBase()
Definition: TetPolynomialBase.cpp:95
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::TetPolynomialBase::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition: TetPolynomialBase.cpp:85
MoFEM::TetPolynomialBase::swichCacheHDivBaseFaceDemkowicz
static bool swichCacheHDivBaseFaceDemkowicz(const void *ptr)
Definition: TetPolynomialBase.cpp:2002
TetBaseCache::HDivBaseCacheItem::nb_gauss_pts
int nb_gauss_pts
Definition: TetPolynomialBase.cpp:22
MoFEM::TetPolynomialBase::getValueHdivAinsworthBrokenBase
MoFEMErrorCode getValueHdivAinsworthBrokenBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:1076
MoFEM::TetPolynomialBase
Calculate base functions on tetrahedral.
Definition: TetPolynomialBase.hpp:19
MoFEM::EntPolynomialBaseCtx::basePolynomialsType0
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Definition: EntPolynomialBaseCtx.hpp:27
Generated by Doxygen 1.8.17 and hosted at University of Glasgow