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