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  auto nb_dofs_on_face =
547  3 * NBFACETRI_AINSWORTH_EDGE_HDIV(
548  AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(face_order)) +
549  NBFACETRI_AINSWORTH_FACE_HDIV(
550  AinsworthOrderHooks::broken_nbfacetri_face_hdiv(face_order));
551  if (nb_dofs_on_face) {
552  auto face_edge_dofs = NBFACETRI_AINSWORTH_EDGE_HDIV(
553  AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(face_order));
554  // three edges on face
555  for (int ee = 0; ee < 3; ee++) {
556  N_face_edge(0, ee).resize(nb_gauss_pts, 3 * face_edge_dofs, false);
557  phi_f_e[ee] = &((N_face_edge(0, ee))(0, 0));
558  }
559  auto face_bubble_dofs = NBFACETRI_AINSWORTH_FACE_HDIV(
560  AinsworthOrderHooks::broken_nbfacetri_face_hdiv(face_order));
561  N_face_bubble[0].resize(nb_gauss_pts, 3 * face_bubble_dofs, false);
562  phi_f = &*(N_face_bubble[0].data().begin());
563 
564  int face_nodes[3] = {0, 1, 2};
565  if (face_edge_dofs)
566  CHKERR Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(
567  face_nodes,
568  AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(face_order),
569  &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0), NULL, phi_f_e,
570  NULL, nb_gauss_pts, 3, base_polynomials);
571  if (face_bubble_dofs)
572  CHKERR Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(
573  face_nodes,
574  AinsworthOrderHooks::broken_nbfacetri_face_hdiv(face_order),
575  &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0), NULL, phi_f, NULL,
576  nb_gauss_pts, 3, base_polynomials);
577 
578  // set shape functions into data structure
579  if (data.dataOnEntities[MBTRI].size() != 1) {
580  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
581  }
582  data.dataOnEntities[MBTRI][0].getN(base).resize(
583  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_HDIV(face_order), false);
584 
585  auto max_face_order =
586  std::max(face_order,
587  AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(face_order));
588  auto max_face_oder =
589  std::max(max_face_order,
590  AinsworthOrderHooks::broken_nbfacetri_face_hdiv(face_order));
591 
592  int col = 0;
593  for (int oo = 0; oo < max_face_order; oo++) {
594  if (oo < AinsworthOrderHooks::broken_nbfacetri_edge_hdiv(face_order))
595  for (int ee = 0; ee < 3; ee++) {
596  for (int dd = 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
597  dd < 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++, col++) {
598  for (int gg = 0; gg < nb_gauss_pts; gg++) {
599  data.dataOnEntities[MBTRI][0].getN(base)(gg, col) =
600  N_face_edge(0, ee)(gg, dd);
601  }
602  }
603  }
604 
605  if (oo < AinsworthOrderHooks::broken_nbfacetri_face_hdiv(face_order))
606  for (int dd = 3 * NBFACETRI_AINSWORTH_FACE_HDIV(oo);
607  dd < 3 * NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++, col++) {
608  for (int gg = 0; gg < nb_gauss_pts; gg++) {
609  data.dataOnEntities[MBTRI][0].getN(base)(gg, col) =
610  N_face_bubble[0](gg, dd);
611  }
612  }
613  }
614  }
615 
616  MoFEMFunctionReturn(0);
617 }
618 
619 MoFEMErrorCode TriPolynomialBase::getValueHdivDemkowiczBase(MatrixDouble &pts) {
620  MoFEMFunctionBegin;
621 
622  EntitiesFieldData &data = cTx->dAta;
623  // set shape functions into data structure
624  if (data.dataOnEntities[MBTRI].size() != 1) {
625  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
626  }
627  const FieldApproximationBase base = cTx->bAse;
628  if (base != DEMKOWICZ_JACOBI_BASE) {
629  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
630  "This should be used only with DEMKOWICZ_JACOBI_BASE "
631  "but base is %s",
632  ApproximationBaseNames[base]);
633  }
634  int order = data.dataOnEntities[MBTRI][0].getOrder();
635  int nb_gauss_pts = pts.size2();
636  data.dataOnEntities[MBTRI][0].getN(base).resize(
637  nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
638  double *phi_f = &*data.dataOnEntities[MBTRI][0].getN(base).data().begin();
639  if (NBFACETRI_DEMKOWICZ_HDIV(order) == 0)
640  MoFEMFunctionReturnHot(0);
641  int face_nodes[3] = {0, 1, 2};
642  CHKERR Hdiv_Demkowicz_Face_MBTET_ON_FACE(
643  face_nodes, order,
644  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
645  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(), phi_f,
646  NULL, nb_gauss_pts, 3);
647 
648  MoFEMFunctionReturn(0);
649 }
650 
651 MoFEMErrorCode TriPolynomialBase::getValueHdiv(MatrixDouble &pts) {
652  MoFEMFunctionBegin;
653 
654  switch (cTx->bAse) {
655  case AINSWORTH_LEGENDRE_BASE:
656  case AINSWORTH_LOBATTO_BASE:
657  return getValueHdivAinsworthBase(pts);
658  case DEMKOWICZ_JACOBI_BASE:
659  return getValueHdivDemkowiczBase(pts);
660  default:
661  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
662  }
663 
664  MoFEMFunctionReturn(0);
665 }
666 
667 MoFEMErrorCode
668 TriPolynomialBase::getValueHcurlAinsworthBase(MatrixDouble &pts) {
669  MoFEMFunctionBegin;
670 
671  EntitiesFieldData &data = cTx->dAta;
672  const FieldApproximationBase base = cTx->bAse;
673  if (data.dataOnEntities[MBTRI].size() != 1)
674  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
675 
676  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
677  double *diffL, const int dim) =
678  cTx->basePolynomialsType0;
679 
680  int nb_gauss_pts = pts.size2();
681 
682  // Calculation H-curl on triangle faces
683  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
684  if (data.dataOnEntities[MBEDGE].size() != 3) {
685  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
686  }
687  int sense[3], order[3];
688  double *hcurl_edge_n[3];
689  double *diff_hcurl_edge_n[3];
690  for (int ee = 0; ee < 3; ee++) {
691  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
692  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
693  "data inconsistency");
694  }
695  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
696  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
697  int nb_dofs =
698  NBEDGE_AINSWORTH_HCURL(data.dataOnEntities[MBEDGE][ee].getOrder());
699  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
700  3 * nb_dofs, false);
701  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(
702  nb_gauss_pts, 2 * 3 * nb_dofs, false);
703  hcurl_edge_n[ee] =
704  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
705  diff_hcurl_edge_n[ee] =
706  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
707  }
708  CHKERR Hcurl_Ainsworth_EdgeBaseFunctions_MBTET_ON_FACE(
709  sense, order,
710  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
711  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
712  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
713  } else {
714  for (int ee = 0; ee < 3; ee++) {
715  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
716  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
717  false);
718  }
719  }
720 
721  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
722 
723  // cerr << data.dataOnEntities[MBVERTEX][0].getN(base) << endl;
724  // cerr << data.dataOnEntities[MBVERTEX][0].getDiffN(base) << endl;
725  //
726  // face
727  if (data.dataOnEntities[MBTRI].size() != 1) {
728  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
729  }
730  int order = data.dataOnEntities[MBTRI][0].getOrder();
731  int nb_dofs = NBFACETRI_AINSWORTH_HCURL(order);
732  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
733  false);
734  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
735  3 * 2 * nb_dofs, false);
736  // cerr << data.dataOnEntities[MBVERTEX][0].getDiffN(base) << endl;
737  int face_nodes[] = {0, 1, 2};
738  CHKERR Hcurl_Ainsworth_FaceFunctions_MBTET_ON_FACE(
739  face_nodes, order,
740  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
741  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
742  &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
743  &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
744  nb_gauss_pts, base_polynomials);
745  // cerr << data.dataOnEntities[MBTRI][0].getN(base) << endl;
746  } else {
747  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 0, false);
748  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
749  }
750 
751  MoFEMFunctionReturn(0);
752 }
753 
754 MoFEMErrorCode
755 TriPolynomialBase::getValueHcurlAinsworthBrokenBase(MatrixDouble &pts) {
756  MoFEMFunctionBegin;
757 
758  EntitiesFieldData &data = cTx->dAta;
759  const FieldApproximationBase base = cTx->bAse;
760  if (data.dataOnEntities[MBTRI].size() != 1)
761  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
762 
763  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
764  double *diffL, const int dim) =
765  cTx->basePolynomialsType0;
766 
767  int nb_gauss_pts = pts.size2();
768 
769  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
770 
771  // face
772  if (data.dataOnEntities[MBTRI].size() != 1) {
773  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
774  }
775 
776  int order = data.dataOnEntities[MBTRI][0].getOrder();
777  auto nb_edge_dofs = NBEDGE_AINSWORTH_HCURL(order);
778  int nb_dofs_face = NBFACETRI_AINSWORTH_HCURL(order);
779 
780  // three faces on triangle
781  int nb_dofs = 3 * nb_edge_dofs + nb_dofs_face;
782  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
783  false);
784  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
785  3 * 2 * nb_dofs, false);
786 
787  MatrixDouble edge_bases[] = {
788  MatrixDouble(nb_gauss_pts, 3 * nb_edge_dofs, false),
789  MatrixDouble(nb_gauss_pts, 3 * nb_edge_dofs, false),
790  MatrixDouble(nb_gauss_pts, 3 * nb_edge_dofs, false)};
791  MatrixDouble diff_edge_bases[] = {
792  MatrixDouble(nb_gauss_pts, 3 * 2 * nb_edge_dofs, false),
793  MatrixDouble(nb_gauss_pts, 3 * 2 * nb_edge_dofs, false),
794  MatrixDouble(nb_gauss_pts, 3 * 2 * nb_edge_dofs, false)};
795 
796  std::array<double *, 3> phi{&*edge_bases[0].data().begin(),
797  &*edge_bases[1].data().begin(),
798  &*edge_bases[2].data().begin()};
799  std::array<double *, 3> diff_phi{&*diff_edge_bases[0].data().begin(),
800  &*diff_edge_bases[1].data().begin(),
801  &*diff_edge_bases[2].data().begin()};
802 
803  std::array<int, 3> edge_order{order, order, order};
804  std::array<int, 3> edge_sense{1, 1, 1};
805  CHKERR Hcurl_Ainsworth_EdgeBaseFunctions_MBTET_ON_FACE(
806  edge_sense.data(), edge_order.data(),
807  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
808  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
809  phi.data(), diff_phi.data(), nb_gauss_pts, base_polynomials);
810 
811  MatrixDouble face_bases(nb_gauss_pts, 3 * nb_dofs_face, false);
812  MatrixDouble diff_face_bases(nb_gauss_pts, 3 * 2 * nb_dofs_face, false);
813  int face_nodes[] = {0, 1, 2};
814  CHKERR Hcurl_Ainsworth_FaceFunctions_MBTET_ON_FACE(
815  face_nodes, order,
816  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
817  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
818  &*face_bases.data().begin(), &*diff_face_bases.data().begin(),
819  nb_gauss_pts, base_polynomials);
820 
821  // edges
822  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_edge_base[] = {
823  getFTensor1FromPtr<3>(&*edge_bases[0].data().begin()),
824  getFTensor1FromPtr<3>(&*edge_bases[1].data().begin()),
825  getFTensor1FromPtr<3>(&*edge_bases[2].data().begin())};
826  FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> t_edge_diff_base[] = {
827  getFTensor2HVecFromPtr<3, 2>(&*diff_edge_bases[0].data().begin()),
828  getFTensor2HVecFromPtr<3, 2>(&*diff_edge_bases[1].data().begin()),
829  getFTensor2HVecFromPtr<3, 2>(&*diff_edge_bases[2].data().begin())};
830 
831  // face
832  auto t_vol_base = getFTensor1FromPtr<3>(&*face_bases.data().begin());
833  auto t_vol_diff_base =
834  getFTensor2HVecFromPtr<3, 2>(&*diff_face_bases.data().begin());
835 
836  auto t_base = getFTensor1FromPtr<3>(
837  &*data.dataOnEntities[MBTRI][0].getN(base).data().begin());
838  auto t_diff_base = getFTensor2HVecFromPtr<3, 2>(
839  &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin());
840 
841  FTENSOR_INDEX(3, i);
842  FTENSOR_INDEX(2, j);
843 
844  for (int gg = 0; gg != nb_gauss_pts; gg++) {
845  for (int oo = 0; oo < order; oo++) {
846  // edges
847  for (int dd = NBEDGE_AINSWORTH_HCURL(oo);
848  dd != NBEDGE_AINSWORTH_HCURL(oo + 1); ++dd) {
849  for (int ee = 0; ee != 3; ++ee) {
850  t_base(i) = t_edge_base[ee](i);
851  ++t_base;
852  ++t_edge_base[ee];
853  t_diff_base(i, j) = t_edge_diff_base[ee](i, j);
854  ++t_diff_base;
855  ++t_edge_diff_base[ee];
856  }
857  }
858  // face
859  for (int dd = NBFACETRI_AINSWORTH_HCURL(oo);
860  dd != NBFACETRI_AINSWORTH_HCURL(oo + 1); ++dd) {
861  t_base(i) = t_vol_base(i);
862  ++t_base;
863  ++t_vol_base;
864  t_diff_base(i, j) = t_vol_diff_base(i, j);
865  ++t_diff_base;
866  ++t_vol_diff_base;
867  }
868  }
869  }
870 
871  } else {
872  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 0, false);
873  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
874  }
875 
876  MoFEMFunctionReturn(0);
877 }
878 
879 MoFEMErrorCode
880 TriPolynomialBase::getValueHcurlDemkowiczBase(MatrixDouble &pts) {
881  MoFEMFunctionBegin;
882 
883  EntitiesFieldData &data = cTx->dAta;
884  const FieldApproximationBase base = cTx->bAse;
885 
886  int nb_gauss_pts = pts.size2();
887 
888  // Calculation H-curl on triangle faces
889  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
890 
891  if (data.dataOnEntities[MBEDGE].size() != 3)
892  SETERRQ1(
893  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
894  "wrong number of data structures on edges, should be three but is %d",
895  data.dataOnEntities[MBEDGE].size());
896 
897  int sense[3], order[3];
898  double *hcurl_edge_n[3];
899  double *diff_hcurl_edge_n[3];
900 
901  for (int ee = 0; ee != 3; ++ee) {
902 
903  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0)
904  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
905  "orientation (sense) of edge is not set");
906 
907  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
908  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
909  int nb_dofs =
910  NBEDGE_DEMKOWICZ_HCURL(data.dataOnEntities[MBEDGE][ee].getOrder());
911  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
912  3 * nb_dofs, false);
913  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(
914  nb_gauss_pts, 2 * 3 * nb_dofs, false);
915 
916  hcurl_edge_n[ee] =
917  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
918  diff_hcurl_edge_n[ee] =
919  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
920  }
921 
922  CHKERR Hcurl_Demkowicz_EdgeBaseFunctions_MBTRI(
923  sense, order,
924  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
925  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
926  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts);
927 
928  } else {
929 
930  // No DOFs on faces, resize base function matrices, indicating that no
931  // dofs on them.
932  for (int ee = 0; ee != 3; ++ee) {
933  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
934  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
935  false);
936  }
937  }
938 
939  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
940 
941  // face
942  if (data.dataOnEntities[MBTRI].size() != 1)
943  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
944  "No data struture to keep base functions on face");
945 
946  int order = data.dataOnEntities[MBTRI][0].getOrder();
947  int nb_dofs = NBFACETRI_DEMKOWICZ_HCURL(order);
948  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
949  false);
950  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
951  3 * 2 * nb_dofs, false);
952 
953  int face_nodes[] = {0, 1, 2};
954  CHKERR Hcurl_Demkowicz_FaceBaseFunctions_MBTRI(
955  face_nodes, order,
956  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
957  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
958  &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
959  &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
960  nb_gauss_pts);
961 
962  } else {
963 
964  // No DOFs on faces, resize base function matrices, indicating that no
965  // dofs on them.
966  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 0, false);
967  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
968  }
969 
970  MoFEMFunctionReturn(0);
971 }
972 
973 MoFEMErrorCode
974 TriPolynomialBase::getValueHcurlDemkowiczBrokenBase(MatrixDouble &pts) {
975  MoFEMFunctionBegin;
976 
977  EntitiesFieldData &data = cTx->dAta;
978  const FieldApproximationBase base = cTx->bAse;
979 
980  int nb_gauss_pts = pts.size2();
981 
982  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
983 
984  // face
985  if (data.dataOnEntities[MBTRI].size() != 1)
986  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
987  "No data struture to keep base functions on face");
988 
989  int order = data.dataOnEntities[MBTRI][0].getOrder();
990  int nb_edge_dofs = NBEDGE_DEMKOWICZ_HCURL(order);
991  int nb_volume_dofs = NBFACETRI_DEMKOWICZ_HCURL(order);
992  int nb_dofs = 3 * nb_edge_dofs + nb_volume_dofs;
993  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
994  false);
995  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
996  3 * 2 * nb_dofs, false);
997 
998  MatrixDouble edge_bases[] = {
999  MatrixDouble(nb_gauss_pts, 3 * nb_edge_dofs, false),
1000  MatrixDouble(nb_gauss_pts, 3 * nb_edge_dofs, false),
1001  MatrixDouble(nb_gauss_pts, 3 * nb_edge_dofs, false)};
1002  MatrixDouble diff_edge_bases[] = {
1003  MatrixDouble(nb_gauss_pts, 3 * 2 * nb_edge_dofs, false),
1004  MatrixDouble(nb_gauss_pts, 3 * 2 * nb_edge_dofs, false),
1005  MatrixDouble(nb_gauss_pts, 3 * 2 * nb_edge_dofs, false)};
1006 
1007  std::array<int, 3> edge_order{order, order, order};
1008  std::array<int, 3> edge_sense{1, 1, 1};
1009  std::array<double *, 3> phi{&*edge_bases[0].data().begin(),
1010  &*edge_bases[1].data().begin(),
1011  &*edge_bases[2].data().begin()};
1012  std::array<double *, 3> diff_phi{&*diff_edge_bases[0].data().begin(),
1013  &*diff_edge_bases[1].data().begin(),
1014  &*diff_edge_bases[2].data().begin()};
1015  CHKERR Hcurl_Demkowicz_EdgeBaseFunctions_MBTRI(
1016  edge_sense.data(), edge_order.data(),
1017  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1018  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1019  phi.data(), diff_phi.data(), nb_gauss_pts);
1020 
1021  MatrixDouble face_bases(nb_gauss_pts, 3 * nb_volume_dofs, false);
1022  MatrixDouble diff_face_bases(nb_gauss_pts, 3 * 2 * nb_volume_dofs, false);
1023  int face_nodes[] = {0, 1, 2};
1024  CHKERR Hcurl_Demkowicz_FaceBaseFunctions_MBTRI(
1025  face_nodes, order,
1026  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1027  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1028  &*face_bases.data().begin(), &*diff_face_bases.data().begin(),
1029  nb_gauss_pts);
1030 
1031  // edges
1032  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_edge_base[] = {
1033  getFTensor1FromPtr<3>(&*edge_bases[0].data().begin()),
1034  getFTensor1FromPtr<3>(&*edge_bases[1].data().begin()),
1035  getFTensor1FromPtr<3>(&*edge_bases[2].data().begin())};
1036  FTensor::Tensor2<FTensor::PackPtr<double *, 6>, 3, 2> t_edge_diff_base[] = {
1037  getFTensor2HVecFromPtr<3, 2>(&*diff_edge_bases[0].data().begin()),
1038  getFTensor2HVecFromPtr<3, 2>(&*diff_edge_bases[1].data().begin()),
1039  getFTensor2HVecFromPtr<3, 2>(&*diff_edge_bases[2].data().begin())};
1040  // face
1041  auto t_vol_base = getFTensor1FromPtr<3>(&*face_bases.data().begin());
1042  auto t_vol_diff_base =
1043  getFTensor2HVecFromPtr<3, 2>(&*diff_face_bases.data().begin());
1044 
1045  auto t_base = getFTensor1FromPtr<3>(
1046  &*data.dataOnEntities[MBTRI][0].getN(base).data().begin());
1047  auto t_diff_base = getFTensor2HVecFromPtr<3, 2>(
1048  &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin());
1049 
1050  FTENSOR_INDEX(3, i);
1051  FTENSOR_INDEX(2, j);
1052 
1053 
1054  for (int gg = 0; gg != nb_gauss_pts; gg++) {
1055  for (int oo = 0; oo < order; oo++) {
1056  // edges
1057  for (int dd = NBEDGE_DEMKOWICZ_HCURL(oo);
1058  dd != NBEDGE_DEMKOWICZ_HCURL(oo + 1); ++dd) {
1059  for (int ee = 0; ee != 3; ++ee) {
1060  t_base(i) = t_edge_base[ee](i);
1061  ++t_base;
1062  ++t_edge_base[ee];
1063  t_diff_base(i, j) = t_edge_diff_base[ee](i, j);
1064  ++t_diff_base;
1065  ++t_edge_diff_base[ee];
1066  }
1067  }
1068  // faces
1069  for (int dd = NBFACETRI_DEMKOWICZ_HCURL(oo);
1070  dd != NBFACETRI_DEMKOWICZ_HCURL(oo + 1); ++dd) {
1071  t_base(i) = t_vol_base(i);
1072  ++t_base;
1073  ++t_vol_base;
1074  t_diff_base(i, j) = t_vol_diff_base(i, j);
1075  ++t_diff_base;
1076  ++t_vol_diff_base;
1077  }
1078  }
1079  }
1080 
1081  } else {
1082 
1083  // No DOFs on faces, resize base function matrices, indicating that no
1084  // dofs on them.
1085  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 0, false);
1086  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
1087  }
1088 
1089  MoFEMFunctionReturn(0);
1090 }
1091 
1092 MoFEMErrorCode TriPolynomialBase::getValueHcurl(MatrixDouble &pts) {
1093  MoFEMFunctionBegin;
1094  switch (cTx->spaceContinuity) {
1095  case CONTINUOUS:
1096 
1097  switch (cTx->bAse) {
1098  case AINSWORTH_LEGENDRE_BASE:
1099  case AINSWORTH_LOBATTO_BASE:
1100  CHKERR getValueHcurlAinsworthBase(pts);
1101  break;
1102  case DEMKOWICZ_JACOBI_BASE:
1103  CHKERR getValueHcurlDemkowiczBase(pts);
1104  break;
1105  default:
1106  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1107  }
1108 
1109  break;
1110 
1111  case DISCONTINUOUS:
1112  switch (cTx->bAse) {
1113  case AINSWORTH_LEGENDRE_BASE:
1114  case AINSWORTH_LOBATTO_BASE:
1115  return getValueHcurlAinsworthBrokenBase(pts);
1116  case DEMKOWICZ_JACOBI_BASE:
1117  return getValueHcurlDemkowiczBrokenBase(pts);
1118  default:
1119  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1120  }
1121  break;
1122 
1123  default:
1124  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown continuity");
1125  }
1126 
1127  MoFEMFunctionReturn(0);
1128 }
1129 
1130 MoFEMErrorCode
1131 TriPolynomialBase::getValue(MatrixDouble &pts,
1132  boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
1133  MoFEMFunctionBegin;
1134 
1135  cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
1136 
1137  int nb_gauss_pts = pts.size2();
1138  if (!nb_gauss_pts) {
1139  MoFEMFunctionReturnHot(0);
1140  }
1141 
1142  if (pts.size1() < 2)
1143  SETERRQ(
1144  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1145  "Wrong dimension of pts, should be at least 3 rows with coordinates");
1146 
1147  const FieldApproximationBase base = cTx->bAse;
1148  EntitiesFieldData &data = cTx->dAta;
1149 
1150  if (base != AINSWORTH_BERNSTEIN_BEZIER_BASE) {
1151  if (cTx->copyNodeBase == LASTBASE) {
1152  data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 3,
1153  false);
1154  CHKERR ShapeMBTRI(
1155  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1156  &pts(0, 0), &pts(1, 0), nb_gauss_pts);
1157  data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(3, 2, false);
1158  CHKERR ShapeDiffMBTRI(
1159  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
1160  } else {
1161  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
1162  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
1163  data.dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(base) =
1164  data.dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(cTx->copyNodeBase);
1165  }
1166  if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
1167  (unsigned int)nb_gauss_pts) {
1168  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1169  "Base functions or nodes has wrong number of integration points "
1170  "for base %s",
1171  ApproximationBaseNames[base]);
1172  }
1173  }
1174  auto set_node_derivative_for_all_gauss_pts = [&]() {
1175  MoFEMFunctionBegin;
1176  // In linear geometry derivatives are constant,
1177  // this in expense of efficiency makes implementation
1178  // consistent between vertices and other types of entities
1179  data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(nb_gauss_pts, 6,
1180  false);
1181  for (int gg = 0; gg != nb_gauss_pts; ++gg)
1182  std::copy(Tools::diffShapeFunMBTRI.begin(),
1183  Tools::diffShapeFunMBTRI.end(),
1184  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(gg, 0));
1185  MoFEMFunctionReturn(0);
1186  };
1187 
1188  CHKERR set_node_derivative_for_all_gauss_pts();
1189 
1190  switch (cTx->spaceContinuity) {
1191  case CONTINUOUS:
1192 
1193  switch (cTx->sPace) {
1194  case H1:
1195  CHKERR getValueH1(pts);
1196  break;
1197  case HDIV:
1198  CHKERR getValueHdiv(pts);
1199  break;
1200  case HCURL:
1201  CHKERR getValueHcurl(pts);
1202  break;
1203  case L2:
1204  CHKERR getValueL2(pts);
1205  break;
1206  default:
1207  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
1208  }
1209 
1210  break;
1211 
1212  case DISCONTINUOUS:
1213  switch (cTx->sPace) {
1214  case HCURL:
1215  CHKERR getValueHcurl(pts);
1216  break;
1217  default:
1218  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
1219  }
1220  break;
1221 
1222  default:
1223  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown continuity");
1224  }
1225 
1226  MoFEMFunctionReturn(0);
1227 }
1228 
1229 MoFEMErrorCode TriPolynomialBase::setDofsSideMap(
1230 
1231  const FieldSpace space, const FieldContinuity continuity,
1232  const FieldApproximationBase base, DofsSideMap &dofs_side_map
1233 
1234 ) {
1235  MoFEMFunctionBegin;
1236 
1237  switch (continuity) {
1238  case DISCONTINUOUS:
1239 
1240  switch (space) {
1241  case HCURL:
1242  CHKERR setDofsSideMapHcurl(space, continuity, base, dofs_side_map);
1243  break;
1244  default:
1245  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space %s",
1246  FieldSpaceNames[space]);
1247  }
1248  break;
1249 
1250  default:
1251  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1252  "Unknown (or not implemented) continuity");
1253  }
1254 
1255  MoFEMFunctionReturn(0);
1256 }
1257 
1258 MoFEMErrorCode TriPolynomialBase::setDofsSideMapHcurl(
1259  const FieldSpace space, const FieldContinuity continuity,
1260  const FieldApproximationBase base, DofsSideMap &dofs_side_map) {
1261  MoFEMFunctionBegin;
1262 
1263  // That has to be consistent with implementation of getValueHdiv for
1264  // particular base functions.
1265 
1266  auto set_ainsworth = [&dofs_side_map]() {
1267  MoFEMFunctionBegin;
1268 
1269  dofs_side_map.clear();
1270 
1271  int dof = 0;
1272  for (int oo = 0; oo < Field::maxBrokenDofsOrder; oo++) {
1273 
1274  // edges
1275  for (int dd = NBEDGE_AINSWORTH_HCURL(oo);
1276  dd != NBEDGE_AINSWORTH_HCURL(oo + 1); ++dd) {
1277  for (int ee = 0; ee != 3; ++ee) {
1278  dofs_side_map.insert(DofsSideMapData{MBEDGE, ee, dof});
1279  ++dof;
1280  }
1281  }
1282  // face
1283  for (int dd = NBFACETRI_AINSWORTH_HCURL(oo);
1284  dd != NBFACETRI_AINSWORTH_HCURL(oo + 1); ++dd) {
1285  dofs_side_map.insert(DofsSideMapData{MBTRI, 0, dof});
1286  ++dof;
1287  }
1288  }
1289 
1290  MoFEMFunctionReturn(0);
1291  };
1292 
1293  auto set_demkowicz = [&dofs_side_map]() {
1294  MoFEMFunctionBegin;
1295 
1296  dofs_side_map.clear();
1297 
1298  int dof = 0;
1299  for (int oo = 0; oo < Field::maxBrokenDofsOrder; oo++) {
1300 
1301  // edges
1302  for (int dd = NBEDGE_DEMKOWICZ_HCURL(oo);
1303  dd != NBEDGE_DEMKOWICZ_HCURL(oo + 1); ++dd) {
1304  for (int ee = 0; ee != 3; ++ee) {
1305  dofs_side_map.insert(DofsSideMapData{MBEDGE, ee, dof});
1306  ++dof;
1307  }
1308  }
1309  // faces
1310  for (int dd = NBFACETRI_DEMKOWICZ_HCURL(oo);
1311  dd != NBFACETRI_DEMKOWICZ_HCURL(oo + 1); ++dd) {
1312  dofs_side_map.insert(DofsSideMapData{MBTRI, 0, dof});
1313  ++dof;
1314  }
1315  }
1316 
1317  MoFEMFunctionReturn(0);
1318  };
1319 
1320  switch (continuity) {
1321  case DISCONTINUOUS:
1322  switch (base) {
1323  case AINSWORTH_LEGENDRE_BASE:
1324  case AINSWORTH_LOBATTO_BASE:
1325  return set_ainsworth();
1326  case DEMKOWICZ_JACOBI_BASE:
1327  return set_demkowicz();
1328  default:
1329  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1330  }
1331  break;
1332 
1333  default:
1334  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
1335  "Unknown (or not implemented) continuity");
1336  }
1337 
1338  MoFEMFunctionReturn(0);
1339 }
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:47
H1
@ H1
continuous field
Definition: definitions.h:85
NBEDGE_H1
#define NBEDGE_H1(P)
Number of base function on edge for H1 space.
Definition: h1_hdiv_hcurl_l2.h:55
MoFEM::Hdiv_Demkowicz_Face_MBTET_ON_FACE
MoFEMErrorCode Hdiv_Demkowicz_Face_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb)
Definition: Hdiv.cpp:634
LASTBASE
@ LASTBASE
Definition: definitions.h:69
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
MoFEM::TriPolynomialBase::getValueHdivDemkowiczBase
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)
Definition: TriPolynomialBase.cpp:619
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
MoFEM::BaseFunction::DofsSideMap
multi_index_container< DofsSideMapData, indexed_by< ordered_non_unique< tag< TypeSide_mi_tag >, composite_key< DofsSideMapData, member< DofsSideMapData, EntityType, &DofsSideMapData::type >, member< DofsSideMapData, int, &DofsSideMapData::side > >>, ordered_unique< tag< EntDofIdx_mi_tag >, member< DofsSideMapData, int, &DofsSideMapData::dof > > > > DofsSideMap
Map entity stype and side to element/entity dof index.
Definition: BaseFunction.hpp:73
NOBASE
@ NOBASE
Definition: definitions.h:59
MoFEM::TriPolynomialBase::N_face_edge
ublas::matrix< MatrixDouble > N_face_edge
Definition: TriPolynomialBase.hpp:57
MoFEM::Field::maxBrokenDofsOrder
static constexpr int maxBrokenDofsOrder
max number of broken dofs
Definition: FieldMultiIndices.hpp:282
ApproximationBaseNames
const static char *const ApproximationBaseNames[]
Definition: definitions.h:72
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
NBFACETRI_AINSWORTH_EDGE_HDIV
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:130
MoFEM::Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE
MoFEMErrorCode Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f_e[3], double *diff_phi_f_e[3], int gdim, int nb, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base functions, Edge-based face functions by Ainsworth .
Definition: Hdiv.cpp:47
phi
static double phi
Definition: poisson_2d_dis_galerkin.cpp:30
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:651
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
FTENSOR_INDEX
#define FTENSOR_INDEX(DIM, I)
Definition: Templates.hpp:2013
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
MoFEM::AinsworthOrderHooks::broken_nbfacetri_edge_hdiv
static boost::function< int(int)> broken_nbfacetri_edge_hdiv
Definition: Hdiv.hpp:25
ShapeDiffMBTRI
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:194
FieldSpace
FieldSpace
approximation spaces
Definition: definitions.h:82
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:880
MoFEM::TriPolynomialBase
Calculate base functions on triangle.
Definition: TriPolynomialBase.hpp:16
FieldContinuity
FieldContinuity
Field continuity.
Definition: definitions.h:99
NBFACETRI_DEMKOWICZ_HDIV
#define NBFACETRI_DEMKOWICZ_HDIV(P)
Definition: h1_hdiv_hcurl_l2.h:139
MoFEM::Types::MatrixInt
UBlasMatrix< int > MatrixInt
Definition: Types.hpp:76
MoFEM::TriPolynomialBase::getValueL2BernsteinBezierBase
MoFEMErrorCode getValueL2BernsteinBezierBase(MatrixDouble &pts)
Definition: TriPolynomialBase.cpp:373
MoFEM::TriPolynomialBase::setDofsSideMap
static MoFEMErrorCode setDofsSideMap(const FieldSpace space, const FieldContinuity continuity, const FieldApproximationBase base, DofsSideMap &)
Set map of dof to side number.
Definition: TriPolynomialBase.cpp:1229
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:1131
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:755
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:174
MoFEM::AinsworthOrderHooks::broken_nbfacetri_face_hdiv
static boost::function< int(int)> broken_nbfacetri_face_hdiv
Definition: Hdiv.hpp:26
MoFEM::TriPolynomialBase::getValueHcurl
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
Definition: TriPolynomialBase.cpp:1092
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
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:974
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
FieldSpaceNames
const static char *const FieldSpaceNames[]
Definition: definitions.h:92
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::TriPolynomialBase::getValueH1BernsteinBezierBase
MoFEMErrorCode getValueH1BernsteinBezierBase(MatrixDouble &pts)
Definition: TriPolynomialBase.cpp:35
MoFEM::TriPolynomialBase::setDofsSideMapHcurl
static MoFEMErrorCode setDofsSideMapHcurl(const FieldSpace space, const FieldContinuity continuity, const FieldApproximationBase base, DofsSideMap &dofs_side_map)
Definition: TriPolynomialBase.cpp:1258
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:957
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::BaseFunction::DofsSideMapData
Definition: BaseFunction.hpp:42
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:58
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:668
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