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