v0.9.1
TriPolynomialBase.cpp
Go to the documentation of this file.
1 /** \file TriPolynomialBase.cpp
2 \brief Implementation of bases on triangle
3 
4 */
5 
6 /* This file is part of MoFEM.
7  * MoFEM is free software: you can redistribute it and/or modify it under
8  * the terms of the GNU Lesser General Public License as published by the
9  * Free Software Foundation, either version 3 of the License, or (at your
10  * option) any later version.
11  *
12  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
13  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  * License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
19 
20 using namespace MoFEM;
21 
24 
27  BaseFunctionUnknownInterface **iface) const {
29  *iface = NULL;
30  if (uuid == IDD_TET_BASE_FUNCTION) {
31  *iface = const_cast<TriPolynomialBase *>(this);
33  } else {
34  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "wrong interference");
35  }
38 }
39 
42 
43  switch (cTx->bAse) {
47  break;
50  break;
51  default:
52  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
53  }
54 
56 }
57 
62  const std::string &field_name = cTx->fieldName;
63  int nb_gauss_pts = pts.size2();
64 
65  auto get_alpha = [field_name](auto &data) -> MatrixInt & {
66  auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
67  if (!ptr)
68  ptr.reset(new MatrixInt());
69  return *ptr;
70  };
71 
72  auto get_base = [field_name](auto &data) -> MatrixDouble & {
73  auto &ptr = data.getBBNSharedPtr(field_name);
74  if (!ptr)
75  ptr.reset(new MatrixDouble());
76  return *ptr;
77  };
78 
79  auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
80  auto &ptr = data.getBBDiffNSharedPtr(field_name);
81  if (!ptr)
82  ptr.reset(new MatrixDouble());
83  return *ptr;
84  };
85 
86  auto get_alpha_by_name_ptr =
87  [](auto &data,
88  const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
89  return data.getBBAlphaIndicesSharedPtr(field_name);
90  };
91 
92  auto get_base_by_name_ptr =
93  [](auto &data,
94  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
95  return data.getBBNSharedPtr(field_name);
96  };
97 
98  auto get_diff_base_by_name_ptr =
99  [](auto &data,
100  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
101  return data.getBBDiffNSharedPtr(field_name);
102  };
103 
104  auto get_alpha_by_order_ptr =
105  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
106  return data.getBBAlphaIndicesByOrderSharedPtr(o);
107  };
108 
109  auto get_base_by_order_ptr =
110  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
111  return data.getBBNByOrderSharedPtr(o);
112  };
113 
114  auto get_diff_base_by_order_ptr =
115  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
116  return data.getBBDiffNByOrderSharedPtr(o);
117  };
118 
119  auto &vert_get_n = get_base(data.dataOnEntities[MBVERTEX][0]);
120  auto &vert_get_diff_n = get_diff_base(data.dataOnEntities[MBVERTEX][0]);
121  vert_get_n.resize(nb_gauss_pts, 3, false);
122  vert_get_diff_n.resize(nb_gauss_pts, 6, false);
123 
124  if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
125  (unsigned int)nb_gauss_pts)
126  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
127  "Base functions or nodes has wrong number of integration points "
128  "for base %s",
130  auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
131 
132  auto &vertex_alpha = get_alpha(data.dataOnEntities[MBVERTEX][0]);
133  vertex_alpha.resize(3, 3, false);
134  vertex_alpha.clear();
135  for (int n = 0; n != 3; ++n)
136  vertex_alpha(n, n) = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n];
137 
139  1, lambda.size1(), vertex_alpha.size1(), &vertex_alpha(0, 0),
140  &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &vert_get_n(0, 0),
141  &vert_get_diff_n(0, 0));
142  for (int n = 0; n != 3; ++n) {
143  const int f = boost::math::factorial<double>(
144  data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n]);
145  for (int g = 0; g != nb_gauss_pts; ++g) {
146  vert_get_n(g, n) *= f;
147  for (int d = 0; d != 2; ++d)
148  vert_get_diff_n(g, 2 * n + d) *= f;
149  }
150  }
151 
152  // edges
153  if (data.spacesOnEntities[MBEDGE].test(H1)) {
154  if (data.dataOnEntities[MBEDGE].size() != 3)
155  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
156  "Wrong size of ent data");
157 
158  constexpr int edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
159  for (int ee = 0; ee != 3; ++ee) {
160  auto &ent_data = data.dataOnEntities[MBEDGE][ee];
161 
162  if (ent_data.getSense() == 0)
163  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
164  "Sense of the edge unknown");
165 
166  const int sense = ent_data.getSense();
167  const int order = ent_data.getDataOrder();
168  const int nb_dofs = NBEDGE_H1(ent_data.getDataOrder());
169 
170  if (nb_dofs) {
171  if (get_alpha_by_order_ptr(ent_data, order)) {
172  get_alpha_by_name_ptr(ent_data, field_name) =
173  get_alpha_by_order_ptr(ent_data, order);
174  get_base_by_name_ptr(ent_data, field_name) =
175  get_base_by_order_ptr(ent_data, order);
176  get_diff_base_by_name_ptr(ent_data, field_name) =
177  get_diff_base_by_order_ptr(ent_data, order);
178  } else {
179  auto &get_n = get_base(ent_data);
180  auto &get_diff_n = get_diff_base(ent_data);
181  get_n.resize(nb_gauss_pts, nb_dofs, false);
182  get_diff_n.resize(nb_gauss_pts, 2 * nb_dofs, false);
183 
184  auto &edge_alpha = get_alpha(data.dataOnEntities[MBEDGE][ee]);
185  edge_alpha.resize(nb_dofs, 3, false);
187  &edge_alpha(0, 0));
188  if (sense == -1) {
189  for (int i = 0; i != edge_alpha.size1(); ++i) {
190  int a = edge_alpha(i, edges_nodes[ee][0]);
191  edge_alpha(i, edges_nodes[ee][0]) =
192  edge_alpha(i, edges_nodes[ee][1]);
193  edge_alpha(i, edges_nodes[ee][1]) = a;
194  }
195  }
197  order, lambda.size1(), edge_alpha.size1(), &edge_alpha(0, 0),
198  &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &get_n(0, 0),
199  &get_diff_n(0, 0));
200 
201  get_alpha_by_order_ptr(ent_data, order) =
202  get_alpha_by_name_ptr(ent_data, field_name);
203  get_base_by_order_ptr(ent_data, order) =
204  get_base_by_name_ptr(ent_data, field_name);
205  get_diff_base_by_order_ptr(ent_data, order) =
206  get_diff_base_by_name_ptr(ent_data, field_name);
207  }
208  }
209  }
210  } else {
211  for (int ee = 0; ee != 3; ++ee) {
212  auto &ent_data = data.dataOnEntities[MBEDGE][ee];
213  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
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, 0, false);
217  get_diff_n.resize(nb_gauss_pts, 0, false);
218  }
219  }
220 
221  if (data.spacesOnEntities[MBTRI].test(H1)) {
222  if (data.dataOnEntities[MBTRI].size() != 1)
223  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
224  "Wrong size ent of ent data");
225 
226  auto &ent_data = data.dataOnEntities[MBTRI][0];
227  const int order = ent_data.getDataOrder();
228  const int nb_dofs = NBFACETRI_H1(order);
229  if (get_alpha_by_order_ptr(ent_data, order)) {
230  get_alpha_by_name_ptr(ent_data, field_name) =
231  get_alpha_by_order_ptr(ent_data, order);
232  get_base_by_name_ptr(ent_data, field_name) =
233  get_base_by_order_ptr(ent_data, order);
234  get_diff_base_by_name_ptr(ent_data, field_name) =
235  get_diff_base_by_order_ptr(ent_data, order);
236  } else {
237 
238  auto &get_n = get_base(ent_data);
239  auto &get_diff_n = get_diff_base(ent_data);
240  get_n.resize(nb_gauss_pts, nb_dofs, false);
241  get_diff_n.resize(nb_gauss_pts, 2 * nb_dofs, false);
242  if (nb_dofs) {
243  auto &tri_alpha = get_alpha(ent_data);
244  tri_alpha.resize(nb_dofs, 3, false);
245 
248  order, lambda.size1(), tri_alpha.size1(), &tri_alpha(0, 0),
249  &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &get_n(0, 0),
250  &get_diff_n(0, 0));
251 
252  get_alpha_by_order_ptr(ent_data, order) =
253  get_alpha_by_name_ptr(ent_data, field_name);
254  get_base_by_order_ptr(ent_data, order) =
255  get_base_by_name_ptr(ent_data, field_name);
256  get_diff_base_by_order_ptr(ent_data, order) =
257  get_diff_base_by_name_ptr(ent_data, field_name);
258  }
259  }
260  } else {
261  auto &ent_data = data.dataOnEntities[MBTRI][0];
262  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
263  auto &get_n = get_base(ent_data);
264  auto &get_diff_n = get_diff_base(ent_data);
265  get_n.resize(nb_gauss_pts, 0, false);
266  get_diff_n.resize(nb_gauss_pts, 0, false);
267  }
268 
270 }
271 
274 
276  const FieldApproximationBase base = cTx->bAse;
277  if (cTx->basePolynomialsType0 == NULL)
278  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
279  "Polynomial type not set");
280  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
281  double *diffL, const int dim) =
283 
284  int nb_gauss_pts = pts.size2();
285 
286  if (data.spacesOnEntities[MBEDGE].test(H1)) {
287  // edges
288  if (data.dataOnEntities[MBEDGE].size() != 3)
289  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
290 
291  int sense[3], order[3];
292  double *H1edgeN[3], *diffH1edgeN[3];
293  for (int ee = 0; ee < 3; ee++) {
294  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0)
295  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
296  "sense of the edge unknown");
297 
298  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
299  order[ee] = data.dataOnEntities[MBEDGE][ee].getDataOrder();
300  int nb_dofs = NBEDGE_H1(data.dataOnEntities[MBEDGE][ee].getDataOrder());
301  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
302  false);
303  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
304  2 * nb_dofs, false);
305  H1edgeN[ee] = &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
306  diffH1edgeN[ee] =
307  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
308  }
310  sense, order,
311  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
312  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
313  H1edgeN, diffH1edgeN, nb_gauss_pts, base_polynomials);
314  }
315 
316  if (data.spacesOnEntities[MBTRI].test(H1)) {
317  // face
318  if (data.dataOnEntities[MBTRI].size() != 1) {
319  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
320  }
321  int nb_dofs = NBFACETRI_H1(data.dataOnEntities[MBTRI][0].getDataOrder());
322  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, nb_dofs,
323  false);
324  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
325  2 * nb_dofs, false);
326  const int face_nodes[] = {0, 1, 2};
328  face_nodes, data.dataOnEntities[MBTRI][0].getDataOrder(),
329  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
330  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
331  &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
332  &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
333  nb_gauss_pts, base_polynomials);
334  }
335 
337 }
338 
341 
342  switch (cTx->bAse) {
346  break;
349  break;
350  default:
351  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
352  }
353 
355 }
356 
359 
361  const FieldApproximationBase base = cTx->bAse;
362  if (cTx->basePolynomialsType0 == NULL)
363  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
364  "Polynomial type not set");
365  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
366  double *diffL, const int dim) =
368 
369  int nb_gauss_pts = pts.size2();
370 
371  data.dataOnEntities[MBTRI][0].getN(base).resize(
372  nb_gauss_pts, NBFACETRI_L2(data.dataOnEntities[MBTRI][0].getDataOrder()),
373  false);
374  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(
375  nb_gauss_pts,
376  2 * NBFACETRI_L2(data.dataOnEntities[MBTRI][0].getDataOrder()), false);
377 
379  data.dataOnEntities[MBTRI][0].getDataOrder(),
380  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
381  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
382  &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
383  &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
384  nb_gauss_pts, base_polynomials);
385 
387 }
388 
392 
394  const std::string &field_name = cTx->fieldName;
395  int nb_gauss_pts = pts.size2();
396 
397  auto get_alpha = [field_name](auto &data) -> MatrixInt & {
398  auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
399  if (!ptr)
400  ptr.reset(new MatrixInt());
401  return *ptr;
402  };
403 
404  auto get_base = [field_name](auto &data) -> MatrixDouble & {
405  auto &ptr = data.getBBNSharedPtr(field_name);
406  if (!ptr)
407  ptr.reset(new MatrixDouble());
408  return *ptr;
409  };
410 
411  auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
412  auto &ptr = data.getBBDiffNSharedPtr(field_name);
413  if (!ptr)
414  ptr.reset(new MatrixDouble());
415  return *ptr;
416  };
417 
418  auto get_alpha_by_name_ptr =
419  [](auto &data,
420  const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
421  return data.getBBAlphaIndicesSharedPtr(field_name);
422  };
423 
424  auto get_base_by_name_ptr =
425  [](auto &data,
426  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
427  return data.getBBNSharedPtr(field_name);
428  };
429 
430  auto get_diff_base_by_name_ptr =
431  [](auto &data,
432  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
433  return data.getBBDiffNSharedPtr(field_name);
434  };
435 
436  auto get_alpha_by_order_ptr =
437  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
438  return data.getBBAlphaIndicesByOrderSharedPtr(o);
439  };
440 
441  auto get_base_by_order_ptr =
442  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
443  return data.getBBNByOrderSharedPtr(o);
444  };
445 
446  auto get_diff_base_by_order_ptr =
447  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
448  return data.getBBDiffNByOrderSharedPtr(o);
449  };
450 
451  if (data.spacesOnEntities[MBTRI].test(L2)) {
452  if (data.dataOnEntities[MBTRI].size() != 1)
453  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
454  "Wrong size ent of ent data");
455 
456  if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
457  (unsigned int)nb_gauss_pts)
458  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
459  "Base functions or nodes has wrong number of integration points "
460  "for base %s",
462  auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
463 
464  auto &ent_data = data.dataOnEntities[MBTRI][0];
465  const int order = ent_data.getDataOrder();
466  const int nb_dofs = NBFACETRI_L2(order);
467 
468  if (get_alpha_by_order_ptr(ent_data, order)) {
469  get_alpha_by_name_ptr(ent_data, field_name) =
470  get_alpha_by_order_ptr(ent_data, order);
471  get_base_by_name_ptr(ent_data, field_name) =
472  get_base_by_order_ptr(ent_data, order);
473  get_diff_base_by_name_ptr(ent_data, field_name) =
474  get_diff_base_by_order_ptr(ent_data, order);
475  } else {
476 
477  auto &get_n = get_base(ent_data);
478  auto &get_diff_n = get_diff_base(ent_data);
479  get_n.resize(nb_gauss_pts, nb_dofs, false);
480  get_diff_n.resize(nb_gauss_pts, 2 * nb_dofs, false);
481  if (nb_dofs) {
482 
483  if (order == 0) {
484 
485  if (nb_dofs != 1)
486  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
487  "Inconsistent number of DOFs");
488 
489  auto &tri_alpha = get_alpha(ent_data);
490  tri_alpha.clear();
491  get_n(0, 0) = 1;
492  get_diff_n.clear();
493 
494  } else {
495 
496  if (nb_dofs != 3 + 3 * NBEDGE_H1(order) + NBFACETRI_H1(order))
497  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
498  "Inconsistent number of DOFs");
499 
500  auto &tri_alpha = get_alpha(ent_data);
501  tri_alpha.resize(nb_dofs, 3, false);
502 
504  &tri_alpha(0, 0));
505 
506  if (order > 1) {
507  std::array<int, 3> face_n{order, order, order};
508  std::array<int *, 3> face_edge_ptr{
509  &tri_alpha(3, 0), &tri_alpha(3 + NBEDGE_H1(order), 0),
510  &tri_alpha(3 + 2 * NBEDGE_H1(order), 0)};
512  face_n.data(), face_edge_ptr.data());
513  if (order > 2)
515  order, &tri_alpha(3 + 3 * NBEDGE_H1(order), 0));
516  }
518  order, lambda.size1(), tri_alpha.size1(), &tri_alpha(0, 0),
519  &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &get_n(0, 0),
520  &get_diff_n(0, 0));
521  }
522 
523  get_alpha_by_order_ptr(ent_data, order) =
524  get_alpha_by_name_ptr(ent_data, field_name);
525  get_base_by_order_ptr(ent_data, order) =
526  get_base_by_name_ptr(ent_data, field_name);
527  get_diff_base_by_order_ptr(ent_data, order) =
528  get_diff_base_by_name_ptr(ent_data, field_name);
529  }
530  }
531  } else {
532  auto &ent_data = data.dataOnEntities[MBTRI][0];
533  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
534  auto &get_n = get_base(ent_data);
535  auto &get_diff_n = get_diff_base(ent_data);
536  get_n.resize(nb_gauss_pts, 0, false);
537  get_diff_n.resize(nb_gauss_pts, 0, false);
538  }
539 
541 }
542 
545 
547  const FieldApproximationBase base = cTx->bAse;
548  if (cTx->basePolynomialsType0 == NULL)
549  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
550  "Polynomial type not set");
551  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
552  double *diffL, const int dim) =
554 
555  int nb_gauss_pts = pts.size2();
556 
557  double *PHI_f_e[3];
558  double *PHI_f;
559 
560  N_face_edge.resize(1, 3, false);
561  N_face_bubble.resize(1, false);
562  int face_order = data.dataOnEntities[MBTRI][0].getDataOrder();
563  // three edges on face
564  for (int ee = 0; ee < 3; ee++) {
565  N_face_edge(0, ee).resize(
566  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(face_order), false);
567  PHI_f_e[ee] = &((N_face_edge(0, ee))(0, 0));
568  }
569  N_face_bubble[0].resize(nb_gauss_pts,
570  3 * NBFACETRI_AINSWORTH_FACE_HDIV(face_order), false);
571  PHI_f = &*(N_face_bubble[0].data().begin());
572 
573  int face_nodes[3] = {0, 1, 2};
575  face_nodes, face_order,
576  &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0), NULL, PHI_f_e, NULL,
577  nb_gauss_pts, 3, base_polynomials);
579  face_nodes, face_order,
580  &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0), NULL, PHI_f, NULL,
581  nb_gauss_pts, 3, base_polynomials);
582 
583  // set shape functions into data structure
584  if (data.dataOnEntities[MBTRI].size() != 1) {
585  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
586  }
587  data.dataOnEntities[MBTRI][0].getN(base).resize(
588  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_HDIV(face_order), false);
589  int col = 0;
590  for (int oo = 0; oo < face_order; oo++) {
591  for (int ee = 0; ee < 3; ee++) {
592  for (int dd = 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
593  dd < 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); 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 = 3 * NBFACETRI_AINSWORTH_FACE_HDIV(oo);
601  dd < 3 * NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++, col++) {
602  for (int gg = 0; gg < nb_gauss_pts; gg++) {
603  data.dataOnEntities[MBTRI][0].getN(base)(gg, col) =
604  N_face_bubble[0](gg, dd);
605  }
606  }
607  }
608 
610 }
611 
614 
616  // set shape functions into data structure
617  if (data.dataOnEntities[MBTRI].size() != 1) {
618  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
619  }
620  const FieldApproximationBase base = cTx->bAse;
621  if (base != DEMKOWICZ_JACOBI_BASE) {
622  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
623  "This should be used only with DEMKOWICZ_JACOBI_BASE "
624  "but base is %s",
625  ApproximationBaseNames[base]);
626  }
627  int order = data.dataOnEntities[MBTRI][0].getDataOrder();
628  int nb_gauss_pts = pts.size2();
629  data.dataOnEntities[MBTRI][0].getN(base).resize(
630  nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
631  double *phi_f = &*data.dataOnEntities[MBTRI][0].getN(base).data().begin();
634  int face_nodes[3] = {0, 1, 2};
636  face_nodes, order,
637  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
638  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(), phi_f,
639  NULL, nb_gauss_pts, 3);
640 
642 }
643 
646 
647  switch (cTx->bAse) {
650  return getValueHdivAinsworthBase(pts);
652  return getValueHdivDemkowiczBase(pts);
653  default:
654  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
655  }
656 
658 }
659 
663 
665  const FieldApproximationBase base = cTx->bAse;
666  if (data.dataOnEntities[MBTRI].size() != 1)
667  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
668 
669  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
670  double *diffL, const int dim) =
672 
673  int nb_gauss_pts = pts.size2();
674 
675  // Calculation H-curl on triangle faces
676  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
677  if (data.dataOnEntities[MBEDGE].size() != 3) {
678  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
679  }
680  int sense[3], order[3];
681  double *hcurl_edge_n[3];
682  double *diff_hcurl_edge_n[3];
683  for (int ee = 0; ee < 3; ee++) {
684  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
685  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
686  "data inconsistency");
687  }
688  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
689  order[ee] = data.dataOnEntities[MBEDGE][ee].getDataOrder();
690  int nb_dofs = NBEDGE_AINSWORTH_HCURL(
691  data.dataOnEntities[MBEDGE][ee].getDataOrder());
692  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
693  3 * nb_dofs, false);
694  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(
695  nb_gauss_pts, 2 * 3 * nb_dofs, false);
696  hcurl_edge_n[ee] =
697  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
698  diff_hcurl_edge_n[ee] =
699  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
700  }
702  sense, order,
703  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
704  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
705  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
706  } else {
707  for (int ee = 0; ee < 3; ee++) {
708  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
709  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
710  false);
711  }
712  }
713 
714  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
715 
716  // cerr << data.dataOnEntities[MBVERTEX][0].getN(base) << endl;
717  // cerr << data.dataOnEntities[MBVERTEX][0].getDiffN(base) << endl;
718  //
719  // face
720  if (data.dataOnEntities[MBTRI].size() != 1) {
721  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
722  }
723  int order = data.dataOnEntities[MBTRI][0].getDataOrder();
724  int nb_dofs = NBFACETRI_AINSWORTH_HCURL(order);
725  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
726  false);
727  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
728  3 * 2 * nb_dofs, false);
729  // cerr << data.dataOnEntities[MBVERTEX][0].getDiffN(base) << endl;
730  int face_nodes[] = {0, 1, 2};
732  face_nodes, order,
733  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
734  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
735  &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
736  &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
737  nb_gauss_pts, base_polynomials);
738  // cerr << data.dataOnEntities[MBTRI][0].getN(base) << endl;
739  } else {
740  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 0, false);
741  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
742  }
743 
745 }
746 
750 
752  const FieldApproximationBase base = cTx->bAse;
753 
754  int nb_gauss_pts = pts.size2();
755 
756  // Calculation H-curl on triangle faces
757  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
758 
759  if (data.dataOnEntities[MBEDGE].size() != 3)
760  SETERRQ1(
761  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
762  "wrong number of data structures on edges, should be three but is %d",
763  data.dataOnEntities[MBEDGE].size());
764 
765  int sense[3], order[3];
766  double *hcurl_edge_n[3];
767  double *diff_hcurl_edge_n[3];
768 
769  for (int ee = 0; ee != 3; ++ee) {
770 
771  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0)
772  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
773  "orientation (sense) of edge is not set");
774 
775  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
776  order[ee] = data.dataOnEntities[MBEDGE][ee].getDataOrder();
777  int nb_dofs = NBEDGE_DEMKOWICZ_HCURL(
778  data.dataOnEntities[MBEDGE][ee].getDataOrder());
779  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
780  3 * nb_dofs, false);
781  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(
782  nb_gauss_pts, 2 * 3 * nb_dofs, false);
783 
784  hcurl_edge_n[ee] =
785  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
786  diff_hcurl_edge_n[ee] =
787  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
788  }
789 
791  sense, order,
792  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
793  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
794  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts);
795 
796  } else {
797 
798  // No DOFs on faces, resize base function matrices, indicating that no
799  // dofs on them.
800  for (int ee = 0; ee != 3; ++ee) {
801  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
802  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
803  false);
804  }
805  }
806 
807  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
808 
809  // face
810  if (data.dataOnEntities[MBTRI].size() != 1)
811  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
812  "No data struture to keep base functions on face");
813 
814  int order = data.dataOnEntities[MBTRI][0].getDataOrder();
815  int nb_dofs = NBFACETRI_DEMKOWICZ_HCURL(order);
816  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
817  false);
818  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
819  3 * 2 * nb_dofs, false);
820 
821  int face_nodes[] = {0, 1, 2};
823  face_nodes, order,
824  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
825  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
826  &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
827  &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
828  nb_gauss_pts);
829 
830  } else {
831 
832  // No DOFs on faces, resize base function matrices, indicating that no
833  // dofs on them.
834  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 0, false);
835  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
836  }
837 
839 }
840 
843 
844  switch (cTx->bAse) {
848  break;
851  break;
852  default:
853  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
854  }
855 
857 }
858 
861  boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
863 
865  CHKERR ctx_ptr->query_interface(IDD_TRI_BASE_FUNCTION, &iface);
866  cTx = reinterpret_cast<EntPolynomialBaseCtx *>(iface);
867 
868  int nb_gauss_pts = pts.size2();
869  if (!nb_gauss_pts) {
871  }
872 
873  if (pts.size1() < 2)
874  SETERRQ(
875  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
876  "Wrong dimension of pts, should be at least 3 rows with coordinates");
877 
878  const FieldApproximationBase base = cTx->bAse;
880 
881  if (base != AINSWORTH_BERNSTEIN_BEZIER_BASE) {
882  if (cTx->copyNodeBase == LASTBASE) {
883  data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 3,
884  false);
886  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
887  &pts(0, 0), &pts(1, 0), nb_gauss_pts);
888  data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(3, 2, false);
890  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
891  } else {
892  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
893  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
894  data.dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(base) =
895  data.dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(cTx->copyNodeBase);
896  }
897  if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
898  (unsigned int)nb_gauss_pts) {
899  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
900  "Base functions or nodes has wrong number of integration points "
901  "for base %s",
902  ApproximationBaseNames[base]);
903  }
904  }
905  auto set_node_derivative_for_all_gauss_pts = [&]() {
907  // In linear geometry derivatives are constant,
908  // this in expense of efficiency makes implementation
909  // consistent between vertices and other types of entities
910  data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(nb_gauss_pts, 6,
911  false);
912  for (int gg = 0; gg != nb_gauss_pts; ++gg)
913  std::copy(Tools::diffShapeFunMBTRI.begin(),
915  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(gg, 0));
917  };
918 
919  CHKERR set_node_derivative_for_all_gauss_pts();
920 
921  switch (cTx->sPace) {
922  case H1:
923  CHKERR getValueH1(pts);
924  break;
925  case HDIV:
926  CHKERR getValueHdiv(pts);
927  break;
928  case HCURL:
929  CHKERR getValueHcurl(pts);
930  break;
931  case L2:
932  CHKERR getValueL2(pts);
933  break;
934  default:
935  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
936  }
937 
939 }
MoFEMErrorCode getValueL2AinsworthBase(MatrixDouble &pts)
MoFEMErrorCode Hcurl_Ainsworth_FaceFunctions_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face H-curl functions.
Definition: Hcurl.cpp:1249
MoFEM interface unique ID.
MoFEMErrorCode getValueHcurlAinsworthBase(MatrixDouble &pts)
field with continuous normal traction
Definition: definitions.h:179
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
ublas::vector< MatrixDouble > N_face_bubble
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
data structure for finite element entityIt keeps that about indices of degrees of freedom,...
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:27
MoFEMErrorCode Hdiv_Demkowicz_Face_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb)
Definition: Hdiv.cpp:617
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
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
static MoFEMErrorCode generateIndicesEdgeTri(const int N[], int *alpha[])
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, MoFEM::BaseFunctionUnknownInterface **iface) const
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
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
const FieldApproximationBase bAse
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, BaseFunctionUnknownInterface **iface) const
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:162
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
ublas::matrix< int, ublas::row_major, IntAllocator > MatrixInt
Definition: Types.hpp:74
static const MOFEMuuid IDD_TET_BASE_FUNCTION
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEMErrorCode getValueH1(MatrixDouble &pts)
const FieldApproximationBase copyNodeBase
static const char *const ApproximationBaseNames[]
Definition: definitions.h:164
const int dim
EntPolynomialBaseCtx * cTx
static MoFEMErrorCode generateIndicesVertexTri(const int N, int *alpha)
MoFEMErrorCode getValueL2(MatrixDouble &pts)
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:206
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
#define NBFACETRI_DEMKOWICZ_HDIV(P)
FieldApproximationBase
approximation base
Definition: definitions.h:150
DataForcesAndSourcesCore & dAta
#define NBEDGE_AINSWORTH_HCURL(P)
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
#define NBFACETRI_L2(P)
Number of base functions on triangle for L2 space.
field with continuous tangents
Definition: definitions.h:178
MoFEMErrorCode getValueL2BernsteinBezierBase(MatrixDouble &pts)
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
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)
ublas::matrix< MatrixDouble > N_face_edge
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:87
#define CHKERR
Inline error check.
Definition: definitions.h:602
#define NBFACETRI_DEMKOWICZ_HCURL(P)
MoFEMErrorCode Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f_e[3], double *diff_phi_f_e[3], int gdim, int nb, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base functions, Edge-based face functions by Ainsworth .
Definition: Hdiv.cpp:37
#define NBEDGE_DEMKOWICZ_HCURL(P)
constexpr int order
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition: fem_tools.c:194
MoFEMErrorCode Hdiv_Ainsworth_FaceBubbleShapeFunctions_ON_FACE(int *faces_nodes, int p, double *N, double *diffN, double *phi_f, double *diff_phi_f, int gdim, int nb, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Face bubble functions by Ainsworth .
Definition: Hdiv.cpp:160
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)
#define NBFACETRI_AINSWORTH_FACE_HDIV(P)
MoFEMErrorCode getValueHdivAinsworthBase(MatrixDouble &pts)
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
#define NBFACETRI_AINSWORTH_HCURL(P)
static const MOFEMuuid IDD_TRI_BASE_FUNCTION
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
continuous field
Definition: definitions.h:177
#define NBFACETRI_AINSWORTH_HDIV(P)
MoFEMErrorCode getValueH1BernsteinBezierBase(MatrixDouble &pts)
static MoFEMErrorCode generateIndicesTriTri(const int N, int *alpha)
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
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:29
MoFEMErrorCode getValueH1AinsworthBase(MatrixDouble &pts)
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26
field with C-1 continuity
Definition: definitions.h:180