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;
15  mutable MatrixDouble N;
17  };
18 
20 
21  int order;
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;
33  };
34 
35  using BaseCacheMI = boost::multi_index_container<
37  boost::multi_index::indexed_by<
38 
39  boost::multi_index::hashed_unique<
40 
41  composite_key<
42 
44  member<BaseCacheItem, int, &BaseCacheItem::order>,
45  member<BaseCacheItem, int, &BaseCacheItem::nb_gauss_pts>>>>
46 
47  >;
48 
49  using HDivBaseFaceCacheMI = boost::multi_index_container<
51  boost::multi_index::indexed_by<
52 
53  boost::multi_index::hashed_unique<
54 
55  composite_key<
56 
58  member<HDivBaseCacheItem, int, &HDivBaseCacheItem::order>,
59  member<HDivBaseCacheItem, int,
60  &HDivBaseCacheItem::nb_gauss_pts>,
61  member<HDivBaseCacheItem, int, &HDivBaseCacheItem::n0>,
62  member<HDivBaseCacheItem, int, &HDivBaseCacheItem::n1>,
63  member<HDivBaseCacheItem, int, &HDivBaseCacheItem::n2>>>>
64 
65  >;
66 
67  static std::map<const void *, BaseCacheMI> hdivBaseInteriorDemkowicz;
68  static std::map<const void *, HDivBaseFaceCacheMI> hDivBaseFaceDemkowicz;
69 };
70 
71 std::map<const void *, TetBaseCache::BaseCacheMI>
73 std::map<const void *, TetBaseCache::HDivBaseFaceCacheMI>
75 
77 TetPolynomialBase::query_interface(boost::typeindex::type_index type_index,
78  UnknownInterface **iface) const {
79 
81  *iface = const_cast<TetPolynomialBase *>(this);
83 }
84 
85 TetPolynomialBase::TetPolynomialBase(const void *ptr) : vPtr(ptr) {}
86 
88  if (vPtr) {
95  }
96 }
97 
100 
101  switch (cTx->bAse) {
105  break;
108  break;
109  default:
110  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
111  }
112 
114 }
115 
118 
119  EntitiesFieldData &data = cTx->dAta;
120  const FieldApproximationBase base = cTx->bAse;
121  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
122  double *diffL, const int dim) =
124 
125  int nb_gauss_pts = pts.size2();
126 
127  int sense[6], order[6];
128  if (data.spacesOnEntities[MBEDGE].test(H1)) {
129  // edges
130  if (data.dataOnEntities[MBEDGE].size() != 6) {
131  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
132  }
133  double *h1_edge_n[6], *diff_h1_egde_n[6];
134  for (int ee = 0; ee != 6; ++ee) {
135  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
136  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
137  "data inconsistency");
138  }
139  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
140  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
141  int nb_dofs = NBEDGE_H1(data.dataOnEntities[MBEDGE][ee].getOrder());
142  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
143  false);
144  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
145  3 * nb_dofs, false);
146  h1_edge_n[ee] =
147  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
148  diff_h1_egde_n[ee] =
149  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
150  }
152  sense, order,
153  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
154  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
155  h1_edge_n, diff_h1_egde_n, nb_gauss_pts, base_polynomials);
156  } else {
157  for (int ee = 0; ee != 6; ++ee) {
158  data.dataOnEntities[MBEDGE][ee].getN(base).resize(0, 0, false);
159  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(0, 0, false);
160  }
161  }
162 
163  if (data.spacesOnEntities[MBTRI].test(H1)) {
164  // faces
165  if (data.dataOnEntities[MBTRI].size() != 4) {
166  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
167  }
168  double *h1_face_n[4], *diff_h1_face_n[4];
169  for (int ff = 0; ff != 4; ++ff) {
170  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
171  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
172  "data inconsistency");
173  }
174  int nb_dofs = NBFACETRI_H1(data.dataOnEntities[MBTRI][ff].getOrder());
175  order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
176  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, nb_dofs,
177  false);
178  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
179  3 * nb_dofs, false);
180  h1_face_n[ff] =
181  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
182  diff_h1_face_n[ff] =
183  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
184  }
185  if (data.facesNodes.size1() != 4) {
186  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
187  }
188  if (data.facesNodes.size2() != 3) {
189  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
190  }
192  &*data.facesNodes.data().begin(), order,
193  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
194  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
195  h1_face_n, diff_h1_face_n, nb_gauss_pts, base_polynomials);
196 
197  } else {
198  for (int ff = 0; ff != 4; ++ff) {
199  data.dataOnEntities[MBTRI][ff].getN(base).resize(0, false);
200  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(0, 0, false);
201  }
202  }
203 
204  if (data.spacesOnEntities[MBTET].test(H1)) {
205  // volume
206  int order = data.dataOnEntities[MBTET][0].getOrder();
207  int nb_vol_dofs = NBVOLUMETET_H1(order);
208  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, nb_vol_dofs,
209  false);
210  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
211  3 * nb_vol_dofs, false);
213  data.dataOnEntities[MBTET][0].getOrder(),
214  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
215  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
216  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
217  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
218  nb_gauss_pts, base_polynomials);
219  } else {
220  data.dataOnEntities[MBTET][0].getN(base).resize(0, 0, false);
221  data.dataOnEntities[MBTET][0].getDiffN(base).resize(0, 0, false);
222  }
223 
225 }
226 
230 
231  EntitiesFieldData &data = cTx->dAta;
232  const std::string field_name = cTx->fieldName;
233  const int nb_gauss_pts = pts.size2();
234 
235  if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
236  (unsigned int)nb_gauss_pts)
237  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
238  "Base functions or nodes has wrong number of integration points "
239  "for base %s",
241  auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
242 
243  auto get_alpha = [field_name](auto &data) -> MatrixInt & {
244  auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
245  if (!ptr)
246  ptr.reset(new MatrixInt());
247  return *ptr;
248  };
249 
250  auto get_base = [field_name](auto &data) -> MatrixDouble & {
251  auto &ptr = data.getBBNSharedPtr(field_name);
252  if (!ptr)
253  ptr.reset(new MatrixDouble());
254  return *ptr;
255  };
256 
257  auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
258  auto &ptr = data.getBBDiffNSharedPtr(field_name);
259  if (!ptr)
260  ptr.reset(new MatrixDouble());
261  return *ptr;
262  };
263 
264  auto get_alpha_by_name_ptr =
265  [](auto &data,
266  const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
267  return data.getBBAlphaIndicesSharedPtr(field_name);
268  };
269 
270  auto get_base_by_name_ptr =
271  [](auto &data,
272  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
273  return data.getBBNSharedPtr(field_name);
274  };
275 
276  auto get_diff_base_by_name_ptr =
277  [](auto &data,
278  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
279  return data.getBBDiffNSharedPtr(field_name);
280  };
281 
282  auto get_alpha_by_order_ptr =
283  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
284  return data.getBBAlphaIndicesByOrderSharedPtr(o);
285  };
286 
287  auto get_base_by_order_ptr =
288  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
289  return data.getBBNByOrderSharedPtr(o);
290  };
291 
292  auto get_diff_base_by_order_ptr =
293  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
294  return data.getBBDiffNByOrderSharedPtr(o);
295  };
296 
297  auto &vert_ent_data = data.dataOnEntities[MBVERTEX][0];
298  auto &vertex_alpha = get_alpha(vert_ent_data);
299  vertex_alpha.resize(4, 4, false);
300  vertex_alpha.clear();
301  for (int n = 0; n != 4; ++n)
302  vertex_alpha(n, n) = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n];
303 
304  auto &vert_get_n = get_base(vert_ent_data);
305  auto &vert_get_diff_n = get_diff_base(vert_ent_data);
306  vert_get_n.resize(nb_gauss_pts, 4, false);
307  vert_get_diff_n.resize(nb_gauss_pts, 12, false);
309  1, lambda.size1(), vertex_alpha.size1(), &vertex_alpha(0, 0),
310  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &vert_get_n(0, 0),
311  &vert_get_diff_n(0, 0));
312  for (int n = 0; n != 4; ++n) {
313  const double f = boost::math::factorial<double>(
314  data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n]);
315  for (int g = 0; g != nb_gauss_pts; ++g) {
316  vert_get_n(g, n) *= f;
317  for (int d = 0; d != 3; ++d)
318  vert_get_diff_n(g, 3 * n + d) *= f;
319  }
320  }
321 
322  // edges
323  if (data.spacesOnEntities[MBEDGE].test(H1)) {
324  if (data.dataOnEntities[MBEDGE].size() != 6)
325  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
326  "Wrong size of ent data");
327 
328  constexpr int edges_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
329  {0, 3}, {1, 3}, {2, 3}};
330  for (int ee = 0; ee != 6; ++ee) {
331  auto &ent_data = data.dataOnEntities[MBEDGE][ee];
332  const int sense = ent_data.getSense();
333  if (sense == 0)
334  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
335  "Sense of the edge unknown");
336  const int order = ent_data.getOrder();
337  const int nb_dofs = NBEDGE_H1(order);
338 
339  if (nb_dofs) {
340  if (get_alpha_by_order_ptr(ent_data, order)) {
341  get_alpha_by_name_ptr(ent_data, field_name) =
342  get_alpha_by_order_ptr(ent_data, order);
343  get_base_by_name_ptr(ent_data, field_name) =
344  get_base_by_order_ptr(ent_data, order);
345  get_diff_base_by_name_ptr(ent_data, field_name) =
346  get_diff_base_by_order_ptr(ent_data, order);
347  } else {
348  auto &get_n = get_base(ent_data);
349  auto &get_diff_n = get_diff_base(ent_data);
350  get_n.resize(nb_gauss_pts, nb_dofs, false);
351  get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
352 
353  auto &edge_alpha = get_alpha(data.dataOnEntities[MBEDGE][ee]);
354  edge_alpha.resize(nb_dofs, 4, false);
356  &edge_alpha(0, 0));
357  if (sense == -1) {
358  for (int i = 0; i != edge_alpha.size1(); ++i) {
359  int a = edge_alpha(i, edges_nodes[ee][0]);
360  edge_alpha(i, edges_nodes[ee][0]) =
361  edge_alpha(i, edges_nodes[ee][1]);
362  edge_alpha(i, edges_nodes[ee][1]) = a;
363  }
364  }
366  order, lambda.size1(), edge_alpha.size1(), &edge_alpha(0, 0),
367  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
368  &get_diff_n(0, 0));
369 
370  get_alpha_by_order_ptr(ent_data, order) =
371  get_alpha_by_name_ptr(ent_data, field_name);
372  get_base_by_order_ptr(ent_data, order) =
373  get_base_by_name_ptr(ent_data, field_name);
374  get_diff_base_by_order_ptr(ent_data, order) =
375  get_diff_base_by_name_ptr(ent_data, field_name);
376  }
377  }
378  }
379  } else {
380  for (int ee = 0; ee != 6; ++ee) {
381  auto &ent_data = data.dataOnEntities[MBEDGE][ee];
382  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
383  auto &get_n = get_base(ent_data);
384  auto &get_diff_n = get_diff_base(ent_data);
385  get_n.resize(nb_gauss_pts, 0, false);
386  get_diff_n.resize(nb_gauss_pts, 0, false);
387  }
388  }
389 
390  // face
391  if (data.spacesOnEntities[MBTRI].test(H1)) {
392  if (data.dataOnEntities[MBTRI].size() != 4)
393  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
394  "Wrong size of ent data");
395  if (data.facesNodes.size1() != 4)
396  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
397  if (data.facesNodes.size2() != 3)
398  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
399 
400  for (int ff = 0; ff != 4; ++ff) {
401  auto &ent_data = data.dataOnEntities[MBTRI][ff];
402  const int order = ent_data.getOrder();
403  const int nb_dofs = NBFACETRI_H1(order);
404 
405  if (nb_dofs) {
406  if (get_alpha_by_order_ptr(ent_data, order)) {
407  get_alpha_by_name_ptr(ent_data, field_name) =
408  get_alpha_by_order_ptr(ent_data, order);
409  get_base_by_name_ptr(ent_data, field_name) =
410  get_base_by_order_ptr(ent_data, order);
411  get_diff_base_by_name_ptr(ent_data, field_name) =
412  get_diff_base_by_order_ptr(ent_data, order);
413  } else {
414 
415  auto &get_n = get_base(ent_data);
416  auto &get_diff_n = get_diff_base(ent_data);
417  get_n.resize(nb_gauss_pts, nb_dofs, false);
418  get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
419 
420  auto &face_alpha = get_alpha(ent_data);
421  face_alpha.resize(nb_dofs, 4, false);
422 
424  &face_alpha(0, 0));
425  senseFaceAlpha.resize(face_alpha.size1(), face_alpha.size2(), false);
426  senseFaceAlpha.clear();
427  constexpr int tri_nodes[4][3] = {
428  {0, 1, 3}, {1, 2, 3}, {0, 2, 3}, {0, 1, 2}};
429  for (int d = 0; d != nb_dofs; ++d)
430  for (int n = 0; n != 3; ++n)
431  senseFaceAlpha(d, data.facesNodes(ff, n)) =
432  face_alpha(d, tri_nodes[ff][n]);
433  face_alpha.swap(senseFaceAlpha);
435  order, lambda.size1(), face_alpha.size1(), &face_alpha(0, 0),
436  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
437  &get_diff_n(0, 0));
438 
439  get_alpha_by_order_ptr(ent_data, order) =
440  get_alpha_by_name_ptr(ent_data, field_name);
441  get_base_by_order_ptr(ent_data, order) =
442  get_base_by_name_ptr(ent_data, field_name);
443  get_diff_base_by_order_ptr(ent_data, order) =
444  get_diff_base_by_name_ptr(ent_data, field_name);
445  }
446  }
447  }
448  } else {
449  for (int ff = 0; ff != 4; ++ff) {
450  auto &ent_data = data.dataOnEntities[MBTRI][ff];
451  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
452  auto &get_n = get_base(ent_data);
453  auto &get_diff_n = get_diff_base(ent_data);
454  get_n.resize(nb_gauss_pts, 0, false);
455  get_diff_n.resize(nb_gauss_pts, 0, false);
456  }
457  }
458 
459  if (data.spacesOnEntities[MBTET].test(H1)) {
460  if (data.dataOnEntities[MBTET].size() != 1)
461  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
462  "Wrong size ent of ent data");
463 
464  auto &ent_data = data.dataOnEntities[MBTET][0];
465  const int order = ent_data.getOrder();
466  const int nb_dofs = NBVOLUMETET_H1(order);
467  if (get_alpha_by_order_ptr(ent_data, order)) {
468  get_alpha_by_name_ptr(ent_data, field_name) =
469  get_alpha_by_order_ptr(ent_data, order);
470  get_base_by_name_ptr(ent_data, field_name) =
471  get_base_by_order_ptr(ent_data, order);
472  get_diff_base_by_name_ptr(ent_data, field_name) =
473  get_diff_base_by_order_ptr(ent_data, order);
474  } else {
475 
476  auto &get_n = get_base(ent_data);
477  auto &get_diff_n = get_diff_base(ent_data);
478  get_n.resize(nb_gauss_pts, nb_dofs, false);
479  get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
480  if (nb_dofs) {
481  auto &tet_alpha = get_alpha(ent_data);
482  tet_alpha.resize(nb_dofs, 4, false);
483 
486  order, lambda.size1(), tet_alpha.size1(), &tet_alpha(0, 0),
487  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
488  &get_diff_n(0, 0));
489 
490  get_alpha_by_order_ptr(ent_data, order) =
491  get_alpha_by_name_ptr(ent_data, field_name);
492  get_base_by_order_ptr(ent_data, order) =
493  get_base_by_name_ptr(ent_data, field_name);
494  get_diff_base_by_order_ptr(ent_data, order) =
495  get_diff_base_by_name_ptr(ent_data, field_name);
496  }
497  }
498  } else {
499  auto &ent_data = data.dataOnEntities[MBTET][0];
500  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
501  auto &get_n = get_base(ent_data);
502  auto &get_diff_n = get_diff_base(ent_data);
503  get_n.resize(nb_gauss_pts, 0, false);
504  get_diff_n.resize(nb_gauss_pts, 0, false);
505  }
506 
508 }
509 
512 
513  switch (cTx->bAse) {
517  break;
520  break;
521  default:
522  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
523  }
524 
526 }
527 
530 
531  EntitiesFieldData &data = cTx->dAta;
532  const FieldApproximationBase base = cTx->bAse;
533  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
534  double *diffL, const int dim) =
536 
537  int nb_gauss_pts = pts.size2();
538 
539  data.dataOnEntities[MBTET][0].getN(base).resize(
540  nb_gauss_pts, NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getOrder()),
541  false);
542  data.dataOnEntities[MBTET][0].getDiffN(base).resize(
543  nb_gauss_pts,
544  3 * NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getOrder()), false);
545 
547  data.dataOnEntities[MBTET][0].getOrder(),
548  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
549  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
550  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
551  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
552  nb_gauss_pts, base_polynomials);
553 
555 }
556 
560 
561  EntitiesFieldData &data = cTx->dAta;
562  const std::string field_name = cTx->fieldName;
563  const int nb_gauss_pts = pts.size2();
564 
565  if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
566  (unsigned int)nb_gauss_pts)
567  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
568  "Base functions or nodes has wrong number of integration points "
569  "for base %s",
571  auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
572 
573  auto get_alpha = [field_name](auto &data) -> MatrixInt & {
574  auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
575  if (!ptr)
576  ptr.reset(new MatrixInt());
577  return *ptr;
578  };
579 
580  auto get_base = [field_name](auto &data) -> MatrixDouble & {
581  auto &ptr = data.getBBNSharedPtr(field_name);
582  if (!ptr)
583  ptr.reset(new MatrixDouble());
584  return *ptr;
585  };
586 
587  auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
588  auto &ptr = data.getBBDiffNSharedPtr(field_name);
589  if (!ptr)
590  ptr.reset(new MatrixDouble());
591  return *ptr;
592  };
593 
594  auto get_alpha_by_name_ptr =
595  [](auto &data,
596  const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
597  return data.getBBAlphaIndicesSharedPtr(field_name);
598  };
599 
600  auto get_base_by_name_ptr =
601  [](auto &data,
602  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
603  return data.getBBNSharedPtr(field_name);
604  };
605 
606  auto get_diff_base_by_name_ptr =
607  [](auto &data,
608  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
609  return data.getBBDiffNSharedPtr(field_name);
610  };
611 
612  auto get_alpha_by_order_ptr =
613  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
614  return data.getBBAlphaIndicesByOrderSharedPtr(o);
615  };
616 
617  auto get_base_by_order_ptr =
618  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
619  return data.getBBNByOrderSharedPtr(o);
620  };
621 
622  auto get_diff_base_by_order_ptr =
623  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
624  return data.getBBDiffNByOrderSharedPtr(o);
625  };
626 
627  if (data.spacesOnEntities[MBTET].test(L2)) {
628  if (data.dataOnEntities[MBTET].size() != 1)
629  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
630  "Wrong size ent of ent data");
631 
632  auto &ent_data = data.dataOnEntities[MBTET][0];
633  const int order = ent_data.getOrder();
634  const int nb_dofs = NBVOLUMETET_L2(order);
635 
636  if (get_alpha_by_order_ptr(ent_data, order)) {
637  get_alpha_by_name_ptr(ent_data, field_name) =
638  get_alpha_by_order_ptr(ent_data, order);
639  get_base_by_name_ptr(ent_data, field_name) =
640  get_base_by_order_ptr(ent_data, order);
641  get_diff_base_by_name_ptr(ent_data, field_name) =
642  get_diff_base_by_order_ptr(ent_data, order);
643  } else {
644 
645  auto &get_n = get_base(ent_data);
646  auto &get_diff_n = get_diff_base(ent_data);
647  get_n.resize(nb_gauss_pts, nb_dofs, false);
648  get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
649 
650  if (nb_dofs) {
651 
652  if (order == 0) {
653 
654  if (nb_dofs != 1)
655  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
656  "Inconsistent number of DOFs");
657 
658  auto &tri_alpha = get_alpha(ent_data);
659  tri_alpha.clear();
660  get_n(0, 0) = 1;
661  get_diff_n.clear();
662 
663  } else {
664 
665  if (nb_dofs != 4 + 6 * NBEDGE_H1(order) + 4 * NBFACETRI_H1(order) +
667  nb_dofs != NBVOLUMETET_L2(order))
668  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
669  "Inconsistent number of DOFs");
670 
671  auto &tet_alpha = get_alpha(ent_data);
672  tet_alpha.resize(nb_dofs, 4, false);
673 
675  &tet_alpha(0, 0));
676  if (order > 1) {
677  std::array<int, 6> edge_n{order, order, order, order, order, order};
678  std::array<int *, 6> tet_edge_ptr{
679  &tet_alpha(4, 0),
680  &tet_alpha(4 + 1 * NBEDGE_H1(order), 0),
681  &tet_alpha(4 + 2 * NBEDGE_H1(order), 0),
682  &tet_alpha(4 + 3 * NBEDGE_H1(order), 0),
683  &tet_alpha(4 + 4 * NBEDGE_H1(order), 0),
684  &tet_alpha(4 + 5 * NBEDGE_H1(order), 0)};
686  tet_edge_ptr.data());
687  if (order > 2) {
688  std::array<int, 6> face_n{order, order, order, order};
689  std::array<int *, 6> tet_face_ptr{
690  &tet_alpha(4 + 6 * NBEDGE_H1(order), 0),
691  &tet_alpha(4 + 6 * NBEDGE_H1(order) + 1 * NBFACETRI_H1(order),
692  0),
693  &tet_alpha(4 + 6 * NBEDGE_H1(order) + 2 * NBFACETRI_H1(order),
694  0),
695  &tet_alpha(4 + 6 * NBEDGE_H1(order) + 3 * NBFACETRI_H1(order),
696  0),
697  };
699  face_n.data(), tet_face_ptr.data());
700  if (order > 3)
702  order,
703  &tet_alpha(
704  4 + 6 * NBEDGE_H1(order) + 4 * NBFACETRI_H1(order), 0));
705  }
706  }
707 
709  order, lambda.size1(), tet_alpha.size1(), &tet_alpha(0, 0),
710  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
711  &get_diff_n(0, 0));
712 
713  get_alpha_by_order_ptr(ent_data, order) =
714  get_alpha_by_name_ptr(ent_data, field_name);
715  get_base_by_order_ptr(ent_data, order) =
716  get_base_by_name_ptr(ent_data, field_name);
717  get_diff_base_by_order_ptr(ent_data, order) =
718  get_diff_base_by_name_ptr(ent_data, field_name);
719  }
720  }
721  }
722  } else {
723  auto &ent_data = data.dataOnEntities[MBTET][0];
724  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
725  auto &get_n = get_base(ent_data);
726  auto &get_diff_n = get_diff_base(ent_data);
727  get_n.resize(nb_gauss_pts, 0, false);
728  get_diff_n.resize(nb_gauss_pts, 0, false);
729  }
730 
732 }
733 
736 
737  EntitiesFieldData &data = cTx->dAta;
738  const FieldApproximationBase base = cTx->bAse;
739  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
740  double *diffL, const int dim) =
742 
743  int nb_gauss_pts = pts.size2();
744 
745  // face shape functions
746 
747  double *phi_f_e[4][3];
748  double *phi_f[4];
749  double *diff_phi_f_e[4][3];
750  double *diff_phi_f[4];
751 
752  N_face_edge.resize(4, 3, false);
753  N_face_bubble.resize(4, false);
754  diffN_face_edge.resize(4, 3, false);
755  diffN_face_bubble.resize(4, false);
756 
757  int faces_order[4];
758  for (int ff = 0; ff != 4; ++ff) {
759  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
760  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
761  }
762  faces_order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
763  // three edges on face
764  for (int ee = 0; ee < 3; ee++) {
765  N_face_edge(ff, ee).resize(
766  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(faces_order[ff]),
767  false);
768  diffN_face_edge(ff, ee).resize(
769  nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_EDGE_HDIV(faces_order[ff]),
770  false);
771  phi_f_e[ff][ee] = &*N_face_edge(ff, ee).data().begin();
772  diff_phi_f_e[ff][ee] = &*diffN_face_edge(ff, ee).data().begin();
773  }
774  N_face_bubble[ff].resize(nb_gauss_pts,
775  3 * NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]),
776  false);
777  diffN_face_bubble[ff].resize(
778  nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]),
779  false);
780  phi_f[ff] = &*(N_face_bubble[ff].data().begin());
781  diff_phi_f[ff] = &*(diffN_face_bubble[ff].data().begin());
782  }
783 
785  &data.facesNodes(0, 0), faces_order,
786  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
787  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(), phi_f_e,
788  diff_phi_f_e, nb_gauss_pts, base_polynomials);
789 
791  &data.facesNodes(0, 0), faces_order,
792  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
793  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(), phi_f,
794  diff_phi_f, nb_gauss_pts, base_polynomials);
795 
796  // volume shape functions
797 
798  double *phi_v_e[6];
799  double *phi_v_f[4];
800  double *phi_v;
801  double *diff_phi_v_e[6];
802  double *diff_phi_v_f[4];
803  double *diff_phi_v;
804 
805  int volume_order = data.dataOnEntities[MBTET][0].getOrder();
806 
807  N_volume_edge.resize(6, false);
808  diffN_volume_edge.resize(6, false);
809  for (int ee = 0; ee != 6; ++ee) {
810  N_volume_edge[ee].resize(
811  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_EDGE_HDIV(volume_order), false);
812  diffN_volume_edge[ee].resize(
813  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_EDGE_HDIV(volume_order), false);
814  phi_v_e[ee] = &*(N_volume_edge[ee].data().begin());
815  diff_phi_v_e[ee] = &*(diffN_volume_edge[ee].data().begin());
816  }
818  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
819  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_v_e,
820  diff_phi_v_e, nb_gauss_pts, base_polynomials);
821 
822  N_volume_face.resize(4, false);
823  diffN_volume_face.resize(4, false);
824  for (int ff = 0; ff != 4; ++ff) {
825  N_volume_face[ff].resize(
826  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order), false);
827  diffN_volume_face[ff].resize(
828  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order), false);
829  phi_v_f[ff] = &*(N_volume_face[ff].data().begin());
830  diff_phi_v_f[ff] = &*(diffN_volume_face[ff].data().begin());
831  }
833  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
834  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_v_f,
835  diff_phi_v_f, nb_gauss_pts, base_polynomials);
836 
837  N_volume_bubble.resize(
838  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order), false);
839  diffN_volume_bubble.resize(
840  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order), false);
841  phi_v = &*(N_volume_bubble.data().begin());
842  diff_phi_v = &*(diffN_volume_bubble.data().begin());
844  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
845  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_v, diff_phi_v,
846  nb_gauss_pts, base_polynomials);
847 
848  // Set shape functions into data structure Shape functions hast to be put
849  // in arrays in order which guarantee hierarchical series of degrees of
850  // freedom, i.e. in other words dofs form sub-entities has to be group
851  // by order.
852 
855 
856  // faces
857  if (data.dataOnEntities[MBTRI].size() != 4) {
858  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
859  }
860  for (int ff = 0; ff != 4; ff++) {
861  data.dataOnEntities[MBTRI][ff].getN(base).resize(
862  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_HDIV(faces_order[ff]), false);
863  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
864  nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_HDIV(faces_order[ff]), false);
865  if (NBFACETRI_AINSWORTH_HDIV(faces_order[ff]) == 0)
866  continue;
867  // face
868  double *base_ptr =
869  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
870  FTensor::Tensor1<double *, 3> t_base(base_ptr, &base_ptr[HVEC1],
871  &base_ptr[HVEC2], 3);
872  double *diff_base_ptr =
873  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
875  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
876  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
877  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
878  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
879  &diff_base_ptr[HVEC2_2], 9);
880  // face-face
881  boost::shared_ptr<FTensor::Tensor1<double *, 3>> t_base_f;
882  boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>> t_diff_base_f;
883  if (NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]) > 0) {
884  base_ptr = phi_f[ff];
885  t_base_f = boost::shared_ptr<FTensor::Tensor1<double *, 3>>(
886  new FTensor::Tensor1<double *, 3>(base_ptr, &base_ptr[HVEC1],
887  &base_ptr[HVEC2], 3));
888  diff_base_ptr = diff_phi_f[ff];
889  t_diff_base_f = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>>(
891  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
892  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
893  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
894  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
895  &diff_base_ptr[HVEC2_2], 9));
896  }
897  // edge-face
898  base_ptr = phi_f_e[ff][0];
899  FTensor::Tensor1<double *, 3> t_base_f_e0(base_ptr, &base_ptr[HVEC1],
900  &base_ptr[HVEC2], 3);
901  diff_base_ptr = diff_phi_f_e[ff][0];
902  FTensor::Tensor2<double *, 3, 3> t_diff_base_f_e0(
903  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
904  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
905  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
906  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
907  &diff_base_ptr[HVEC2_2], 9);
908  base_ptr = phi_f_e[ff][1];
909  FTensor::Tensor1<double *, 3> t_base_f_e1(base_ptr, &base_ptr[HVEC1],
910  &base_ptr[HVEC2], 3);
911  diff_base_ptr = diff_phi_f_e[ff][1];
912  FTensor::Tensor2<double *, 3, 3> t_diff_base_f_e1(
913  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
914  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
915  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
916  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
917  &diff_base_ptr[HVEC2_2], 9);
918  base_ptr = phi_f_e[ff][2];
919  FTensor::Tensor1<double *, 3> t_base_f_e2(base_ptr, &base_ptr[HVEC1],
920  &base_ptr[HVEC2], 3);
921  diff_base_ptr = diff_phi_f_e[ff][2];
922  FTensor::Tensor2<double *, 3, 3> t_diff_base_f_e2(
923  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
924  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
925  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
926  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
927  &diff_base_ptr[HVEC2_2], 9);
928  for (int gg = 0; gg != nb_gauss_pts; gg++) {
929  for (int oo = 0; oo != faces_order[ff]; oo++) {
930  for (int dd = NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
931  dd != NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
932  t_base(i) = t_base_f_e0(i);
933  ++t_base;
934  ++t_base_f_e0;
935  t_diff_base(i, j) = t_diff_base_f_e0(i, j);
936  ++t_diff_base;
937  ++t_diff_base_f_e0;
938  t_base(i) = t_base_f_e1(i);
939  ++t_base;
940  ++t_base_f_e1;
941  t_diff_base(i, j) = t_diff_base_f_e1(i, j);
942  ++t_diff_base;
943  ++t_diff_base_f_e1;
944  t_base(i) = t_base_f_e2(i);
945  ++t_base;
946  ++t_base_f_e2;
947  t_diff_base(i, j) = t_diff_base_f_e2(i, j);
948  ++t_diff_base;
949  ++t_diff_base_f_e2;
950  }
951  for (int dd = NBFACETRI_AINSWORTH_FACE_HDIV(oo);
952  dd != NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
953  t_base(i) = (*t_base_f)(i);
954  ++t_base;
955  ++(*t_base_f);
956  t_diff_base(i, j) = (*t_diff_base_f)(i, j);
957  ++t_diff_base;
958  ++(*t_diff_base_f);
959  }
960  }
961  }
962  }
963 
964  // volume
965  data.dataOnEntities[MBTET][0].getN(base).resize(
966  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_HDIV(volume_order), false);
967  data.dataOnEntities[MBTET][0].getDiffN(base).resize(
968  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_HDIV(volume_order), false);
969  if (NBVOLUMETET_AINSWORTH_HDIV(volume_order) > 0) {
970  double *base_ptr =
971  &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
972  FTensor::Tensor1<double *, 3> t_base(base_ptr, &base_ptr[HVEC1],
973  &base_ptr[HVEC2], 3);
974  double *diff_base_ptr =
975  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
977  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
978  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
979  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
980  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
981  &diff_base_ptr[HVEC2_2], 9);
982  // edges
983  std::vector<FTensor::Tensor1<double *, 3>> t_base_v_e;
984  t_base_v_e.reserve(6);
985  std::vector<FTensor::Tensor2<double *, 3, 3>> t_diff_base_v_e;
986  t_diff_base_v_e.reserve(6);
987  for (int ee = 0; ee != 6; ee++) {
988  base_ptr = phi_v_e[ee];
989  diff_base_ptr = diff_phi_v_e[ee];
990  t_base_v_e.push_back(FTensor::Tensor1<double *, 3>(
991  base_ptr, &base_ptr[HVEC1], &base_ptr[HVEC2], 3));
992  t_diff_base_v_e.push_back(FTensor::Tensor2<double *, 3, 3>(
993  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
994  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
995  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
996  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
997  &diff_base_ptr[HVEC2_2], 9));
998  }
999  // faces
1000  std::vector<FTensor::Tensor1<double *, 3>> t_base_v_f;
1001  t_base_v_f.reserve(4);
1002  std::vector<FTensor::Tensor2<double *, 3, 3>> t_diff_base_v_f;
1003  t_diff_base_v_f.reserve(4);
1004  if (NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order) > 0) {
1005  for (int ff = 0; ff != 4; ff++) {
1006  base_ptr = phi_v_f[ff];
1007  diff_base_ptr = diff_phi_v_f[ff];
1008  t_base_v_f.push_back(FTensor::Tensor1<double *, 3>(
1009  base_ptr, &base_ptr[HVEC1], &base_ptr[HVEC2], 3));
1010  t_diff_base_v_f.push_back(FTensor::Tensor2<double *, 3, 3>(
1011  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
1012  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
1013  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
1014  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
1015  &diff_base_ptr[HVEC2_2], 9));
1016  }
1017  }
1018  boost::shared_ptr<FTensor::Tensor1<double *, 3>> t_base_v;
1019  boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>> t_diff_base_v;
1020  if (NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order) > 0) {
1021  base_ptr = phi_v;
1022  t_base_v = boost::shared_ptr<FTensor::Tensor1<double *, 3>>(
1023  new FTensor::Tensor1<double *, 3>(base_ptr, &base_ptr[HVEC1],
1024  &base_ptr[HVEC2], 3));
1025  diff_base_ptr = diff_phi_v;
1026  t_diff_base_v = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>>(
1028  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
1029  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
1030  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
1031  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
1032  &diff_base_ptr[HVEC2_2], 9));
1033  }
1034  for (int gg = 0; gg != nb_gauss_pts; gg++) {
1035  for (int oo = 0; oo < volume_order; oo++) {
1036  for (int dd = NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo);
1037  dd < NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
1038  for (int ee = 0; ee < 6; ee++) {
1039  t_base(i) = t_base_v_e[ee](i);
1040  ++t_base;
1041  ++t_base_v_e[ee];
1042  t_diff_base(i, j) = t_diff_base_v_e[ee](i, j);
1043  ++t_diff_base;
1044  ++t_diff_base_v_e[ee];
1045  }
1046  }
1047  for (int dd = NBVOLUMETET_AINSWORTH_FACE_HDIV(oo);
1048  dd < NBVOLUMETET_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
1049  for (int ff = 0; ff < 4; ff++) {
1050  t_base(i) = t_base_v_f[ff](i);
1051  ++t_base;
1052  ++t_base_v_f[ff];
1053  t_diff_base(i, j) = t_diff_base_v_f[ff](i, j);
1054  ++t_diff_base;
1055  ++t_diff_base_v_f[ff];
1056  }
1057  }
1058  for (int dd = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo);
1059  dd < NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo + 1); dd++) {
1060  t_base(i) = (*t_base_v)(i);
1061  ++t_base;
1062  ++(*t_base_v);
1063  t_diff_base(i, j) = (*t_diff_base_v)(i, j);
1064  ++t_diff_base;
1065  ++(*t_diff_base_v);
1066  }
1067  }
1068  }
1069  }
1070 
1072 }
1073 
1076 
1077  EntitiesFieldData &data = cTx->dAta;
1078  const FieldApproximationBase base = cTx->bAse;
1079  if (base != DEMKOWICZ_JACOBI_BASE) {
1080  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1081  "This should be used only with DEMKOWICZ_JACOBI_BASE "
1082  "but base is %s",
1083  ApproximationBaseNames[base]);
1084  }
1085  int nb_gauss_pts = pts.size2();
1086 
1087  int volume_order = data.dataOnEntities[MBTET][0].getOrder();
1088 
1089  int p_f[4];
1090  double *phi_f[4];
1091  double *diff_phi_f[4];
1092 
1093  auto get_face_cache_ptr = [this]() -> TetBaseCache::HDivBaseFaceCacheMI * {
1094  if (vPtr) {
1095  auto it = TetBaseCache::hDivBaseFaceDemkowicz.find(vPtr);
1096  if (it != TetBaseCache::hDivBaseFaceDemkowicz.end()) {
1097  return &it->second;
1098  }
1099  }
1100  return nullptr;
1101  };
1102 
1103  auto face_cache_ptr = get_face_cache_ptr();
1104 
1105  // Calculate base function on tet faces
1106  for (int ff = 0; ff != 4; ff++) {
1107  int face_order = data.dataOnEntities[MBTRI][ff].getOrder();
1108  int order = volume_order > face_order ? volume_order : face_order;
1109  data.dataOnEntities[MBTRI][ff].getN(base).resize(
1110  nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
1111  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
1112  nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
1113  if (NBFACETRI_DEMKOWICZ_HDIV(order) == 0)
1114  continue;
1115 
1116  if (face_cache_ptr) {
1117  auto it = face_cache_ptr->find(boost::make_tuple(
1118 
1119  face_order, nb_gauss_pts,
1120 
1121  data.facesNodes(ff, 0), data.facesNodes(ff, 1), data.facesNodes(ff, 2)
1122 
1123  ));
1124  if (it != face_cache_ptr->end()) {
1125  noalias(data.dataOnEntities[MBTRI][ff].getN(base)) = it->N;
1126  noalias(data.dataOnEntities[MBTRI][ff].getDiffN(base)) = it->diffN;
1127  continue;
1128  }
1129  }
1130 
1131  p_f[ff] = order;
1132  phi_f[ff] = &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1133  diff_phi_f[ff] =
1134  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1135 
1137  &data.facesNodes(ff, 0), order,
1138  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1139  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1140  phi_f[ff], diff_phi_f[ff], nb_gauss_pts, 4);
1141  if (face_cache_ptr) {
1142  auto p = face_cache_ptr->emplace(TetBaseCache::HDivBaseCacheItem{
1143  face_order, nb_gauss_pts, data.facesNodes(ff, 0),
1144  data.facesNodes(ff, 1), data.facesNodes(ff, 2)});
1145  p.first->N = data.dataOnEntities[MBTRI][ff].getN(base);
1146  p.first->diffN = data.dataOnEntities[MBTRI][ff].getDiffN(base);
1147  }
1148  }
1149 
1150  auto get_interior_cache = [this]() -> TetBaseCache::BaseCacheMI * {
1151  if (vPtr) {
1153  if (it != TetBaseCache::hdivBaseInteriorDemkowicz.end()) {
1154  return &it->second;
1155  }
1156  }
1157  return nullptr;
1158  };
1159 
1160  auto interior_cache_ptr = get_interior_cache();
1161 
1162  // Calculate base functions in tet interior
1163  if (NBVOLUMETET_DEMKOWICZ_HDIV(volume_order) > 0) {
1164  data.dataOnEntities[MBTET][0].getN(base).resize(
1165  nb_gauss_pts, 3 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
1166  data.dataOnEntities[MBTET][0].getDiffN(base).resize(
1167  nb_gauss_pts, 9 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
1168 
1169  for (int v = 0; v != 1; ++v) {
1170  if (interior_cache_ptr) {
1171  auto it = interior_cache_ptr->find(
1172  boost::make_tuple(volume_order, nb_gauss_pts));
1173  if (it != interior_cache_ptr->end()) {
1174  noalias(data.dataOnEntities[MBTET][0].getN(base)) = it->N;
1175  noalias(data.dataOnEntities[MBTET][0].getDiffN(base)) = it->diffN;
1176  continue;
1177  }
1178  }
1179 
1180  double *phi_v = &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
1181  double *diff_phi_v =
1182  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
1183 
1185  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
1186  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), p_f, phi_f,
1187  diff_phi_f, phi_v, diff_phi_v, nb_gauss_pts);
1188  if (interior_cache_ptr) {
1189  auto p = interior_cache_ptr->emplace(
1190  TetBaseCache::BaseCacheItem{volume_order, nb_gauss_pts});
1191  p.first->N = data.dataOnEntities[MBTET][0].getN(base);
1192  p.first->diffN = data.dataOnEntities[MBTET][0].getDiffN(base);
1193  }
1194  }
1195  }
1196 
1197  // Set size of face base correctly
1198  for (int ff = 0; ff != 4; ff++) {
1199  int face_order = data.dataOnEntities[MBTRI][ff].getOrder();
1200  data.dataOnEntities[MBTRI][ff].getN(base).resize(
1201  nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
1202  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
1203  nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
1204  }
1205 
1207 }
1208 
1211 
1212  switch (cTx->bAse) {
1215  return getValueHdivAinsworthBase(pts);
1216  case DEMKOWICZ_JACOBI_BASE:
1217  return getValueHdivDemkowiczBase(pts);
1218  default:
1219  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1220  }
1221 
1223 }
1224 
1228 
1229  EntitiesFieldData &data = cTx->dAta;
1230  const FieldApproximationBase base = cTx->bAse;
1231  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
1232  double *diffL, const int dim) =
1234 
1235  int nb_gauss_pts = pts.size2();
1236 
1237  // edges
1238  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
1239  int sense[6], order[6];
1240  if (data.dataOnEntities[MBEDGE].size() != 6) {
1241  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1242  }
1243  double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
1244  for (int ee = 0; ee != 6; ee++) {
1245  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
1246  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1247  "data inconsistency");
1248  }
1249  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
1250  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
1251  int nb_dofs =
1252  NBEDGE_AINSWORTH_HCURL(data.dataOnEntities[MBEDGE][ee].getOrder());
1253  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
1254  3 * nb_dofs, false);
1255  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
1256  9 * nb_dofs, false);
1257  hcurl_edge_n[ee] =
1258  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
1259  diff_hcurl_edge_n[ee] =
1260  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
1261  }
1263  sense, order,
1264  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1265  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1266  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
1267  } else {
1268  for (int ee = 0; ee != 6; ee++) {
1269  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
1270  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
1271  false);
1272  }
1273  }
1274 
1275  // triangles
1276  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
1277  int order[4];
1278  // faces
1279  if (data.dataOnEntities[MBTRI].size() != 4) {
1280  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1281  }
1282  double *hcurl_base_n[4], *diff_hcurl_base_n[4];
1283  for (int ff = 0; ff != 4; ff++) {
1284  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
1285  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1286  "data inconsistency");
1287  }
1288  order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
1289  int nb_dofs = NBFACETRI_AINSWORTH_HCURL(order[ff]);
1290  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1291  3 * nb_dofs, false);
1292  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1293  9 * nb_dofs, false);
1294  hcurl_base_n[ff] =
1295  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1296  diff_hcurl_base_n[ff] =
1297  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1298  }
1299  if (data.facesNodes.size1() != 4) {
1300  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1301  }
1302  if (data.facesNodes.size2() != 3) {
1303  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
1304  }
1306  &*data.facesNodes.data().begin(), order,
1307  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1308  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1309  hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts, base_polynomials);
1310  } else {
1311  for (int ff = 0; ff != 4; ff++) {
1312  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
1313  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
1314  false);
1315  }
1316  }
1317 
1318  if (data.spacesOnEntities[MBTET].test(HCURL)) {
1319 
1320  // volume
1321  int order = data.dataOnEntities[MBTET][0].getOrder();
1322  int nb_vol_dofs = NBVOLUMETET_AINSWORTH_HCURL(order);
1323  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
1324  3 * nb_vol_dofs, false);
1325  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
1326  9 * nb_vol_dofs, false);
1328  data.dataOnEntities[MBTET][0].getOrder(),
1329  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1330  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1331  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
1332  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
1333  nb_gauss_pts, base_polynomials);
1334 
1335  } else {
1336  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
1337  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
1338  }
1339 
1341 }
1342 
1346 
1347  EntitiesFieldData &data = cTx->dAta;
1348  const FieldApproximationBase base = cTx->bAse;
1349  if (base != DEMKOWICZ_JACOBI_BASE) {
1350  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1351  "This should be used only with DEMKOWICZ_JACOBI_BASE "
1352  "but base is %s",
1353  ApproximationBaseNames[base]);
1354  }
1355 
1356  int nb_gauss_pts = pts.size2();
1357 
1358  // edges
1359  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
1360  int sense[6], order[6];
1361  if (data.dataOnEntities[MBEDGE].size() != 6) {
1362  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1363  "wrong size of data structure, expected space for six edges "
1364  "but is %d",
1365  data.dataOnEntities[MBEDGE].size());
1366  }
1367  double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
1368  for (int ee = 0; ee != 6; ee++) {
1369  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
1370  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1371  "orintation of edges is not set");
1372  }
1373  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
1374  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
1375  int nb_dofs =
1376  NBEDGE_DEMKOWICZ_HCURL(data.dataOnEntities[MBEDGE][ee].getOrder());
1377  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
1378  3 * nb_dofs, false);
1379  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
1380  9 * nb_dofs, false);
1381  hcurl_edge_n[ee] =
1382  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
1383  diff_hcurl_edge_n[ee] =
1384  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
1385  }
1387  sense, order,
1388  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1389  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1390  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts);
1391  } else {
1392  // No DOFs on edges, resize base function matrices, indicating that no
1393  // dofs on them.
1394  for (int ee = 0; ee != 6; ee++) {
1395  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
1396  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
1397  false);
1398  }
1399  }
1400 
1401  // triangles
1402  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
1403  int order[4];
1404  // faces
1405  if (data.dataOnEntities[MBTRI].size() != 4) {
1406  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1407  "data structure for storing face h-curl base have wrong size "
1408  "should be four but is %d",
1409  data.dataOnEntities[MBTRI].size());
1410  }
1411  double *hcurl_base_n[4], *diff_hcurl_base_n[4];
1412  for (int ff = 0; ff != 4; ff++) {
1413  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
1414  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1415  "orintation of face is not set");
1416  }
1417  order[ff] = data.dataOnEntities[MBTRI][ff].getOrder();
1418  int nb_dofs = NBFACETRI_DEMKOWICZ_HCURL(order[ff]);
1419  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1420  3 * nb_dofs, false);
1421  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1422  9 * nb_dofs, false);
1423  hcurl_base_n[ff] =
1424  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1425  diff_hcurl_base_n[ff] =
1426  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1427  }
1428  if (data.facesNodes.size1() != 4) {
1429  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1430  "data inconsistency, should be four faces");
1431  }
1432  if (data.facesNodes.size2() != 3) {
1433  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1434  "data inconsistency, should be three nodes on face");
1435  }
1437  &*data.facesNodes.data().begin(), order,
1438  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1439  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1440  hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts);
1441  } else {
1442  // No DOFs on faces, resize base function matrices, indicating that no
1443  // dofs on them.
1444  for (int ff = 0; ff != 4; ff++) {
1445  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
1446  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
1447  false);
1448  }
1449  }
1450 
1451  if (data.spacesOnEntities[MBTET].test(HCURL)) {
1452  // volume
1453  int order = data.dataOnEntities[MBTET][0].getOrder();
1454  int nb_vol_dofs = NBVOLUMETET_DEMKOWICZ_HCURL(order);
1455  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
1456  3 * nb_vol_dofs, false);
1457  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
1458  9 * nb_vol_dofs, false);
1460  data.dataOnEntities[MBTET][0].getOrder(),
1461  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1462  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1463  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
1464  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
1465  nb_gauss_pts);
1466  } else {
1467  // No DOFs on faces, resize base function matrices, indicating that no
1468  // dofs on them.
1469  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
1470  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
1471  }
1472 
1474 }
1475 
1478 
1479  switch (cTx->bAse) {
1483  break;
1484  case DEMKOWICZ_JACOBI_BASE:
1486  break;
1487  default:
1488  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1489  }
1490 
1492 }
1493 
1496  boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
1498 
1499  cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
1500 
1501  int nb_gauss_pts = pts.size2();
1502  if (!nb_gauss_pts)
1504 
1505  if (pts.size1() < 3)
1506  SETERRQ(
1507  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1508  "Wrong dimension of pts, should be at least 3 rows with coordinates");
1509 
1511  const FieldApproximationBase base = cTx->bAse;
1512  EntitiesFieldData &data = cTx->dAta;
1513  if (cTx->copyNodeBase == LASTBASE) {
1514  data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4,
1515  false);
1517  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1518  &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
1519  } else {
1520  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
1521  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
1522  }
1523  if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
1524  (unsigned int)nb_gauss_pts) {
1525  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1526  "Base functions or nodes has wrong number of integration points "
1527  "for base %s",
1528  ApproximationBaseNames[base]);
1529  }
1530  data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, false);
1531  std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
1532  data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
1533  }
1534 
1535  switch (cTx->sPace) {
1536  case H1:
1537  CHKERR getValueH1(pts);
1538  break;
1539  case HDIV:
1540  CHKERR getValueHdiv(pts);
1541  break;
1542  case HCURL:
1543  CHKERR getValueHcurl(pts);
1544  break;
1545  case L2:
1546  CHKERR getValueL2(pts);
1547  break;
1548  default:
1549  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space");
1550  }
1551 
1553 }
1554 
1556  auto it = TetBaseCache::hDivBaseFaceDemkowicz.find(ptr);
1557  if (it != TetBaseCache::hDivBaseFaceDemkowicz.end()) {
1558  MOFEM_LOG_CHANNEL("WORLD");
1559  MOFEM_TAG_AND_LOG("WORLD", Sev::noisy, "TetPolynomialBase")
1560  << "Cache off hDivBaseFaceDemkowicz: " << it->second.size();
1562  return false;
1563  } else {
1564  MOFEM_LOG_CHANNEL("WORLD");
1565  MOFEM_TAG_AND_LOG("WORLD", Sev::noisy, "TetPolynomialBase")
1566  << "Cache on hDivBaseFaceDemkowicz";
1568  return true;
1569  }
1570 }
1571 
1573  auto it = TetBaseCache::hdivBaseInteriorDemkowicz.find(ptr);
1574  if (it != TetBaseCache::hdivBaseInteriorDemkowicz.end()) {
1575  MOFEM_LOG_CHANNEL("WORLD");
1576  MOFEM_TAG_AND_LOG("WORLD", Sev::noisy, "TetPolynomialBase")
1577  << "Cache off hdivBaseInteriorDemkowicz: " << it->second.size();
1579  return false;
1580  } else {
1581  MOFEM_LOG_CHANNEL("WORLD");
1582  MOFEM_TAG_AND_LOG("WORLD", Sev::noisy, "TetPolynomialBase")
1583  << "Cache on hdivBaseInteriorDemkowicz";
1585  return true;
1586  }
1587 }
1588 
1590  std::vector<void *> v) {
1591  for (auto fe_ptr : v) {
1594  }
1597  }
1598  }
1599 }
1600 
1602  std::vector<void *> v) {
1603  for (auto fe_ptr : v) {
1605  }
1607  }
1608  }
1609 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface 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
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:98
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET
MoFEMErrorCode Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_f[], double *diff_phi_v_f[], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: Hdiv.cpp:368
NBEDGE_H1
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
Definition: h1_hdiv_hcurl_l2.h:55
MoFEM::Hdiv_Demkowicz_Face_MBTET_ON_FACE
MoFEMErrorCode Hdiv_Demkowicz_Face_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb)
Definition: Hdiv.cpp:617
TetBaseCache::BaseCacheItem::nb_gauss_pts
int nb_gauss_pts
Definition: TetPolynomialBase.cpp:14
TetBaseCache::HDivBaseCacheItem::diffN
MatrixDouble diffN
Definition: TetPolynomialBase.cpp:32
LASTBASE
@ LASTBASE
Definition: definitions.h:69
TetBaseCache::HDivBaseCacheItem::n2
int n2
Definition: TetPolynomialBase.cpp:29
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
L2_Ainsworth_ShapeFunctions_MBTET
PetscErrorCode L2_Ainsworth_ShapeFunctions_MBTET(int p, double *N, double *diffN, double *L2N, double *diff_L2N, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Get base functions on tetrahedron for L2 space.
Definition: l2.c:74
HVEC0_2
@ HVEC0_2
Definition: definitions.h:198
MoFEM::TetPolynomialBase::getValueL2
MoFEMErrorCode getValueL2(MatrixDouble &pts)
Get base functions for L2 space.
Definition: TetPolynomialBase.cpp:510
NBVOLUMETET_DEMKOWICZ_HCURL
#define NBVOLUMETET_DEMKOWICZ_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:110
TetBaseCache::HDivBaseCacheItem::N
MatrixDouble N
Definition: TetPolynomialBase.cpp:31
MoFEM::EntPolynomialBaseCtx
Class used to pass element data to calculate base functions on tet,triangle,edge.
Definition: EntPolynomialBaseCtx.hpp:22
MoFEM::TetPolynomialBase::N_face_bubble
ublas::vector< MatrixDouble > N_face_bubble
Definition: TetPolynomialBase.hpp:84
TetBaseCache::hDivBaseFaceDemkowicz
static std::map< const void *, HDivBaseFaceCacheMI > hDivBaseFaceDemkowicz
Definition: TetPolynomialBase.cpp:68
MoFEM::TetPolynomialBase::getValueHcurlAinsworthBase
MoFEMErrorCode getValueHcurlAinsworthBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:1226
MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET
MoFEMErrorCode Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f_e[4][3], double *diff_phi_f_e[4][3], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base functions, Edge-based face functions by Ainsworth .
Definition: Hdiv.cpp:13
NOBASE
@ NOBASE
Definition: definitions.h:59
TetBaseCache::BaseCacheMI
boost::multi_index_container< BaseCacheItem, boost::multi_index::indexed_by< boost::multi_index::hashed_unique< composite_key< BaseCacheItem, member< BaseCacheItem, int, &BaseCacheItem::order >, member< BaseCacheItem, int, &BaseCacheItem::nb_gauss_pts > >> > > BaseCacheMI
Definition: TetPolynomialBase.cpp:47
ApproximationBaseNames
const static char *const ApproximationBaseNames[]
Definition: definitions.h:72
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
NBFACETRI_AINSWORTH_EDGE_HDIV
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:130
MoFEM::Hcurl_Demkowicz_EdgeBaseFunctions_MBTET
MoFEMErrorCode Hcurl_Demkowicz_EdgeBaseFunctions_MBTET(int *sense, int *p, double *n, double *diff_n, double *phi[], double *diff_phi[], int nb_integration_pts)
Edge based H-curl base functions on tetrahedral.
Definition: Hcurl.cpp:2068
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
TetBaseCache
Definition: TetPolynomialBase.cpp:10
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
TetBaseCache::HDivBaseCacheItem::n0
int n0
Definition: TetPolynomialBase.cpp:27
MoFEM::TetPolynomialBase::senseFaceAlpha
MatrixInt senseFaceAlpha
Definition: TetPolynomialBase.hpp:108
HVEC0_1
@ HVEC0_1
Definition: definitions.h:195
MoFEM::EntPolynomialBaseCtx::dAta
EntitiesFieldData & dAta
Definition: EntPolynomialBaseCtx.hpp:36
NBVOLUMETET_AINSWORTH_FACE_HDIV
#define NBVOLUMETET_AINSWORTH_FACE_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:134
MoFEM::Tools::diffShapeFunMBTET
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition: Tools.hpp:271
MoFEM::TetPolynomialBase::N_face_edge
ublas::matrix< MatrixDouble > N_face_edge
Definition: TetPolynomialBase.hpp:83
NBVOLUMETET_H1
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
Definition: h1_hdiv_hcurl_l2.h:75
MoFEM::Hdiv_Demkowicz_Interior_MBTET
MoFEMErrorCode Hdiv_Demkowicz_Interior_MBTET(int p, double *N, double *diffN, int p_face[], double *phi_f[4], double *diff_phi_f[4], double *phi_v, double *diff_phi_v, int gdim)
Definition: Hdiv.cpp:762
MoFEM::TetPolynomialBase::getValueH1AinsworthBase
MoFEMErrorCode getValueH1AinsworthBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:116
MoFEM::TetPolynomialBase::diffN_volume_edge
ublas::vector< MatrixDouble > diffN_volume_edge
Definition: TetPolynomialBase.hpp:91
NBVOLUMETET_L2
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
Definition: h1_hdiv_hcurl_l2.h:27
HVEC0_0
@ HVEC0_0
Definition: definitions.h:192
HVEC1_1
@ HVEC1_1
Definition: definitions.h:196
MoFEM::TetPolynomialBase::getValueL2AinsworthBase
MoFEMErrorCode getValueL2AinsworthBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:528
MoFEM::TetPolynomialBase::TetPolynomialBase
TetPolynomialBase(const void *ptr=nullptr)
Definition: TetPolynomialBase.cpp:85
MoFEM::EntitiesFieldData::spacesOnEntities
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
Definition: EntitiesFieldData.hpp:50
MoFEM::EntitiesFieldData::facesNodes
MatrixInt facesNodes
nodes on finite element faces
Definition: EntitiesFieldData.hpp:46
TetBaseCache::BaseCacheItem::N
MatrixDouble N
Definition: TetPolynomialBase.cpp:15
HVEC1
@ HVEC1
Definition: definitions.h:186
FTensor::Tensor2< double *, 3, 3 >
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
HVEC2_1
@ HVEC2_1
Definition: definitions.h:197
TetBaseCache::hdivBaseInteriorDemkowicz
static std::map< const void *, BaseCacheMI > hdivBaseInteriorDemkowicz
Definition: TetPolynomialBase.cpp:67
MoFEM::EntPolynomialBaseCtx::copyNodeBase
const FieldApproximationBase copyNodeBase
Definition: EntPolynomialBaseCtx.hpp:40
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::TetPolynomialBase::getValueHcurl
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
Get base functions for Hcurl space.
Definition: TetPolynomialBase.cpp:1476
HVEC2_2
@ HVEC2_2
Definition: definitions.h:200
HVEC1_2
@ HVEC1_2
Definition: definitions.h:199
MoFEM::BernsteinBezier::generateIndicesEdgeTet
static MoFEMErrorCode generateIndicesEdgeTet(const int N[], int *alpha[])
Definition: BernsteinBezier.cpp:138
TetBaseCache::HDivBaseCacheItem
Definition: TetPolynomialBase.cpp:19
MoFEM::BernsteinBezier::generateIndicesTriTet
static MoFEMErrorCode generateIndicesTriTet(const int N[], int *alpha[])
Definition: BernsteinBezier.cpp:174
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::TetPolynomialBase::swichCacheHdivBaseInteriorDemkowicz
static bool swichCacheHdivBaseInteriorDemkowicz(const void *ptr)
Definition: TetPolynomialBase.cpp:1572
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:93
MoFEM::TetPolynomialBase::N_volume_edge
ublas::vector< MatrixDouble > N_volume_edge
Definition: TetPolynomialBase.hpp:85
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:36
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:89
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:86
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
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:1209
MoFEM::TetPolynomialBase::getValueHdivAinsworthBase
MoFEMErrorCode getValueHdivAinsworthBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:734
MoFEM::EntPolynomialBaseCtx::bAse
const FieldApproximationBase bAse
Definition: EntPolynomialBaseCtx.hpp:38
TetBaseCache::BaseCacheItem::diffN
MatrixDouble diffN
Definition: TetPolynomialBase.cpp:16
MoFEM::TetPolynomialBase::diffN_volume_face
ublas::vector< MatrixDouble > diffN_volume_face
Definition: TetPolynomialBase.hpp:92
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
MoFEM::TetPolynomialBase::getValueHcurlDemkowiczBase
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:1344
FTensor::Index< 'i', 3 >
NBFACETRI_DEMKOWICZ_HCURL
#define NBFACETRI_DEMKOWICZ_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:109
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
MoFEM::TetPolynomialBase::getValueL2BernsteinBezierBase
MoFEMErrorCode getValueL2BernsteinBezierBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:558
convert.n
n
Definition: convert.py:82
MoFEM::TetPolynomialBase::swichCacheHDivBaseDemkowiczOff
static void swichCacheHDivBaseDemkowiczOff(std::vector< void * > v)
Definition: TetPolynomialBase.cpp:1601
MoFEM::UnknownInterface
base class for all interface classes
Definition: UnknownInterface.hpp:34
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
MoFEM::Tools::shapeFunMBTET
static MoFEMErrorCode shapeFunMBTET(double *shape, const double *ksi, const double *eta, const double *zeta, const double nb)
Calculate shape functions on tetrahedron.
Definition: Tools.hpp:738
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
MoFEM::TetPolynomialBase::swichCacheHDivBaseDemkowiczOn
static void swichCacheHDivBaseDemkowiczOn(std::vector< void * > v)
Definition: TetPolynomialBase.cpp:1589
MoFEM::TetPolynomialBase::cTx
EntPolynomialBaseCtx * cTx
Definition: TetPolynomialBase.hpp:37
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:484
NBVOLUMETET_AINSWORTH_HCURL
#define NBVOLUMETET_AINSWORTH_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:105
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
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:90
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MoFEM::TetPolynomialBase::getValueHdivDemkowiczBase
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:1074
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::Hdiv_Ainsworth_FaceBubbleShapeFunctions
MoFEMErrorCode Hdiv_Ainsworth_FaceBubbleShapeFunctions(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f[], double *diff_phi_f[], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face bubble functions by Ainsworth .
Definition: Hdiv.cpp:137
NBVOLUMETET_AINSWORTH_HDIV
#define NBVOLUMETET_AINSWORTH_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:137
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
MoFEM::TetPolynomialBase::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: TetPolynomialBase.cpp:1495
MoFEM::EntPolynomialBaseCtx::fieldName
const std::string fieldName
Definition: EntPolynomialBaseCtx.hpp:39
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
TetBaseCache::HDivBaseCacheItem::order
int order
Definition: TetPolynomialBase.cpp:21
HVEC2_0
@ HVEC2_0
Definition: definitions.h:194
MoFEM::Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET
MoFEMErrorCode Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_e[6], double *diff_phi_v_e[6], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base function, Edge-based interior (volume) functions by Ainsworth .
Definition: Hdiv.cpp:280
MoFEM::Hcurl_Demkowicz_FaceBaseFunctions_MBTET
MoFEMErrorCode Hcurl_Demkowicz_FaceBaseFunctions_MBTET(int *faces_nodes, int *p, double *n, double *diff_n, double *phi[], double *diff_phi[], int nb_integration_pts)
Face base interior function.
Definition: Hcurl.cpp:2395
MoFEM::TetPolynomialBase::N_volume_bubble
MatrixDouble N_volume_bubble
Definition: TetPolynomialBase.hpp:87
MoFEM::EntitiesFieldData::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition: EntitiesFieldData.hpp:56
NBFACETRI_AINSWORTH_FACE_HDIV
#define NBFACETRI_AINSWORTH_FACE_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:131
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:440
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
HVEC2
@ HVEC2
Definition: definitions.h:186
MoFEM::EntitiesFieldData
data structure for finite element entity
Definition: EntitiesFieldData.hpp:40
MoFEM::TetPolynomialBase::getValueH1BernsteinBezierBase
MoFEMErrorCode getValueH1BernsteinBezierBase(MatrixDouble &pts)
Definition: TetPolynomialBase.cpp:228
NBFACETRI_AINSWORTH_HDIV
#define NBFACETRI_AINSWORTH_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:132
TetBaseCache::BaseCacheItem::order
int order
Definition: TetPolynomialBase.cpp:13
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:416
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:87
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
HVEC1_0
@ HVEC1_0
Definition: definitions.h:193
MoFEM::TetPolynomialBase::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition: TetPolynomialBase.cpp:77
MoFEM::TetPolynomialBase::swichCacheHDivBaseFaceDemkowicz
static bool swichCacheHDivBaseFaceDemkowicz(const void *ptr)
Definition: TetPolynomialBase.cpp:1555
TetBaseCache::HDivBaseCacheItem::nb_gauss_pts
int nb_gauss_pts
Definition: TetPolynomialBase.cpp:22
MoFEM::TetPolynomialBase
Calculate base functions on tetrahedral.
Definition: TetPolynomialBase.hpp:19
MoFEM::EntPolynomialBaseCtx::basePolynomialsType0
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Definition: EntPolynomialBaseCtx.hpp:27