v0.14.0
TriPolynomialBase.cpp
Go to the documentation of this file.
1 /** \file TriPolynomialBase.cpp
2 \brief Implementation of bases on triangle
3 
4 */
5 
6 using namespace MoFEM;
7 
8 MoFEMErrorCode
9 TriPolynomialBase::query_interface(boost::typeindex::type_index type_index,
10  UnknownInterface **iface) const {
11  MoFEMFunctionBegin;
12  *iface = const_cast<TriPolynomialBase *>(this);
13  MoFEMFunctionReturn(0);
14 }
15 
16 MoFEMErrorCode TriPolynomialBase::getValueH1(MatrixDouble &pts) {
17  MoFEMFunctionBegin;
18 
19  switch (cTx->bAse) {
20  case AINSWORTH_LEGENDRE_BASE:
21  case AINSWORTH_LOBATTO_BASE:
22  CHKERR getValueH1AinsworthBase(pts);
23  break;
24  case AINSWORTH_BERNSTEIN_BEZIER_BASE:
25  CHKERR getValueH1BernsteinBezierBase(pts);
26  break;
27  default:
28  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
29  }
30 
31  MoFEMFunctionReturn(0);
32 }
33 
34 MoFEMErrorCode
35 TriPolynomialBase::getValueH1BernsteinBezierBase(MatrixDouble &pts) {
36  MoFEMFunctionBegin;
37  EntitiesFieldData &data = cTx->dAta;
38  const std::string &field_name = cTx->fieldName;
39  int nb_gauss_pts = pts.size2();
40 
41  auto get_alpha = [field_name](auto &data) -> MatrixInt & {
42  auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
43  if (!ptr)
44  ptr.reset(new MatrixInt());
45  return *ptr;
46  };
47 
48  auto get_base = [field_name](auto &data) -> MatrixDouble & {
49  auto &ptr = data.getBBNSharedPtr(field_name);
50  if (!ptr)
51  ptr.reset(new MatrixDouble());
52  return *ptr;
53  };
54 
55  auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
56  auto &ptr = data.getBBDiffNSharedPtr(field_name);
57  if (!ptr)
58  ptr.reset(new MatrixDouble());
59  return *ptr;
60  };
61 
62  auto get_alpha_by_name_ptr =
63  [](auto &data,
64  const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
65  return data.getBBAlphaIndicesSharedPtr(field_name);
66  };
67 
68  auto get_base_by_name_ptr =
69  [](auto &data,
70  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
71  return data.getBBNSharedPtr(field_name);
72  };
73 
74  auto get_diff_base_by_name_ptr =
75  [](auto &data,
76  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
77  return data.getBBDiffNSharedPtr(field_name);
78  };
79 
80  auto get_alpha_by_order_ptr =
81  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
82  return data.getBBAlphaIndicesByOrderSharedPtr(o);
83  };
84 
85  auto get_base_by_order_ptr =
86  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
87  return data.getBBNByOrderSharedPtr(o);
88  };
89 
90  auto get_diff_base_by_order_ptr =
91  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
92  return data.getBBDiffNByOrderSharedPtr(o);
93  };
94 
95  auto &vert_get_n = get_base(data.dataOnEntities[MBVERTEX][0]);
96  auto &vert_get_diff_n = get_diff_base(data.dataOnEntities[MBVERTEX][0]);
97  vert_get_n.resize(nb_gauss_pts, 3, false);
98  vert_get_diff_n.resize(nb_gauss_pts, 6, false);
99 
100  if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
101  (unsigned int)nb_gauss_pts)
102  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
103  "Base functions or nodes has wrong number of integration points "
104  "for base %s",
105  ApproximationBaseNames[NOBASE]);
106  auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
107 
108  auto &vertex_alpha = get_alpha(data.dataOnEntities[MBVERTEX][0]);
109  vertex_alpha.resize(3, 3, false);
110  vertex_alpha.clear();
111  for (int n = 0; n != 3; ++n)
112  vertex_alpha(n, n) = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n];
113 
114  CHKERR BernsteinBezier::baseFunctionsTri(
115  1, lambda.size1(), vertex_alpha.size1(), &vertex_alpha(0, 0),
116  &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &vert_get_n(0, 0),
117  &vert_get_diff_n(0, 0));
118  for (int n = 0; n != 3; ++n) {
119  const int f = boost::math::factorial<double>(
120  data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n]);
121  for (int g = 0; g != nb_gauss_pts; ++g) {
122  vert_get_n(g, n) *= f;
123  for (int d = 0; d != 2; ++d)
124  vert_get_diff_n(g, 2 * n + d) *= f;
125  }
126  }
127 
128  // edges
129  if (data.spacesOnEntities[MBEDGE].test(H1)) {
130  if (data.dataOnEntities[MBEDGE].size() != 3)
131  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
132  "Wrong size of ent data");
133 
134  constexpr int edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
135  for (int ee = 0; ee != 3; ++ee) {
136  auto &ent_data = data.dataOnEntities[MBEDGE][ee];
137 
138  if (ent_data.getSense() == 0)
139  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
140  "Sense of the edge unknown");
141 
142  const int sense = ent_data.getSense();
143  const int order = ent_data.getOrder();
144  const int nb_dofs = NBEDGE_H1(ent_data.getOrder());
145 
146  if (nb_dofs) {
147  if (get_alpha_by_order_ptr(ent_data, order)) {
148  get_alpha_by_name_ptr(ent_data, field_name) =
149  get_alpha_by_order_ptr(ent_data, order);
150  get_base_by_name_ptr(ent_data, field_name) =
151  get_base_by_order_ptr(ent_data, order);
152  get_diff_base_by_name_ptr(ent_data, field_name) =
153  get_diff_base_by_order_ptr(ent_data, order);
154  } else {
155  auto &get_n = get_base(ent_data);
156  auto &get_diff_n = get_diff_base(ent_data);
157  get_n.resize(nb_gauss_pts, nb_dofs, false);
158  get_diff_n.resize(nb_gauss_pts, 2 * nb_dofs, false);
159 
160  auto &edge_alpha = get_alpha(data.dataOnEntities[MBEDGE][ee]);
161  edge_alpha.resize(nb_dofs, 3, false);
162  CHKERR BernsteinBezier::generateIndicesEdgeTri(ee, order,
163  &edge_alpha(0, 0));
164  if (sense == -1) {
165  for (int i = 0; i != edge_alpha.size1(); ++i) {
166  int a = edge_alpha(i, edges_nodes[ee][0]);
167  edge_alpha(i, edges_nodes[ee][0]) =
168  edge_alpha(i, edges_nodes[ee][1]);
169  edge_alpha(i, edges_nodes[ee][1]) = a;
170  }
171  }
172  CHKERR BernsteinBezier::baseFunctionsTri(
173  order, lambda.size1(), edge_alpha.size1(), &edge_alpha(0, 0),
174  &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &get_n(0, 0),
175  &get_diff_n(0, 0));
176 
177  get_alpha_by_order_ptr(ent_data, order) =
178  get_alpha_by_name_ptr(ent_data, field_name);
179  get_base_by_order_ptr(ent_data, order) =
180  get_base_by_name_ptr(ent_data, field_name);
181  get_diff_base_by_order_ptr(ent_data, order) =
182  get_diff_base_by_name_ptr(ent_data, field_name);
183  }
184  }
185  }
186  } else {
187  for (int ee = 0; ee != 3; ++ee) {
188  auto &ent_data = data.dataOnEntities[MBEDGE][ee];
189  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
190  auto &get_n = get_base(ent_data);
191  auto &get_diff_n = get_diff_base(ent_data);
192  get_n.resize(nb_gauss_pts, 0, false);
193  get_diff_n.resize(nb_gauss_pts, 0, false);
194  }
195  }
196 
197  if (data.spacesOnEntities[MBTRI].test(H1)) {
198  if (data.dataOnEntities[MBTRI].size() != 1)
199  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
200  "Wrong size ent of ent data");
201 
202  auto &ent_data = data.dataOnEntities[MBTRI][0];
203  const int order = ent_data.getOrder();
204  const int nb_dofs = NBFACETRI_H1(order);
205  if (get_alpha_by_order_ptr(ent_data, order)) {
206  get_alpha_by_name_ptr(ent_data, field_name) =
207  get_alpha_by_order_ptr(ent_data, order);
208  get_base_by_name_ptr(ent_data, field_name) =
209  get_base_by_order_ptr(ent_data, order);
210  get_diff_base_by_name_ptr(ent_data, field_name) =
211  get_diff_base_by_order_ptr(ent_data, order);
212  } else {
213 
214  auto &get_n = get_base(ent_data);
215  auto &get_diff_n = get_diff_base(ent_data);
216  get_n.resize(nb_gauss_pts, nb_dofs, false);
217  get_diff_n.resize(nb_gauss_pts, 2 * nb_dofs, false);
218  if (nb_dofs) {
219  auto &tri_alpha = get_alpha(ent_data);
220  tri_alpha.resize(nb_dofs, 3, false);
221 
222  CHKERR BernsteinBezier::generateIndicesTriTri(order, &tri_alpha(0, 0));
223  CHKERR BernsteinBezier::baseFunctionsTri(
224  order, lambda.size1(), tri_alpha.size1(), &tri_alpha(0, 0),
225  &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &get_n(0, 0),
226  &get_diff_n(0, 0));
227 
228  get_alpha_by_order_ptr(ent_data, order) =
229  get_alpha_by_name_ptr(ent_data, field_name);
230  get_base_by_order_ptr(ent_data, order) =
231  get_base_by_name_ptr(ent_data, field_name);
232  get_diff_base_by_order_ptr(ent_data, order) =
233  get_diff_base_by_name_ptr(ent_data, field_name);
234  }
235  }
236  } else {
237  auto &ent_data = data.dataOnEntities[MBTRI][0];
238  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
239  auto &get_n = get_base(ent_data);
240  auto &get_diff_n = get_diff_base(ent_data);
241  get_n.resize(nb_gauss_pts, 0, false);
242  get_diff_n.resize(nb_gauss_pts, 0, false);
243  }
244 
245  MoFEMFunctionReturn(0);
246 }
247 
248 MoFEMErrorCode TriPolynomialBase::getValueH1AinsworthBase(MatrixDouble &pts) {
249  MoFEMFunctionBegin;
250 
251  EntitiesFieldData &data = cTx->dAta;
252  const FieldApproximationBase base = cTx->bAse;
253 
254  if (cTx->basePolynomialsType0 == NULL)
255  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
256  "Polynomial type not set");
257 
258  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
259  double *diffL, const int dim) =
260  cTx->basePolynomialsType0;
261 
262  int nb_gauss_pts = pts.size2();
263 
264  if (data.spacesOnEntities[MBEDGE].test(H1)) {
265  // edges
266  if (data.dataOnEntities[MBEDGE].size() != 3)
267  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
268  "expected size of data.dataOnEntities[MBEDGE] is 3");
269 
270  int sense[3], order[3];
271  double *H1edgeN[3], *diffH1edgeN[3];
272  for (int ee = 0; ee < 3; ee++) {
273 
274  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0)
275  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
276  "sense of the edge unknown");
277 
278  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
279  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
280  int nb_dofs = NBEDGE_H1(data.dataOnEntities[MBEDGE][ee].getOrder());
281  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
282  false);
283  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
284  2 * nb_dofs, false);
285  H1edgeN[ee] = &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
286  diffH1edgeN[ee] =
287  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
288  }
289  CHKERR H1_EdgeShapeFunctions_MBTRI(
290  sense, order,
291  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
292  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
293  H1edgeN, diffH1edgeN, nb_gauss_pts, base_polynomials);
294  }
295 
296  if (data.spacesOnEntities[MBTRI].test(H1)) {
297  // face
298  if (data.dataOnEntities[MBTRI].size() != 1) {
299  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
300  "expected that size data.dataOnEntities[MBTRI] is one");
301  }
302 
303  int nb_dofs = NBFACETRI_H1(data.dataOnEntities[MBTRI][0].getOrder());
304  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, nb_dofs,
305  false);
306  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
307  2 * nb_dofs, false);
308  const int face_nodes[] = {0, 1, 2};
309  CHKERR H1_FaceShapeFunctions_MBTRI(
310  face_nodes, data.dataOnEntities[MBTRI][0].getOrder(),
311  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
312  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
313  &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
314  &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
315  nb_gauss_pts, base_polynomials);
316  }
317 
318  MoFEMFunctionReturn(0);
319 }
320 
321 MoFEMErrorCode TriPolynomialBase::getValueL2(MatrixDouble &pts) {
322  MoFEMFunctionBegin;
323 
324  switch (cTx->bAse) {
325  case AINSWORTH_LEGENDRE_BASE:
326  case AINSWORTH_LOBATTO_BASE:
327  case DEMKOWICZ_JACOBI_BASE:
328  CHKERR getValueL2AinsworthBase(pts);
329  break;
330  case AINSWORTH_BERNSTEIN_BEZIER_BASE:
331  CHKERR getValueL2BernsteinBezierBase(pts);
332  break;
333  default:
334  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
335  }
336 
337  MoFEMFunctionReturn(0);
338 }
339 
340 MoFEMErrorCode TriPolynomialBase::getValueL2AinsworthBase(MatrixDouble &pts) {
341  MoFEMFunctionBegin;
342 
343  EntitiesFieldData &data = cTx->dAta;
344  const FieldApproximationBase base = cTx->bAse;
345  if (cTx->basePolynomialsType0 == NULL)
346  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
347  "Polynomial type not set");
348  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
349  double *diffL, const int dim) =
350  cTx->basePolynomialsType0;
351 
352  int nb_gauss_pts = pts.size2();
353 
354  data.dataOnEntities[MBTRI][0].getN(base).resize(
355  nb_gauss_pts, NBFACETRI_L2(data.dataOnEntities[MBTRI][0].getOrder()),
356  false);
357  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(
358  nb_gauss_pts, 2 * NBFACETRI_L2(data.dataOnEntities[MBTRI][0].getOrder()),
359  false);
360 
361  CHKERR L2_Ainsworth_ShapeFunctions_MBTRI(
362  data.dataOnEntities[MBTRI][0].getOrder(),
363  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
364  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
365  &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
366  &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
367  nb_gauss_pts, base_polynomials);
368 
369  MoFEMFunctionReturn(0);
370 }
371 
372 MoFEMErrorCode
373 TriPolynomialBase::getValueL2BernsteinBezierBase(MatrixDouble &pts) {
374  MoFEMFunctionBegin;
375 
376  EntitiesFieldData &data = cTx->dAta;
377  const std::string &field_name = cTx->fieldName;
378  int nb_gauss_pts = pts.size2();
379 
380  auto get_alpha = [field_name](auto &data) -> MatrixInt & {
381  auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
382  if (!ptr)
383  ptr.reset(new MatrixInt());
384  return *ptr;
385  };
386 
387  auto get_base = [field_name](auto &data) -> MatrixDouble & {
388  auto &ptr = data.getBBNSharedPtr(field_name);
389  if (!ptr)
390  ptr.reset(new MatrixDouble());
391  return *ptr;
392  };
393 
394  auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
395  auto &ptr = data.getBBDiffNSharedPtr(field_name);
396  if (!ptr)
397  ptr.reset(new MatrixDouble());
398  return *ptr;
399  };
400 
401  auto get_alpha_by_name_ptr =
402  [](auto &data,
403  const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
404  return data.getBBAlphaIndicesSharedPtr(field_name);
405  };
406 
407  auto get_base_by_name_ptr =
408  [](auto &data,
409  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
410  return data.getBBNSharedPtr(field_name);
411  };
412 
413  auto get_diff_base_by_name_ptr =
414  [](auto &data,
415  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
416  return data.getBBDiffNSharedPtr(field_name);
417  };
418 
419  auto get_alpha_by_order_ptr =
420  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
421  return data.getBBAlphaIndicesByOrderSharedPtr(o);
422  };
423 
424  auto get_base_by_order_ptr =
425  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
426  return data.getBBNByOrderSharedPtr(o);
427  };
428 
429  auto get_diff_base_by_order_ptr =
430  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
431  return data.getBBDiffNByOrderSharedPtr(o);
432  };
433 
434  if (data.spacesOnEntities[MBTRI].test(L2)) {
435  if (data.dataOnEntities[MBTRI].size() != 1)
436  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
437  "Wrong size ent of ent data");
438 
439  if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
440  (unsigned int)nb_gauss_pts)
441  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
442  "Base functions or nodes has wrong number of integration points "
443  "for base %s",
444  ApproximationBaseNames[NOBASE]);
445  auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
446 
447  auto &ent_data = data.dataOnEntities[MBTRI][0];
448  const int order = ent_data.getOrder();
449  const int nb_dofs = NBFACETRI_L2(order);
450 
451  if (get_alpha_by_order_ptr(ent_data, order)) {
452  get_alpha_by_name_ptr(ent_data, field_name) =
453  get_alpha_by_order_ptr(ent_data, order);
454  get_base_by_name_ptr(ent_data, field_name) =
455  get_base_by_order_ptr(ent_data, order);
456  get_diff_base_by_name_ptr(ent_data, field_name) =
457  get_diff_base_by_order_ptr(ent_data, order);
458  } else {
459 
460  auto &get_n = get_base(ent_data);
461  auto &get_diff_n = get_diff_base(ent_data);
462  get_n.resize(nb_gauss_pts, nb_dofs, false);
463  get_diff_n.resize(nb_gauss_pts, 2 * nb_dofs, false);
464  if (nb_dofs) {
465 
466  if (order == 0) {
467 
468  if (nb_dofs != 1)
469  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
470  "Inconsistent number of DOFs");
471 
472  auto &tri_alpha = get_alpha(ent_data);
473  tri_alpha.clear();
474  get_n(0, 0) = 1;
475  get_diff_n.clear();
476 
477  } else {
478 
479  if (nb_dofs != 3 + 3 * NBEDGE_H1(order) + NBFACETRI_H1(order))
480  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
481  "Inconsistent number of DOFs");
482 
483  auto &tri_alpha = get_alpha(ent_data);
484  tri_alpha.resize(nb_dofs, 3, false);
485 
486  CHKERR BernsteinBezier::generateIndicesVertexTri(order,
487  &tri_alpha(0, 0));
488 
489  if (order > 1) {
490  std::array<int, 3> face_n{order, order, order};
491  std::array<int *, 3> face_edge_ptr{
492  &tri_alpha(3, 0), &tri_alpha(3 + NBEDGE_H1(order), 0),
493  &tri_alpha(3 + 2 * NBEDGE_H1(order), 0)};
494  CHKERR BernsteinBezier::generateIndicesEdgeTri(
495  face_n.data(), face_edge_ptr.data());
496  if (order > 2)
497  CHKERR BernsteinBezier::generateIndicesTriTri(
498  order, &tri_alpha(3 + 3 * NBEDGE_H1(order), 0));
499  }
500  CHKERR BernsteinBezier::baseFunctionsTri(
501  order, lambda.size1(), tri_alpha.size1(), &tri_alpha(0, 0),
502  &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &get_n(0, 0),
503  &get_diff_n(0, 0));
504  }
505 
506  get_alpha_by_order_ptr(ent_data, order) =
507  get_alpha_by_name_ptr(ent_data, field_name);
508  get_base_by_order_ptr(ent_data, order) =
509  get_base_by_name_ptr(ent_data, field_name);
510  get_diff_base_by_order_ptr(ent_data, order) =
511  get_diff_base_by_name_ptr(ent_data, field_name);
512  }
513  }
514  } else {
515  auto &ent_data = data.dataOnEntities[MBTRI][0];
516  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
517  auto &get_n = get_base(ent_data);
518  auto &get_diff_n = get_diff_base(ent_data);
519  get_n.resize(nb_gauss_pts, 0, false);
520  get_diff_n.resize(nb_gauss_pts, 0, false);
521  }
522 
523  MoFEMFunctionReturn(0);
524 }
525 
526 MoFEMErrorCode TriPolynomialBase::getValueHdivAinsworthBase(MatrixDouble &pts) {
527  MoFEMFunctionBegin;
528 
529  EntitiesFieldData &data = cTx->dAta;
530  const FieldApproximationBase base = cTx->bAse;
531  if (cTx->basePolynomialsType0 == NULL)
532  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
533  "Polynomial type not set");
534  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
535  double *diffL, const int dim) =
536  cTx->basePolynomialsType0;
537 
538  int nb_gauss_pts = pts.size2();
539 
540  double *phi_f_e[3];
541  double *phi_f;
542 
543  N_face_edge.resize(1, 3, false);
544  N_face_bubble.resize(1, false);
545  int face_order = data.dataOnEntities[MBTRI][0].getOrder();
546  if (face_order > 0) {
547  // three edges on face
548  for (int ee = 0; ee < 3; ee++) {
549  N_face_edge(0, ee).resize(
550  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(face_order), false);
551  phi_f_e[ee] = &((N_face_edge(0, ee))(0, 0));
552  }
553  N_face_bubble[0].resize(
554  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_FACE_HDIV(face_order), false);
555  phi_f = &*(N_face_bubble[0].data().begin());
556 
557  int face_nodes[3] = {0, 1, 2};
558  CHKERR Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(
559  face_nodes, face_order,
560  &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0), NULL, phi_f_e, NULL,
561  nb_gauss_pts, 3, base_polynomials);
562  CHKERR Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(
563  face_nodes, face_order,
564  &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0), NULL, phi_f, NULL,
565  nb_gauss_pts, 3, base_polynomials);
566 
567  // set shape functions into data structure
568  if (data.dataOnEntities[MBTRI].size() != 1) {
569  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
570  }
571  data.dataOnEntities[MBTRI][0].getN(base).resize(
572  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_HDIV(face_order), false);
573  int col = 0;
574  for (int oo = 0; oo < face_order; oo++) {
575  for (int ee = 0; ee < 3; ee++) {
576  for (int dd = 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
577  dd < 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++, col++) {
578  for (int gg = 0; gg < nb_gauss_pts; gg++) {
579  data.dataOnEntities[MBTRI][0].getN(base)(gg, col) =
580  N_face_edge(0, ee)(gg, dd);
581  }
582  }
583  }
584  for (int dd = 3 * NBFACETRI_AINSWORTH_FACE_HDIV(oo);
585  dd < 3 * NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++, col++) {
586  for (int gg = 0; gg < nb_gauss_pts; gg++) {
587  data.dataOnEntities[MBTRI][0].getN(base)(gg, col) =
588  N_face_bubble[0](gg, dd);
589  }
590  }
591  }
592  }
593 
594  MoFEMFunctionReturn(0);
595 }
596 
597 MoFEMErrorCode TriPolynomialBase::getValueHdivDemkowiczBase(MatrixDouble &pts) {
598  MoFEMFunctionBegin;
599 
600  EntitiesFieldData &data = cTx->dAta;
601  // set shape functions into data structure
602  if (data.dataOnEntities[MBTRI].size() != 1) {
603  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
604  }
605  const FieldApproximationBase base = cTx->bAse;
606  if (base != DEMKOWICZ_JACOBI_BASE) {
607  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
608  "This should be used only with DEMKOWICZ_JACOBI_BASE "
609  "but base is %s",
610  ApproximationBaseNames[base]);
611  }
612  int order = data.dataOnEntities[MBTRI][0].getOrder();
613  int nb_gauss_pts = pts.size2();
614  data.dataOnEntities[MBTRI][0].getN(base).resize(
615  nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
616  double *phi_f = &*data.dataOnEntities[MBTRI][0].getN(base).data().begin();
617  if (NBFACETRI_DEMKOWICZ_HDIV(order) == 0)
618  MoFEMFunctionReturnHot(0);
619  int face_nodes[3] = {0, 1, 2};
620  CHKERR Hdiv_Demkowicz_Face_MBTET_ON_FACE(
621  face_nodes, order,
622  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
623  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(), phi_f,
624  NULL, nb_gauss_pts, 3);
625 
626  MoFEMFunctionReturn(0);
627 }
628 
629 MoFEMErrorCode TriPolynomialBase::getValueHdiv(MatrixDouble &pts) {
630  MoFEMFunctionBegin;
631 
632  switch (cTx->bAse) {
633  case AINSWORTH_LEGENDRE_BASE:
634  case AINSWORTH_LOBATTO_BASE:
635  return getValueHdivAinsworthBase(pts);
636  case DEMKOWICZ_JACOBI_BASE:
637  return getValueHdivDemkowiczBase(pts);
638  default:
639  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
640  }
641 
642  MoFEMFunctionReturn(0);
643 }
644 
645 MoFEMErrorCode
646 TriPolynomialBase::getValueHcurlAinsworthBase(MatrixDouble &pts) {
647  MoFEMFunctionBegin;
648 
649  EntitiesFieldData &data = cTx->dAta;
650  const FieldApproximationBase base = cTx->bAse;
651  if (data.dataOnEntities[MBTRI].size() != 1)
652  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
653 
654  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
655  double *diffL, const int dim) =
656  cTx->basePolynomialsType0;
657 
658  int nb_gauss_pts = pts.size2();
659 
660  // Calculation H-curl on triangle faces
661  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
662  if (data.dataOnEntities[MBEDGE].size() != 3) {
663  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
664  }
665  int sense[3], order[3];
666  double *hcurl_edge_n[3];
667  double *diff_hcurl_edge_n[3];
668  for (int ee = 0; ee < 3; ee++) {
669  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
670  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
671  "data inconsistency");
672  }
673  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
674  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
675  int nb_dofs =
676  NBEDGE_AINSWORTH_HCURL(data.dataOnEntities[MBEDGE][ee].getOrder());
677  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
678  3 * nb_dofs, false);
679  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(
680  nb_gauss_pts, 2 * 3 * nb_dofs, false);
681  hcurl_edge_n[ee] =
682  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
683  diff_hcurl_edge_n[ee] =
684  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
685  }
686  CHKERR Hcurl_Ainsworth_EdgeBaseFunctions_MBTET_ON_FACE(
687  sense, order,
688  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
689  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
690  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
691  } else {
692  for (int ee = 0; ee < 3; ee++) {
693  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
694  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
695  false);
696  }
697  }
698 
699  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
700 
701  // cerr << data.dataOnEntities[MBVERTEX][0].getN(base) << endl;
702  // cerr << data.dataOnEntities[MBVERTEX][0].getDiffN(base) << endl;
703  //
704  // face
705  if (data.dataOnEntities[MBTRI].size() != 1) {
706  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
707  }
708  int order = data.dataOnEntities[MBTRI][0].getOrder();
709  int nb_dofs = NBFACETRI_AINSWORTH_HCURL(order);
710  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
711  false);
712  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
713  3 * 2 * nb_dofs, false);
714  // cerr << data.dataOnEntities[MBVERTEX][0].getDiffN(base) << endl;
715  int face_nodes[] = {0, 1, 2};
716  CHKERR Hcurl_Ainsworth_FaceFunctions_MBTET_ON_FACE(
717  face_nodes, order,
718  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
719  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
720  &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
721  &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
722  nb_gauss_pts, base_polynomials);
723  // cerr << data.dataOnEntities[MBTRI][0].getN(base) << endl;
724  } else {
725  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 0, false);
726  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
727  }
728 
729  MoFEMFunctionReturn(0);
730 }
731 
732 MoFEMErrorCode
733 TriPolynomialBase::getValueHcurlAinsworthBrokenBase(MatrixDouble &pts) {
734  MoFEMFunctionBegin;
735 
736  EntitiesFieldData &data = cTx->dAta;
737  const FieldApproximationBase base = cTx->bAse;
738  if (data.dataOnEntities[MBTRI].size() != 1)
739  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
740 
741  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
742  double *diffL, const int dim) =
743  cTx->basePolynomialsType0;
744 
745  int nb_gauss_pts = pts.size2();
746 
747  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
748 
749  // face
750  if (data.dataOnEntities[MBTRI].size() != 1) {
751  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
752  }
753 
754  int order = data.dataOnEntities[MBTRI][0].getOrder();
755  auto nb_edge_dofs = NBEDGE_AINSWORTH_HCURL(order);
756  int nb_dofs_face = NBFACETRI_AINSWORTH_HCURL(order);
757 
758  // three faces on triangle
759  int nb_dofs = 3 * nb_edge_dofs + nb_dofs_face;
760  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
761  false);
762  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
763  3 * 2 * nb_dofs, false);
764 
765  MatrixDouble edge_bases[] = {
766  MatrixDouble(nb_gauss_pts, 3 * nb_edge_dofs, false),
767  MatrixDouble(nb_gauss_pts, 3 * nb_edge_dofs, false),
768  MatrixDouble(nb_gauss_pts, 3 * nb_edge_dofs, false)};
769  MatrixDouble diff_edge_bases[] = {
770  MatrixDouble(nb_gauss_pts, 3 * 2 * nb_edge_dofs, false),
771  MatrixDouble(nb_gauss_pts, 3 * 2 * nb_edge_dofs, false),
772  MatrixDouble(nb_gauss_pts, 3 * 2 * nb_edge_dofs, false)};
773 
774  std::array<double *, 3> phi{&*edge_bases[0].data().begin(),
775  &*edge_bases[1].data().begin(),
776  &*edge_bases[2].data().begin()};
777  std::array<double *, 3> diff_phi{&*diff_edge_bases[0].data().begin(),
778  &*diff_edge_bases[1].data().begin(),
779  &*diff_edge_bases[2].data().begin()};
780 
781  std::array<int, 3> edge_order{order, order, order};
782  std::array<int, 3> edge_sense{1, 1, 1};
783  CHKERR Hcurl_Ainsworth_EdgeBaseFunctions_MBTET_ON_FACE(
784  edge_sense.data(), edge_order.data(),
785  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
786  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
787  phi.data(), diff_phi.data(), nb_gauss_pts, base_polynomials);
788 
789  MatrixDouble face_bases(nb_gauss_pts, 3 * nb_dofs_face, false);
790  MatrixDouble diff_face_bases(nb_gauss_pts, 3 * 2 * nb_dofs_face, false);
791  int face_nodes[] = {0, 1, 2};
792  CHKERR Hcurl_Ainsworth_FaceFunctions_MBTET_ON_FACE(
793  face_nodes, order,
794  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
795  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
796  &*face_bases.data().begin(), &*diff_face_bases.data().begin(),
797  nb_gauss_pts, base_polynomials);
798 
799  // edges
800  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_edge_base[] = {
801  getFTensor1FromPtr<3>(&*edge_bases[0].data().begin()),
802  getFTensor1FromPtr<3>(&*edge_bases[1].data().begin()),
803  getFTensor1FromPtr<3>(&*edge_bases[2].data().begin())};
804  FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> t_edge_diff_base[] = {
805  getFTensor2HVecFromPtr<3, 2>(&*diff_edge_bases[0].data().begin()),
806  getFTensor2HVecFromPtr<3, 2>(&*diff_edge_bases[1].data().begin()),
807  getFTensor2HVecFromPtr<3, 2>(&*diff_edge_bases[2].data().begin())};
808 
809  // face
810  auto t_vol_base = getFTensor1FromPtr<3>(&*face_bases.data().begin());
811  auto t_vol_diff_base =
812  getFTensor2HVecFromPtr<3, 2>(&*diff_face_bases.data().begin());
813 
814  auto t_base = getFTensor1FromPtr<3>(
815  &*data.dataOnEntities[MBTRI][0].getN(base).data().begin());
816  auto t_diff_base = getFTensor2HVecFromPtr<3, 2>(
817  &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin());
818 
819  i_FTIndex<3> i;
820  j_FTIndex<2> j;
821 
822  for (int gg = 0; gg != nb_gauss_pts; gg++) {
823  for (int oo = 0; oo < order; oo++) {
824  // edges
825  for (int dd = NBEDGE_AINSWORTH_HCURL(oo);
826  dd != NBEDGE_AINSWORTH_HCURL(oo + 1); ++dd) {
827  for (int ee = 0; ee != 3; ++ee) {
828  t_base(i) = t_edge_base[ee](i);
829  ++t_base;
830  ++t_edge_base[ee];
831  t_diff_base(i, j) = t_edge_diff_base[ee](i, j);
832  ++t_diff_base;
833  ++t_edge_diff_base[ee];
834  }
835  }
836  // face
837  for (int dd = NBFACETRI_AINSWORTH_HCURL(oo);
838  dd != NBFACETRI_AINSWORTH_HCURL(oo + 1); ++dd) {
839  t_base(i) = t_vol_base(i);
840  ++t_base;
841  ++t_vol_base;
842  t_diff_base(i, j) = t_vol_diff_base(i, j);
843  ++t_diff_base;
844  ++t_vol_diff_base;
845  }
846  }
847  }
848 
849  } else {
850  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 0, false);
851  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
852  }
853 
854  MoFEMFunctionReturn(0);
855 }
856 
857 MoFEMErrorCode
858 TriPolynomialBase::getValueHcurlDemkowiczBase(MatrixDouble &pts) {
859  MoFEMFunctionBegin;
860 
861  EntitiesFieldData &data = cTx->dAta;
862  const FieldApproximationBase base = cTx->bAse;
863 
864  int nb_gauss_pts = pts.size2();
865 
866  // Calculation H-curl on triangle faces
867  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
868 
869  if (data.dataOnEntities[MBEDGE].size() != 3)
870  SETERRQ1(
871  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
872  "wrong number of data structures on edges, should be three but is %d",
873  data.dataOnEntities[MBEDGE].size());
874 
875  int sense[3], order[3];
876  double *hcurl_edge_n[3];
877  double *diff_hcurl_edge_n[3];
878 
879  for (int ee = 0; ee != 3; ++ee) {
880 
881  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0)
882  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
883  "orientation (sense) of edge is not set");
884 
885  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
886  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
887  int nb_dofs =
888  NBEDGE_DEMKOWICZ_HCURL(data.dataOnEntities[MBEDGE][ee].getOrder());
889  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
890  3 * nb_dofs, false);
891  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(
892  nb_gauss_pts, 2 * 3 * nb_dofs, false);
893 
894  hcurl_edge_n[ee] =
895  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
896  diff_hcurl_edge_n[ee] =
897  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
898  }
899 
900  CHKERR Hcurl_Demkowicz_EdgeBaseFunctions_MBTRI(
901  sense, order,
902  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
903  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
904  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts);
905 
906  } else {
907 
908  // No DOFs on faces, resize base function matrices, indicating that no
909  // dofs on them.
910  for (int ee = 0; ee != 3; ++ee) {
911  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
912  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
913  false);
914  }
915  }
916 
917  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
918 
919  // face
920  if (data.dataOnEntities[MBTRI].size() != 1)
921  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
922  "No data struture to keep base functions on face");
923 
924  int order = data.dataOnEntities[MBTRI][0].getOrder();
925  int nb_dofs = NBFACETRI_DEMKOWICZ_HCURL(order);
926  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
927  false);
928  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
929  3 * 2 * nb_dofs, false);
930 
931  int face_nodes[] = {0, 1, 2};
932  CHKERR Hcurl_Demkowicz_FaceBaseFunctions_MBTRI(
933  face_nodes, order,
934  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
935  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
936  &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
937  &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
938  nb_gauss_pts);
939 
940  } else {
941 
942  // No DOFs on faces, resize base function matrices, indicating that no
943  // dofs on them.
944  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 0, false);
945  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
946  }
947 
948  MoFEMFunctionReturn(0);
949 }
950 
951 MoFEMErrorCode
952 TriPolynomialBase::getValueHcurlDemkowiczBrokenBase(MatrixDouble &pts) {
953  MoFEMFunctionBegin;
954 
955  EntitiesFieldData &data = cTx->dAta;
956  const FieldApproximationBase base = cTx->bAse;
957 
958  int nb_gauss_pts = pts.size2();
959 
960  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
961 
962  // face
963  if (data.dataOnEntities[MBTRI].size() != 1)
964  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
965  "No data struture to keep base functions on face");
966 
967  int order = data.dataOnEntities[MBTRI][0].getOrder();
968  int nb_edge_dofs = NBEDGE_DEMKOWICZ_HCURL(order);
969  int nb_volume_dofs = NBFACETRI_DEMKOWICZ_HCURL(order);
970  int nb_dofs = 3 * nb_edge_dofs + nb_volume_dofs;
971  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
972  false);
973  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
974  3 * 2 * nb_dofs, false);
975 
976  MatrixDouble edge_bases[] = {
977  MatrixDouble(nb_gauss_pts, 3 * nb_edge_dofs, false),
978  MatrixDouble(nb_gauss_pts, 3 * nb_edge_dofs, false),
979  MatrixDouble(nb_gauss_pts, 3 * nb_edge_dofs, false)};
980  MatrixDouble diff_edge_bases[] = {
981  MatrixDouble(nb_gauss_pts, 3 * 2 * nb_edge_dofs, false),
982  MatrixDouble(nb_gauss_pts, 3 * 2 * nb_edge_dofs, false),
983  MatrixDouble(nb_gauss_pts, 3 * 2 * nb_edge_dofs, false)};
984 
985  std::array<int, 3> edge_order{order, order, order};
986  std::array<int, 3> edge_sense{1, 1, 1};
987  std::array<double *, 3> phi{&*edge_bases[0].data().begin(),
988  &*edge_bases[1].data().begin(),
989  &*edge_bases[2].data().begin()};
990  std::array<double *, 3> diff_phi{&*diff_edge_bases[0].data().begin(),
991  &*diff_edge_bases[1].data().begin(),
992  &*diff_edge_bases[2].data().begin()};
993  CHKERR Hcurl_Demkowicz_EdgeBaseFunctions_MBTRI(
994  edge_sense.data(), edge_order.data(),
995  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
996  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
997  phi.data(), diff_phi.data(), nb_gauss_pts);
998 
999  MatrixDouble face_bases(nb_gauss_pts, 3 * nb_volume_dofs, false);
1000  MatrixDouble diff_face_bases(nb_gauss_pts, 3 * 2 * nb_volume_dofs, false);
1001  int face_nodes[] = {0, 1, 2};
1002  CHKERR Hcurl_Demkowicz_FaceBaseFunctions_MBTRI(
1003  face_nodes, order,
1004  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1005  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1006  &*face_bases.data().begin(), &*diff_face_bases.data().begin(),
1007  nb_gauss_pts);
1008 
1009  // edges
1010  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_edge_base[] = {
1011  getFTensor1FromPtr<3>(&*edge_bases[0].data().begin()),
1012  getFTensor1FromPtr<3>(&*edge_bases[1].data().begin()),
1013  getFTensor1FromPtr<3>(&*edge_bases[2].data().begin())};
1014  FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> t_edge_diff_base[] = {
1015  getFTensor2HVecFromPtr<3, 2>(&*diff_edge_bases[0].data().begin()),
1016  getFTensor2HVecFromPtr<3, 2>(&*diff_edge_bases[1].data().begin()),
1017  getFTensor2HVecFromPtr<3, 2>(&*diff_edge_bases[2].data().begin())};
1018  // face
1019  auto t_vol_base = getFTensor1FromPtr<3>(&*face_bases.data().begin());
1020  auto t_vol_diff_base =
1021  getFTensor2HVecFromPtr<3, 2>(&*diff_face_bases.data().begin());
1022 
1023  auto t_base = getFTensor1FromPtr<3>(
1024  &*data.dataOnEntities[MBTRI][0].getN(base).data().begin());
1025  auto t_diff_base = getFTensor2HVecFromPtr<3, 2>(
1026  &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin());
1027 
1028  i_FTIndex<3> i;
1029  j_FTIndex<2> j;
1030 
1031  for (int gg = 0; gg != nb_gauss_pts; gg++) {
1032  for (int oo = 0; oo < order; oo++) {
1033  // edges
1034  for (int dd = NBEDGE_DEMKOWICZ_HCURL(oo);
1035  dd != NBEDGE_DEMKOWICZ_HCURL(oo + 1); ++dd) {
1036  for (int ee = 0; ee != 3; ++ee) {
1037  t_base(i) = t_edge_base[ee](i);
1038  ++t_base;
1039  ++t_edge_base[ee];
1040  t_diff_base(i, j) = t_edge_diff_base[ee](i, j);
1041  ++t_diff_base;
1042  ++t_edge_diff_base[ee];
1043  }
1044  }
1045  // faces
1046  for (int dd = NBFACETRI_DEMKOWICZ_HCURL(oo);
1047  dd != NBFACETRI_DEMKOWICZ_HCURL(oo + 1); ++dd) {
1048  t_base(i) = t_vol_base(i);
1049  ++t_base;
1050  ++t_vol_base;
1051  t_diff_base(i, j) = t_vol_diff_base(i, j);
1052  ++t_diff_base;
1053  ++t_vol_diff_base;
1054  }
1055  }
1056  }
1057 
1058  } else {
1059 
1060  // No DOFs on faces, resize base function matrices, indicating that no
1061  // dofs on them.
1062  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 0, false);
1063  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
1064  }
1065 
1066  MoFEMFunctionReturn(0);
1067 }
1068 
1069 MoFEMErrorCode TriPolynomialBase::getValueHcurl(MatrixDouble &pts) {
1070  MoFEMFunctionBegin;
1071  switch (cTx->spaceContinuity) {
1072  case CONTINUOUS:
1073 
1074  switch (cTx->bAse) {
1075  case AINSWORTH_LEGENDRE_BASE:
1076  case AINSWORTH_LOBATTO_BASE:
1077  CHKERR getValueHcurlAinsworthBase(pts);
1078  break;
1079  case DEMKOWICZ_JACOBI_BASE:
1080  CHKERR getValueHcurlDemkowiczBase(pts);
1081  break;
1082  default:
1083  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1084  }
1085 
1086  break;
1087 
1088  case DISCONTINUOUS:
1089  switch (cTx->bAse) {
1090  case AINSWORTH_LEGENDRE_BASE:
1091  case AINSWORTH_LOBATTO_BASE:
1092  return getValueHcurlAinsworthBrokenBase(pts);
1093  case DEMKOWICZ_JACOBI_BASE:
1094  return getValueHcurlDemkowiczBrokenBase(pts);
1095  default:
1096  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1097  }
1098  break;
1099 
1100  default:
1101  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown continuity");
1102  }
1103 
1104  MoFEMFunctionReturn(0);
1105 }
1106 
1107 MoFEMErrorCode
1108 TriPolynomialBase::getValue(MatrixDouble &pts,
1109  boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
1110  MoFEMFunctionBegin;
1111 
1112  cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
1113 
1114  int nb_gauss_pts = pts.size2();
1115  if (!nb_gauss_pts) {
1116  MoFEMFunctionReturnHot(0);
1117  }
1118 
1119  if (pts.size1() < 2)
1120  SETERRQ(
1121  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1122  "Wrong dimension of pts, should be at least 3 rows with coordinates");
1123 
1124  const FieldApproximationBase base = cTx->bAse;
1125  EntitiesFieldData &data = cTx->dAta;
1126 
1127  if (base != AINSWORTH_BERNSTEIN_BEZIER_BASE) {
1128  if (cTx->copyNodeBase == LASTBASE) {
1129  data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 3,
1130  false);
1131  CHKERR ShapeMBTRI(
1132  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1133  &pts(0, 0), &pts(1, 0), nb_gauss_pts);
1134  data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(3, 2, false);
1135  CHKERR ShapeDiffMBTRI(
1136  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
1137  } else {
1138  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
1139  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
1140  data.dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(base) =
1141  data.dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(cTx->copyNodeBase);
1142  }
1143  if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
1144  (unsigned int)nb_gauss_pts) {
1145  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1146  "Base functions or nodes has wrong number of integration points "
1147  "for base %s",
1148  ApproximationBaseNames[base]);
1149  }
1150  }
1151  auto set_node_derivative_for_all_gauss_pts = [&]() {
1152  MoFEMFunctionBegin;
1153  // In linear geometry derivatives are constant,
1154  // this in expense of efficiency makes implementation
1155  // consistent between vertices and other types of entities
1156  data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(nb_gauss_pts, 6,
1157  false);
1158  for (int gg = 0; gg != nb_gauss_pts; ++gg)
1159  std::copy(Tools::diffShapeFunMBTRI.begin(),
1160  Tools::diffShapeFunMBTRI.end(),
1161  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(gg, 0));
1162  MoFEMFunctionReturn(0);
1163  };
1164 
1165  CHKERR set_node_derivative_for_all_gauss_pts();
1166 
1167  switch (cTx->spaceContinuity) {
1168  case CONTINUOUS:
1169 
1170  switch (cTx->sPace) {
1171  case H1:
1172  CHKERR getValueH1(pts);
1173  break;
1174  case HDIV:
1175  CHKERR getValueHdiv(pts);
1176  break;
1177  case HCURL:
1178  CHKERR getValueHcurl(pts);
1179  break;
1180  case L2:
1181  CHKERR getValueL2(pts);
1182  break;
1183  default:
1184  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
1185  }
1186 
1187  break;
1188 
1189  case DISCONTINUOUS:
1190  switch (cTx->sPace) {
1191  case HCURL:
1192  CHKERR getValueHcurl(pts);
1193  break;
1194  default:
1195  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
1196  }
1197  break;
1198 
1199  default:
1200  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown continuity");
1201  }
1202 
1203  MoFEMFunctionReturn(0);
1204 }
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
g
constexpr double g
Definition: shallow_wave.cpp:63
MoFEM::TriPolynomialBase::cTx
EntPolynomialBaseCtx * cTx
Definition: TriPolynomialBase.hpp:28
H1
@ H1
continuous field
Definition: definitions.h:85
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
LASTBASE
@ LASTBASE
Definition: definitions.h:69
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MoFEM::TriPolynomialBase::getValueHdivDemkowiczBase
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)
Definition: TriPolynomialBase.cpp:597
L2_Ainsworth_ShapeFunctions_MBTRI
PetscErrorCode L2_Ainsworth_ShapeFunctions_MBTRI(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 triangle for L2 space.
Definition: l2.c:19
MoFEM::EntPolynomialBaseCtx
Class used to pass element data to calculate base functions on tet,triangle,edge.
Definition: EntPolynomialBaseCtx.hpp:22
H1_EdgeShapeFunctions_MBTRI
PetscErrorCode H1_EdgeShapeFunctions_MBTRI(int *sense, int *p, double *N, double *diffN, double *edgeN[3], double *diff_edgeN[3], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
H1_EdgeShapeFunctions_MBTRI.
Definition: h1.c:17
NOBASE
@ NOBASE
Definition: definitions.h:59
MoFEM::TriPolynomialBase::N_face_edge
ublas::matrix< MatrixDouble > N_face_edge
Definition: TriPolynomialBase.hpp:38
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:37
phi
static double phi
Definition: poisson_2d_dis_galerkin.cpp:29
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::EntPolynomialBaseCtx::dAta
EntitiesFieldData & dAta
Definition: EntPolynomialBaseCtx.hpp:36
MoFEM::TriPolynomialBase::getValueHdiv
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
Definition: TriPolynomialBase.cpp:629
MoFEM::BernsteinBezier::generateIndicesTriTri
static MoFEMErrorCode generateIndicesTriTri(const int N, int *alpha)
Definition: BernsteinBezier.cpp:124
MoFEM::Hcurl_Ainsworth_FaceFunctions_MBTET_ON_FACE
MoFEMErrorCode Hcurl_Ainsworth_FaceFunctions_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, 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:1249
MoFEM::EntitiesFieldData::spacesOnEntities
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
Definition: EntitiesFieldData.hpp:49
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
order
constexpr int order
Definition: dg_projection.cpp:18
ShapeDiffMBTRI
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:194
MoFEM::Hcurl_Demkowicz_EdgeBaseFunctions_MBTRI
MoFEMErrorCode Hcurl_Demkowicz_EdgeBaseFunctions_MBTRI(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 teriangle.
Definition: Hcurl.cpp:2105
MoFEM::EntPolynomialBaseCtx::copyNodeBase
const FieldApproximationBase copyNodeBase
Definition: EntPolynomialBaseCtx.hpp:41
MoFEM::TriPolynomialBase::getValueHdivAinsworthBase
MoFEMErrorCode getValueHdivAinsworthBase(MatrixDouble &pts)
Definition: TriPolynomialBase.cpp:526
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::Hcurl_Ainsworth_EdgeBaseFunctions_MBTET_ON_FACE
MoFEMErrorCode Hcurl_Ainsworth_EdgeBaseFunctions_MBTET_ON_FACE(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 face.
Definition: Hcurl.cpp:237
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
H1_FaceShapeFunctions_MBTRI
PetscErrorCode H1_FaceShapeFunctions_MBTRI(const int *face_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:191
a
constexpr double a
Definition: approx_sphere.cpp:30
MoFEM::TriPolynomialBase::getValueH1
MoFEMErrorCode getValueH1(MatrixDouble &pts)
Definition: TriPolynomialBase.cpp:16
MoFEM::TriPolynomialBase::getValueHcurlDemkowiczBase
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)
Definition: TriPolynomialBase.cpp:858
MoFEM::TriPolynomialBase
Calculate base functions on triangle.
Definition: TriPolynomialBase.hpp:16
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::TriPolynomialBase::getValueL2BernsteinBezierBase
MoFEMErrorCode getValueL2BernsteinBezierBase(MatrixDouble &pts)
Definition: TriPolynomialBase.cpp:373
MoFEM::BernsteinBezier::generateIndicesEdgeTri
static MoFEMErrorCode generateIndicesEdgeTri(const int N[], int *alpha[])
Definition: BernsteinBezier.cpp:99
MoFEM::EntPolynomialBaseCtx::sPace
const FieldSpace sPace
Definition: EntPolynomialBaseCtx.hpp:37
MoFEM::TriPolynomialBase::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: TriPolynomialBase.cpp:1108
MoFEM::TriPolynomialBase::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition: TriPolynomialBase.cpp:9
NBEDGE_DEMKOWICZ_HCURL
#define NBEDGE_DEMKOWICZ_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:108
MoFEM::TriPolynomialBase::getValueHcurlAinsworthBrokenBase
MoFEMErrorCode getValueHcurlAinsworthBrokenBase(MatrixDouble &pts)
Definition: TriPolynomialBase.cpp:733
MoFEM::L
VectorDouble L
Definition: Projection10NodeCoordsOnField.cpp:124
MoFEM::TriPolynomialBase::getValueL2AinsworthBase
MoFEMErrorCode getValueL2AinsworthBase(MatrixDouble &pts)
Definition: TriPolynomialBase.cpp:340
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
MoFEM::Tools::diffShapeFunMBTRI
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:104
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:160
MoFEM::TriPolynomialBase::getValueHcurl
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
Definition: TriPolynomialBase.cpp:1069
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::EntPolynomialBaseCtx::bAse
const FieldApproximationBase bAse
Definition: EntPolynomialBaseCtx.hpp:39
MoFEM::TriPolynomialBase::getValueH1AinsworthBase
MoFEMErrorCode getValueH1AinsworthBase(MatrixDouble &pts)
Definition: TriPolynomialBase.cpp:248
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
FTensor::Index< 'i', DIM >
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
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
MoFEM::TriPolynomialBase::getValueHcurlDemkowiczBrokenBase
MoFEMErrorCode getValueHcurlDemkowiczBrokenBase(MatrixDouble &pts)
Definition: TriPolynomialBase.cpp:952
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
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::TriPolynomialBase::getValueH1BernsteinBezierBase
MoFEMErrorCode getValueH1BernsteinBezierBase(MatrixDouble &pts)
Definition: TriPolynomialBase.cpp:35
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
MoFEM::getFTensor2HVecFromPtr< 3, 2 >
FTensor::Tensor2< FTensor::PackPtr< double *, 6 >, 3, 2 > getFTensor2HVecFromPtr< 3, 2 >(double *ptr)
Definition: Templates.hpp:935
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
NBFACETRI_L2
#define NBFACETRI_L2(P)
Number of base functions on triangle for L2 space.
Definition: h1_hdiv_hcurl_l2.h: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
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::BernsteinBezier::baseFunctionsTri
static MoFEMErrorCode baseFunctionsTri(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:431
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
MoFEM::EntPolynomialBaseCtx::fieldName
const std::string fieldName
Definition: EntPolynomialBaseCtx.hpp:40
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
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
MoFEM::TriPolynomialBase::N_face_bubble
ublas::vector< MatrixDouble > N_face_bubble
Definition: TriPolynomialBase.hpp:39
NBFACETRI_AINSWORTH_HCURL
#define NBFACETRI_AINSWORTH_HCURL(P)
Definition: h1_hdiv_hcurl_l2.h:100
MoFEM::Hcurl_Demkowicz_FaceBaseFunctions_MBTRI
MoFEMErrorCode Hcurl_Demkowicz_FaceBaseFunctions_MBTRI(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:2433
MoFEM::EntitiesFieldData
data structure for finite element entity
Definition: EntitiesFieldData.hpp:40
MoFEM::BernsteinBezier::generateIndicesVertexTri
static MoFEMErrorCode generateIndicesVertexTri(const int N, int *alpha)
Definition: BernsteinBezier.cpp:90
MoFEM::TriPolynomialBase::getValueL2
MoFEMErrorCode getValueL2(MatrixDouble &pts)
Definition: TriPolynomialBase.cpp:321
NBFACETRI_AINSWORTH_HDIV
#define NBFACETRI_AINSWORTH_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:132
MoFEM::TriPolynomialBase::getValueHcurlAinsworthBase
MoFEMErrorCode getValueHcurlAinsworthBase(MatrixDouble &pts)
Definition: TriPolynomialBase.cpp:646
CONTINUOUS
@ CONTINUOUS
Regular field.
Definition: definitions.h:100
convert.int
int
Definition: convert.py:64
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
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::EntPolynomialBaseCtx::basePolynomialsType0
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Definition: EntPolynomialBaseCtx.hpp:27
ShapeMBTRI
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition: fem_tools.c:182
Generated by Doxygen 1.8.17 and hosted at University of Glasgow