v0.13.0
Public Member Functions | Private Member Functions | Private Attributes | List of all members
MoFEM::TriPolynomialBase Struct Reference

Calculate base functions on triangle. More...

#include <src/approximation/TriPolynomialBase.hpp>

Inheritance diagram for MoFEM::TriPolynomialBase:
[legend]
Collaboration diagram for MoFEM::TriPolynomialBase:
[legend]

Public Member Functions

MoFEMErrorCode query_interface (boost::typeindex::type_index type_index, UnknownInterface **iface) const
 
 TriPolynomialBase ()
 
 ~TriPolynomialBase ()
 
MoFEMErrorCode getValue (MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
 
- Public Member Functions inherited from MoFEM::BaseFunction
virtual MoFEMErrorCode getValue (MatrixDouble &pts_x, MatrixDouble &pts_t, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
 
- Public Member Functions inherited from MoFEM::BaseFunctionUnknownInterface
virtual ~BaseFunctionUnknownInterface ()=default
 
- Public Member Functions inherited from MoFEM::UnknownInterface
template<class IFACE >
MoFEMErrorCode registerInterface (bool error_if_registration_failed=true)
 Register interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE *&iface) const
 Get interface refernce to pointer of interface. More...
 
template<class IFACE >
MoFEMErrorCode getInterface (IFACE **const iface) const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_pointer< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get interface pointer to pointer of interface. More...
 
template<class IFACE , typename boost::enable_if< boost::is_reference< IFACE >, int >::type = 0>
IFACE getInterface () const
 Get reference to interface. More...
 
template<class IFACE >
IFACE * getInterface () const
 Function returning pointer to interface. More...
 
virtual ~UnknownInterface ()=default
 

Private Member Functions

MoFEMErrorCode getValueH1 (MatrixDouble &pts)
 
MoFEMErrorCode getValueH1AinsworthBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueH1BernsteinBezierBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueL2 (MatrixDouble &pts)
 
MoFEMErrorCode getValueL2AinsworthBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueL2BernsteinBezierBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHdiv (MatrixDouble &pts)
 
MoFEMErrorCode getValueHcurl (MatrixDouble &pts)
 
MoFEMErrorCode getValueHdivAinsworthBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHcurlAinsworthBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHdivDemkowiczBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueHcurlDemkowiczBase (MatrixDouble &pts)
 

Private Attributes

EntPolynomialBaseCtxcTx
 
ublas::matrix< MatrixDoubleN_face_edge
 
ublas::vector< MatrixDoubleN_face_bubble
 
ublas::matrix< MatrixDoublediffN_face_edge
 
ublas::vector< MatrixDoublediffN_face_bubble
 

Additional Inherited Members

- Static Public Member Functions inherited from MoFEM::UnknownInterface
static MoFEMErrorCode getLibVersion (Version &version)
 Get library version. More...
 
static MoFEMErrorCode getFileVersion (moab::Interface &moab, Version &version)
 Get database major version. More...
 
static MoFEMErrorCode setFileVersion (moab::Interface &moab, Version version=Version(MoFEM_VERSION_MAJOR, MoFEM_VERSION_MINOR, MoFEM_VERSION_BUILD))
 Get database major version. More...
 
static MoFEMErrorCode getInterfaceVersion (Version &version)
 Get database major version. More...
 

Detailed Description

Calculate base functions on triangle.

Definition at line 30 of file TriPolynomialBase.hpp.

Constructor & Destructor Documentation

◆ TriPolynomialBase()

TriPolynomialBase::TriPolynomialBase ( )

Definition at line 22 of file TriPolynomialBase.cpp.

22 {}

◆ ~TriPolynomialBase()

TriPolynomialBase::~TriPolynomialBase ( )

Definition at line 23 of file TriPolynomialBase.cpp.

23 {}

Member Function Documentation

◆ getValue()

MoFEMErrorCode TriPolynomialBase::getValue ( MatrixDouble pts,
boost::shared_ptr< BaseFunctionCtx ctx_ptr 
)
virtual

Reimplemented from MoFEM::BaseFunction.

Definition at line 859 of file TriPolynomialBase.cpp.

860  {
862 
863  cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
864 
865  int nb_gauss_pts = pts.size2();
866  if (!nb_gauss_pts) {
868  }
869 
870  if (pts.size1() < 2)
871  SETERRQ(
872  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
873  "Wrong dimension of pts, should be at least 3 rows with coordinates");
874 
875  const FieldApproximationBase base = cTx->bAse;
876  EntitiesFieldData &data = cTx->dAta;
877 
878  if (base != AINSWORTH_BERNSTEIN_BEZIER_BASE) {
879  if (cTx->copyNodeBase == LASTBASE) {
880  data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 3,
881  false);
883  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
884  &pts(0, 0), &pts(1, 0), nb_gauss_pts);
885  data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(3, 2, false);
887  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
888  } else {
889  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
890  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
891  data.dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(base) =
892  data.dataOnEntities[MBVERTEX][0].getDiffNSharedPtr(cTx->copyNodeBase);
893  }
894  if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
895  (unsigned int)nb_gauss_pts) {
896  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
897  "Base functions or nodes has wrong number of integration points "
898  "for base %s",
899  ApproximationBaseNames[base]);
900  }
901  }
902  auto set_node_derivative_for_all_gauss_pts = [&]() {
904  // In linear geometry derivatives are constant,
905  // this in expense of efficiency makes implementation
906  // consistent between vertices and other types of entities
907  data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(nb_gauss_pts, 6,
908  false);
909  for (int gg = 0; gg != nb_gauss_pts; ++gg)
910  std::copy(Tools::diffShapeFunMBTRI.begin(),
912  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(gg, 0));
914  };
915 
916  CHKERR set_node_derivative_for_all_gauss_pts();
917 
918  switch (cTx->sPace) {
919  case H1:
920  CHKERR getValueH1(pts);
921  break;
922  case HDIV:
923  CHKERR getValueHdiv(pts);
924  break;
925  case HCURL:
926  CHKERR getValueHcurl(pts);
927  break;
928  case L2:
929  CHKERR getValueL2(pts);
930  break;
931  default:
932  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
933  }
934 
936 }
FieldApproximationBase
approximation base
Definition: definitions.h:71
@ LASTBASE
Definition: definitions.h:82
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:77
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ H1
continuous field
Definition: definitions.h:98
@ HCURL
field with continuous tangents
Definition: definitions.h:99
@ HDIV
field with continuous normal traction
Definition: definitions.h:100
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:45
static const char *const ApproximationBaseNames[]
Definition: definitions.h:85
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:206
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition: fem_tools.c:194
Class used to pass element data to calculate base functions on tet,triangle,edge.
const FieldApproximationBase copyNodeBase
const FieldApproximationBase bAse
data structure for finite element entity
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:83
MoFEMErrorCode getValueH1(MatrixDouble &pts)
MoFEMErrorCode getValueL2(MatrixDouble &pts)
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
EntPolynomialBaseCtx * cTx
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

◆ getValueH1()

MoFEMErrorCode TriPolynomialBase::getValueH1 ( MatrixDouble pts)
private

Definition at line 33 of file TriPolynomialBase.cpp.

33  {
35 
36  switch (cTx->bAse) {
40  break;
43  break;
44  default:
45  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
46  }
47 
49 }
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:75
MoFEMErrorCode getValueH1BernsteinBezierBase(MatrixDouble &pts)
MoFEMErrorCode getValueH1AinsworthBase(MatrixDouble &pts)

◆ getValueH1AinsworthBase()

MoFEMErrorCode TriPolynomialBase::getValueH1AinsworthBase ( MatrixDouble pts)
private

Definition at line 265 of file TriPolynomialBase.cpp.

265  {
267 
268  EntitiesFieldData &data = cTx->dAta;
269  const FieldApproximationBase base = cTx->bAse;
270 
271  if (cTx->basePolynomialsType0 == NULL)
272  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
273  "Polynomial type not set");
274 
275  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
276  double *diffL, const int dim) =
278 
279  int nb_gauss_pts = pts.size2();
280 
281  if (data.spacesOnEntities[MBEDGE].test(H1)) {
282  // edges
283  if (data.dataOnEntities[MBEDGE].size() != 3)
284  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
285  "expected size of data.dataOnEntities[MBEDGE] is 3");
286 
287  int sense[3], order[3];
288  double *H1edgeN[3], *diffH1edgeN[3];
289  for (int ee = 0; ee < 3; ee++) {
290 
291  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0)
292  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
293  "sense of the edge unknown");
294 
295  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
296  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
297  int nb_dofs = NBEDGE_H1(data.dataOnEntities[MBEDGE][ee].getOrder());
298  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
299  false);
300  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
301  2 * nb_dofs, false);
302  H1edgeN[ee] = &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
303  diffH1edgeN[ee] =
304  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
305  }
307  sense, order,
308  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
309  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
310  H1edgeN, diffH1edgeN, nb_gauss_pts, base_polynomials);
311  }
312 
313  if (data.spacesOnEntities[MBTRI].test(H1)) {
314  // face
315  if (data.dataOnEntities[MBTRI].size() != 1) {
316  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
317  "expected that size data.dataOnEntities[MBTRI] is one");
318  }
319 
320  int nb_dofs = NBFACETRI_H1(data.dataOnEntities[MBTRI][0].getOrder());
321  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, nb_dofs,
322  false);
323  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
324  2 * nb_dofs, false);
325  const int face_nodes[] = {0, 1, 2};
327  face_nodes, data.dataOnEntities[MBTRI][0].getOrder(),
328  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
329  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
330  &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
331  &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
332  nb_gauss_pts, base_polynomials);
333  }
334 
336 }
static Index< 'p', 3 > p
const int dim
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
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
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:201
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
PetscErrorCode(* basePolynomialsType0)(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types

◆ getValueH1BernsteinBezierBase()

MoFEMErrorCode TriPolynomialBase::getValueH1BernsteinBezierBase ( MatrixDouble pts)
private

Definition at line 52 of file TriPolynomialBase.cpp.

52  {
54  EntitiesFieldData &data = cTx->dAta;
55  const std::string &field_name = cTx->fieldName;
56  int nb_gauss_pts = pts.size2();
57 
58  auto get_alpha = [field_name](auto &data) -> MatrixInt & {
59  auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
60  if (!ptr)
61  ptr.reset(new MatrixInt());
62  return *ptr;
63  };
64 
65  auto get_base = [field_name](auto &data) -> MatrixDouble & {
66  auto &ptr = data.getBBNSharedPtr(field_name);
67  if (!ptr)
68  ptr.reset(new MatrixDouble());
69  return *ptr;
70  };
71 
72  auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
73  auto &ptr = data.getBBDiffNSharedPtr(field_name);
74  if (!ptr)
75  ptr.reset(new MatrixDouble());
76  return *ptr;
77  };
78 
79  auto get_alpha_by_name_ptr =
80  [](auto &data,
81  const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
82  return data.getBBAlphaIndicesSharedPtr(field_name);
83  };
84 
85  auto get_base_by_name_ptr =
86  [](auto &data,
87  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
88  return data.getBBNSharedPtr(field_name);
89  };
90 
91  auto get_diff_base_by_name_ptr =
92  [](auto &data,
93  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
94  return data.getBBDiffNSharedPtr(field_name);
95  };
96 
97  auto get_alpha_by_order_ptr =
98  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
99  return data.getBBAlphaIndicesByOrderSharedPtr(o);
100  };
101 
102  auto get_base_by_order_ptr =
103  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
104  return data.getBBNByOrderSharedPtr(o);
105  };
106 
107  auto get_diff_base_by_order_ptr =
108  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
109  return data.getBBDiffNByOrderSharedPtr(o);
110  };
111 
112  auto &vert_get_n = get_base(data.dataOnEntities[MBVERTEX][0]);
113  auto &vert_get_diff_n = get_diff_base(data.dataOnEntities[MBVERTEX][0]);
114  vert_get_n.resize(nb_gauss_pts, 3, false);
115  vert_get_diff_n.resize(nb_gauss_pts, 6, false);
116 
117  if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
118  (unsigned int)nb_gauss_pts)
119  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
120  "Base functions or nodes has wrong number of integration points "
121  "for base %s",
123  auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
124 
125  auto &vertex_alpha = get_alpha(data.dataOnEntities[MBVERTEX][0]);
126  vertex_alpha.resize(3, 3, false);
127  vertex_alpha.clear();
128  for (int n = 0; n != 3; ++n)
129  vertex_alpha(n, n) = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n];
130 
132  1, lambda.size1(), vertex_alpha.size1(), &vertex_alpha(0, 0),
133  &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &vert_get_n(0, 0),
134  &vert_get_diff_n(0, 0));
135  for (int n = 0; n != 3; ++n) {
136  const int f = boost::math::factorial<double>(
137  data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n]);
138  for (int g = 0; g != nb_gauss_pts; ++g) {
139  vert_get_n(g, n) *= f;
140  for (int d = 0; d != 2; ++d)
141  vert_get_diff_n(g, 2 * n + d) *= f;
142  }
143  }
144 
145  // edges
146  if (data.spacesOnEntities[MBEDGE].test(H1)) {
147  if (data.dataOnEntities[MBEDGE].size() != 3)
148  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
149  "Wrong size of ent data");
150 
151  constexpr int edges_nodes[3][2] = {{0, 1}, {1, 2}, {2, 0}};
152  for (int ee = 0; ee != 3; ++ee) {
153  auto &ent_data = data.dataOnEntities[MBEDGE][ee];
154 
155  if (ent_data.getSense() == 0)
156  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
157  "Sense of the edge unknown");
158 
159  const int sense = ent_data.getSense();
160  const int order = ent_data.getOrder();
161  const int nb_dofs = NBEDGE_H1(ent_data.getOrder());
162 
163  if (nb_dofs) {
164  if (get_alpha_by_order_ptr(ent_data, order)) {
165  get_alpha_by_name_ptr(ent_data, field_name) =
166  get_alpha_by_order_ptr(ent_data, order);
167  get_base_by_name_ptr(ent_data, field_name) =
168  get_base_by_order_ptr(ent_data, order);
169  get_diff_base_by_name_ptr(ent_data, field_name) =
170  get_diff_base_by_order_ptr(ent_data, order);
171  } else {
172  auto &get_n = get_base(ent_data);
173  auto &get_diff_n = get_diff_base(ent_data);
174  get_n.resize(nb_gauss_pts, nb_dofs, false);
175  get_diff_n.resize(nb_gauss_pts, 2 * nb_dofs, false);
176 
177  auto &edge_alpha = get_alpha(data.dataOnEntities[MBEDGE][ee]);
178  edge_alpha.resize(nb_dofs, 3, false);
180  &edge_alpha(0, 0));
181  if (sense == -1) {
182  for (int i = 0; i != edge_alpha.size1(); ++i) {
183  int a = edge_alpha(i, edges_nodes[ee][0]);
184  edge_alpha(i, edges_nodes[ee][0]) =
185  edge_alpha(i, edges_nodes[ee][1]);
186  edge_alpha(i, edges_nodes[ee][1]) = a;
187  }
188  }
190  order, lambda.size1(), edge_alpha.size1(), &edge_alpha(0, 0),
191  &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &get_n(0, 0),
192  &get_diff_n(0, 0));
193 
194  get_alpha_by_order_ptr(ent_data, order) =
195  get_alpha_by_name_ptr(ent_data, field_name);
196  get_base_by_order_ptr(ent_data, order) =
197  get_base_by_name_ptr(ent_data, field_name);
198  get_diff_base_by_order_ptr(ent_data, order) =
199  get_diff_base_by_name_ptr(ent_data, field_name);
200  }
201  }
202  }
203  } else {
204  for (int ee = 0; ee != 3; ++ee) {
205  auto &ent_data = data.dataOnEntities[MBEDGE][ee];
206  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
207  auto &get_n = get_base(ent_data);
208  auto &get_diff_n = get_diff_base(ent_data);
209  get_n.resize(nb_gauss_pts, 0, false);
210  get_diff_n.resize(nb_gauss_pts, 0, false);
211  }
212  }
213 
214  if (data.spacesOnEntities[MBTRI].test(H1)) {
215  if (data.dataOnEntities[MBTRI].size() != 1)
216  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
217  "Wrong size ent of ent data");
218 
219  auto &ent_data = data.dataOnEntities[MBTRI][0];
220  const int order = ent_data.getOrder();
221  const int nb_dofs = NBFACETRI_H1(order);
222  if (get_alpha_by_order_ptr(ent_data, order)) {
223  get_alpha_by_name_ptr(ent_data, field_name) =
224  get_alpha_by_order_ptr(ent_data, order);
225  get_base_by_name_ptr(ent_data, field_name) =
226  get_base_by_order_ptr(ent_data, order);
227  get_diff_base_by_name_ptr(ent_data, field_name) =
228  get_diff_base_by_order_ptr(ent_data, order);
229  } else {
230 
231  auto &get_n = get_base(ent_data);
232  auto &get_diff_n = get_diff_base(ent_data);
233  get_n.resize(nb_gauss_pts, nb_dofs, false);
234  get_diff_n.resize(nb_gauss_pts, 2 * nb_dofs, false);
235  if (nb_dofs) {
236  auto &tri_alpha = get_alpha(ent_data);
237  tri_alpha.resize(nb_dofs, 3, false);
238 
241  order, lambda.size1(), tri_alpha.size1(), &tri_alpha(0, 0),
242  &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &get_n(0, 0),
243  &get_diff_n(0, 0));
244 
245  get_alpha_by_order_ptr(ent_data, order) =
246  get_alpha_by_name_ptr(ent_data, field_name);
247  get_base_by_order_ptr(ent_data, order) =
248  get_base_by_name_ptr(ent_data, field_name);
249  get_diff_base_by_order_ptr(ent_data, order) =
250  get_diff_base_by_name_ptr(ent_data, field_name);
251  }
252  }
253  } else {
254  auto &ent_data = data.dataOnEntities[MBTRI][0];
255  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
256  auto &get_n = get_base(ent_data);
257  auto &get_diff_n = get_diff_base(ent_data);
258  get_n.resize(nb_gauss_pts, 0, false);
259  get_diff_n.resize(nb_gauss_pts, 0, false);
260  }
261 
263 }
static Index< 'n', 3 > n
static Index< 'o', 3 > o
constexpr double a
@ NOBASE
Definition: definitions.h:72
constexpr double lambda
FTensor::Index< 'i', SPACE_DIM > i
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
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasMatrix< int > MatrixInt
Definition: Types.hpp:87
constexpr double g
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)
static MoFEMErrorCode generateIndicesTriTri(const int N, int *alpha)
static MoFEMErrorCode generateIndicesEdgeTri(const int N[], int *alpha[])

◆ getValueHcurl()

MoFEMErrorCode TriPolynomialBase::getValueHcurl ( MatrixDouble pts)
private

Definition at line 840 of file TriPolynomialBase.cpp.

840  {
842 
843  switch (cTx->bAse) {
847  break;
850  break;
851  default:
852  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
853  }
854 
856 }
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)
MoFEMErrorCode getValueHcurlAinsworthBase(MatrixDouble &pts)

◆ getValueHcurlAinsworthBase()

MoFEMErrorCode TriPolynomialBase::getValueHcurlAinsworthBase ( MatrixDouble pts)
private

Definition at line 660 of file TriPolynomialBase.cpp.

660  {
662 
663  EntitiesFieldData &data = cTx->dAta;
664  const FieldApproximationBase base = cTx->bAse;
665  if (data.dataOnEntities[MBTRI].size() != 1)
666  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
667 
668  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
669  double *diffL, const int dim) =
671 
672  int nb_gauss_pts = pts.size2();
673 
674  // Calculation H-curl on triangle faces
675  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
676  if (data.dataOnEntities[MBEDGE].size() != 3) {
677  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
678  }
679  int sense[3], order[3];
680  double *hcurl_edge_n[3];
681  double *diff_hcurl_edge_n[3];
682  for (int ee = 0; ee < 3; ee++) {
683  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
684  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
685  "data inconsistency");
686  }
687  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
688  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
689  int nb_dofs = NBEDGE_AINSWORTH_HCURL(
690  data.dataOnEntities[MBEDGE][ee].getOrder());
691  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
692  3 * nb_dofs, false);
693  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(
694  nb_gauss_pts, 2 * 3 * nb_dofs, false);
695  hcurl_edge_n[ee] =
696  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
697  diff_hcurl_edge_n[ee] =
698  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
699  }
701  sense, order,
702  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
703  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
704  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
705  } else {
706  for (int ee = 0; ee < 3; ee++) {
707  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
708  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
709  false);
710  }
711  }
712 
713  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
714 
715  // cerr << data.dataOnEntities[MBVERTEX][0].getN(base) << endl;
716  // cerr << data.dataOnEntities[MBVERTEX][0].getDiffN(base) << endl;
717  //
718  // face
719  if (data.dataOnEntities[MBTRI].size() != 1) {
720  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
721  }
722  int order = data.dataOnEntities[MBTRI][0].getOrder();
723  int nb_dofs = NBFACETRI_AINSWORTH_HCURL(order);
724  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
725  false);
726  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
727  3 * 2 * nb_dofs, false);
728  // cerr << data.dataOnEntities[MBVERTEX][0].getDiffN(base) << endl;
729  int face_nodes[] = {0, 1, 2};
731  face_nodes, order,
732  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
733  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
734  &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
735  &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
736  nb_gauss_pts, base_polynomials);
737  // cerr << data.dataOnEntities[MBTRI][0].getN(base) << endl;
738  } else {
739  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 0, false);
740  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
741  }
742 
744 }
#define NBFACETRI_AINSWORTH_HCURL(P)
#define NBEDGE_AINSWORTH_HCURL(P)
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
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

◆ getValueHcurlDemkowiczBase()

MoFEMErrorCode TriPolynomialBase::getValueHcurlDemkowiczBase ( MatrixDouble pts)
private

Definition at line 747 of file TriPolynomialBase.cpp.

747  {
749 
750  EntitiesFieldData &data = cTx->dAta;
751  const FieldApproximationBase base = cTx->bAse;
752 
753  int nb_gauss_pts = pts.size2();
754 
755  // Calculation H-curl on triangle faces
756  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
757 
758  if (data.dataOnEntities[MBEDGE].size() != 3)
759  SETERRQ1(
760  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
761  "wrong number of data structures on edges, should be three but is %d",
762  data.dataOnEntities[MBEDGE].size());
763 
764  int sense[3], order[3];
765  double *hcurl_edge_n[3];
766  double *diff_hcurl_edge_n[3];
767 
768  for (int ee = 0; ee != 3; ++ee) {
769 
770  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0)
771  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
772  "orientation (sense) of edge is not set");
773 
774  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
775  order[ee] = data.dataOnEntities[MBEDGE][ee].getOrder();
776  int nb_dofs = NBEDGE_DEMKOWICZ_HCURL(
777  data.dataOnEntities[MBEDGE][ee].getOrder());
778  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
779  3 * nb_dofs, false);
780  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(
781  nb_gauss_pts, 2 * 3 * nb_dofs, false);
782 
783  hcurl_edge_n[ee] =
784  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
785  diff_hcurl_edge_n[ee] =
786  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
787  }
788 
790  sense, order,
791  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
792  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
793  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts);
794 
795  } else {
796 
797  // No DOFs on faces, resize base function matrices, indicating that no
798  // dofs on them.
799  for (int ee = 0; ee != 3; ++ee) {
800  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
801  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
802  false);
803  }
804  }
805 
806  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
807 
808  // face
809  if (data.dataOnEntities[MBTRI].size() != 1)
810  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
811  "No data struture to keep base functions on face");
812 
813  int order = data.dataOnEntities[MBTRI][0].getOrder();
814  int nb_dofs = NBFACETRI_DEMKOWICZ_HCURL(order);
815  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 3 * nb_dofs,
816  false);
817  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts,
818  3 * 2 * nb_dofs, false);
819 
820  int face_nodes[] = {0, 1, 2};
822  face_nodes, order,
823  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
824  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
825  &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
826  &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
827  nb_gauss_pts);
828 
829  } else {
830 
831  // No DOFs on faces, resize base function matrices, indicating that no
832  // dofs on them.
833  data.dataOnEntities[MBTRI][0].getN(base).resize(nb_gauss_pts, 0, false);
834  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
835  }
836 
838 }
#define NBEDGE_DEMKOWICZ_HCURL(P)
#define NBFACETRI_DEMKOWICZ_HCURL(P)
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
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

◆ getValueHdiv()

MoFEMErrorCode TriPolynomialBase::getValueHdiv ( MatrixDouble pts)
private

Definition at line 643 of file TriPolynomialBase.cpp.

643  {
645 
646  switch (cTx->bAse) {
649  return getValueHdivAinsworthBase(pts);
651  return getValueHdivDemkowiczBase(pts);
652  default:
653  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
654  }
655 
657 }
MoFEMErrorCode getValueHdivAinsworthBase(MatrixDouble &pts)
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)

◆ getValueHdivAinsworthBase()

MoFEMErrorCode TriPolynomialBase::getValueHdivAinsworthBase ( MatrixDouble pts)
private

Definition at line 542 of file TriPolynomialBase.cpp.

542  {
544 
545  EntitiesFieldData &data = cTx->dAta;
546  const FieldApproximationBase base = cTx->bAse;
547  if (cTx->basePolynomialsType0 == NULL)
548  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
549  "Polynomial type not set");
550  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
551  double *diffL, const int dim) =
553 
554  int nb_gauss_pts = pts.size2();
555 
556  double *PHI_f_e[3];
557  double *PHI_f;
558 
559  N_face_edge.resize(1, 3, false);
560  N_face_bubble.resize(1, false);
561  int face_order = data.dataOnEntities[MBTRI][0].getOrder();
562  // three edges on face
563  for (int ee = 0; ee < 3; ee++) {
564  N_face_edge(0, ee).resize(
565  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(face_order), false);
566  PHI_f_e[ee] = &((N_face_edge(0, ee))(0, 0));
567  }
568  N_face_bubble[0].resize(nb_gauss_pts,
569  3 * NBFACETRI_AINSWORTH_FACE_HDIV(face_order), false);
570  PHI_f = &*(N_face_bubble[0].data().begin());
571 
572  int face_nodes[3] = {0, 1, 2};
574  face_nodes, face_order,
575  &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0), NULL, PHI_f_e, NULL,
576  nb_gauss_pts, 3, base_polynomials);
578  face_nodes, face_order,
579  &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0), NULL, PHI_f, NULL,
580  nb_gauss_pts, 3, base_polynomials);
581 
582  // set shape functions into data structure
583  if (data.dataOnEntities[MBTRI].size() != 1) {
584  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
585  }
586  data.dataOnEntities[MBTRI][0].getN(base).resize(
587  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_HDIV(face_order), false);
588  int col = 0;
589  for (int oo = 0; oo < face_order; oo++) {
590  for (int ee = 0; ee < 3; ee++) {
591  for (int dd = 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
592  dd < 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++, col++) {
593  for (int gg = 0; gg < nb_gauss_pts; gg++) {
594  data.dataOnEntities[MBTRI][0].getN(base)(gg, col) =
595  N_face_edge(0, ee)(gg, dd);
596  }
597  }
598  }
599  for (int dd = 3 * NBFACETRI_AINSWORTH_FACE_HDIV(oo);
600  dd < 3 * NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++, col++) {
601  for (int gg = 0; gg < nb_gauss_pts; gg++) {
602  data.dataOnEntities[MBTRI][0].getN(base)(gg, col) =
603  N_face_bubble[0](gg, dd);
604  }
605  }
606  }
607 
609 }
#define NBFACETRI_AINSWORTH_FACE_HDIV(P)
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)
#define NBFACETRI_AINSWORTH_HDIV(P)
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
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
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
ublas::matrix< MatrixDouble > N_face_edge
ublas::vector< MatrixDouble > N_face_bubble

◆ getValueHdivDemkowiczBase()

MoFEMErrorCode TriPolynomialBase::getValueHdivDemkowiczBase ( MatrixDouble pts)
private

Definition at line 611 of file TriPolynomialBase.cpp.

611  {
613 
614  EntitiesFieldData &data = cTx->dAta;
615  // set shape functions into data structure
616  if (data.dataOnEntities[MBTRI].size() != 1) {
617  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
618  }
619  const FieldApproximationBase base = cTx->bAse;
620  if (base != DEMKOWICZ_JACOBI_BASE) {
621  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
622  "This should be used only with DEMKOWICZ_JACOBI_BASE "
623  "but base is %s",
624  ApproximationBaseNames[base]);
625  }
626  int order = data.dataOnEntities[MBTRI][0].getOrder();
627  int nb_gauss_pts = pts.size2();
628  data.dataOnEntities[MBTRI][0].getN(base).resize(
629  nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
630  double *phi_f = &*data.dataOnEntities[MBTRI][0].getN(base).data().begin();
633  int face_nodes[3] = {0, 1, 2};
635  face_nodes, order,
636  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
637  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(), phi_f,
638  NULL, nb_gauss_pts, 3);
639 
641 }
#define NBFACETRI_DEMKOWICZ_HDIV(P)
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

◆ getValueL2()

MoFEMErrorCode TriPolynomialBase::getValueL2 ( MatrixDouble pts)
private

Definition at line 338 of file TriPolynomialBase.cpp.

338  {
340 
341  switch (cTx->bAse) {
345  break;
348  break;
349  default:
350  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
351  }
352 
354 }
MoFEMErrorCode getValueL2BernsteinBezierBase(MatrixDouble &pts)
MoFEMErrorCode getValueL2AinsworthBase(MatrixDouble &pts)

◆ getValueL2AinsworthBase()

MoFEMErrorCode TriPolynomialBase::getValueL2AinsworthBase ( MatrixDouble pts)
private

Definition at line 356 of file TriPolynomialBase.cpp.

356  {
358 
359  EntitiesFieldData &data = cTx->dAta;
360  const FieldApproximationBase base = cTx->bAse;
361  if (cTx->basePolynomialsType0 == NULL)
362  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
363  "Polynomial type not set");
364  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
365  double *diffL, const int dim) =
367 
368  int nb_gauss_pts = pts.size2();
369 
370  data.dataOnEntities[MBTRI][0].getN(base).resize(
371  nb_gauss_pts, NBFACETRI_L2(data.dataOnEntities[MBTRI][0].getOrder()),
372  false);
373  data.dataOnEntities[MBTRI][0].getDiffN(base).resize(
374  nb_gauss_pts,
375  2 * NBFACETRI_L2(data.dataOnEntities[MBTRI][0].getOrder()), false);
376 
378  data.dataOnEntities[MBTRI][0].getOrder(),
379  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
380  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
381  &*data.dataOnEntities[MBTRI][0].getN(base).data().begin(),
382  &*data.dataOnEntities[MBTRI][0].getDiffN(base).data().begin(),
383  nb_gauss_pts, base_polynomials);
384 
386 }
#define NBFACETRI_L2(P)
Number of base functions on triangle for L2 space.
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

◆ getValueL2BernsteinBezierBase()

MoFEMErrorCode TriPolynomialBase::getValueL2BernsteinBezierBase ( MatrixDouble pts)
private

Definition at line 389 of file TriPolynomialBase.cpp.

389  {
391 
392  EntitiesFieldData &data = cTx->dAta;
393  const std::string &field_name = cTx->fieldName;
394  int nb_gauss_pts = pts.size2();
395 
396  auto get_alpha = [field_name](auto &data) -> MatrixInt & {
397  auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
398  if (!ptr)
399  ptr.reset(new MatrixInt());
400  return *ptr;
401  };
402 
403  auto get_base = [field_name](auto &data) -> MatrixDouble & {
404  auto &ptr = data.getBBNSharedPtr(field_name);
405  if (!ptr)
406  ptr.reset(new MatrixDouble());
407  return *ptr;
408  };
409 
410  auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
411  auto &ptr = data.getBBDiffNSharedPtr(field_name);
412  if (!ptr)
413  ptr.reset(new MatrixDouble());
414  return *ptr;
415  };
416 
417  auto get_alpha_by_name_ptr =
418  [](auto &data,
419  const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
420  return data.getBBAlphaIndicesSharedPtr(field_name);
421  };
422 
423  auto get_base_by_name_ptr =
424  [](auto &data,
425  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
426  return data.getBBNSharedPtr(field_name);
427  };
428 
429  auto get_diff_base_by_name_ptr =
430  [](auto &data,
431  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
432  return data.getBBDiffNSharedPtr(field_name);
433  };
434 
435  auto get_alpha_by_order_ptr =
436  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
437  return data.getBBAlphaIndicesByOrderSharedPtr(o);
438  };
439 
440  auto get_base_by_order_ptr =
441  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
442  return data.getBBNByOrderSharedPtr(o);
443  };
444 
445  auto get_diff_base_by_order_ptr =
446  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
447  return data.getBBDiffNByOrderSharedPtr(o);
448  };
449 
450  if (data.spacesOnEntities[MBTRI].test(L2)) {
451  if (data.dataOnEntities[MBTRI].size() != 1)
452  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
453  "Wrong size ent of ent data");
454 
455  if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
456  (unsigned int)nb_gauss_pts)
457  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
458  "Base functions or nodes has wrong number of integration points "
459  "for base %s",
461  auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
462 
463  auto &ent_data = data.dataOnEntities[MBTRI][0];
464  const int order = ent_data.getOrder();
465  const int nb_dofs = NBFACETRI_L2(order);
466 
467  if (get_alpha_by_order_ptr(ent_data, order)) {
468  get_alpha_by_name_ptr(ent_data, field_name) =
469  get_alpha_by_order_ptr(ent_data, order);
470  get_base_by_name_ptr(ent_data, field_name) =
471  get_base_by_order_ptr(ent_data, order);
472  get_diff_base_by_name_ptr(ent_data, field_name) =
473  get_diff_base_by_order_ptr(ent_data, order);
474  } else {
475 
476  auto &get_n = get_base(ent_data);
477  auto &get_diff_n = get_diff_base(ent_data);
478  get_n.resize(nb_gauss_pts, nb_dofs, false);
479  get_diff_n.resize(nb_gauss_pts, 2 * nb_dofs, false);
480  if (nb_dofs) {
481 
482  if (order == 0) {
483 
484  if (nb_dofs != 1)
485  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
486  "Inconsistent number of DOFs");
487 
488  auto &tri_alpha = get_alpha(ent_data);
489  tri_alpha.clear();
490  get_n(0, 0) = 1;
491  get_diff_n.clear();
492 
493  } else {
494 
495  if (nb_dofs != 3 + 3 * NBEDGE_H1(order) + NBFACETRI_H1(order))
496  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
497  "Inconsistent number of DOFs");
498 
499  auto &tri_alpha = get_alpha(ent_data);
500  tri_alpha.resize(nb_dofs, 3, false);
501 
503  &tri_alpha(0, 0));
504 
505  if (order > 1) {
506  std::array<int, 3> face_n{order, order, order};
507  std::array<int *, 3> face_edge_ptr{
508  &tri_alpha(3, 0), &tri_alpha(3 + NBEDGE_H1(order), 0),
509  &tri_alpha(3 + 2 * NBEDGE_H1(order), 0)};
511  face_n.data(), face_edge_ptr.data());
512  if (order > 2)
514  order, &tri_alpha(3 + 3 * NBEDGE_H1(order), 0));
515  }
517  order, lambda.size1(), tri_alpha.size1(), &tri_alpha(0, 0),
518  &lambda(0, 0), Tools::diffShapeFunMBTRI.data(), &get_n(0, 0),
519  &get_diff_n(0, 0));
520  }
521 
522  get_alpha_by_order_ptr(ent_data, order) =
523  get_alpha_by_name_ptr(ent_data, field_name);
524  get_base_by_order_ptr(ent_data, order) =
525  get_base_by_name_ptr(ent_data, field_name);
526  get_diff_base_by_order_ptr(ent_data, order) =
527  get_diff_base_by_name_ptr(ent_data, field_name);
528  }
529  }
530  } else {
531  auto &ent_data = data.dataOnEntities[MBTRI][0];
532  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
533  auto &get_n = get_base(ent_data);
534  auto &get_diff_n = get_diff_base(ent_data);
535  get_n.resize(nb_gauss_pts, 0, false);
536  get_diff_n.resize(nb_gauss_pts, 0, false);
537  }
538 
540 }
static MoFEMErrorCode generateIndicesVertexTri(const int N, int *alpha)

◆ query_interface()

MoFEMErrorCode TriPolynomialBase::query_interface ( boost::typeindex::type_index  type_index,
UnknownInterface **  iface 
) const
virtual

Reimplemented from MoFEM::BaseFunction.

Definition at line 26 of file TriPolynomialBase.cpp.

27  {
29  *iface = const_cast<TriPolynomialBase *>(this);
31 }
Calculate base functions on triangle.

Member Data Documentation

◆ cTx

EntPolynomialBaseCtx* MoFEM::TriPolynomialBase::cTx
private

Definition at line 41 of file TriPolynomialBase.hpp.

◆ diffN_face_bubble

ublas::vector<MatrixDouble> MoFEM::TriPolynomialBase::diffN_face_bubble
private

Definition at line 54 of file TriPolynomialBase.hpp.

◆ diffN_face_edge

ublas::matrix<MatrixDouble> MoFEM::TriPolynomialBase::diffN_face_edge
private

Definition at line 53 of file TriPolynomialBase.hpp.

◆ N_face_bubble

ublas::vector<MatrixDouble> MoFEM::TriPolynomialBase::N_face_bubble
private

Definition at line 52 of file TriPolynomialBase.hpp.

◆ N_face_edge

ublas::matrix<MatrixDouble> MoFEM::TriPolynomialBase::N_face_edge
private

Definition at line 51 of file TriPolynomialBase.hpp.


The documentation for this struct was generated from the following files: