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

Calculate base functions on tetrahedral. More...

#include <src/approximation/TetPolynomialBase.hpp>

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

Public Member Functions

MoFEMErrorCode query_interface (const MOFEMuuid &uuid, BaseFunctionUnknownInterface **iface) const
 
 TetPolynomialBase ()
 
 ~TetPolynomialBase ()
 
MoFEMErrorCode getValue (MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
 
- Public Member Functions inherited from MoFEM::BaseFunction
 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
 

Private Member Functions

MoFEMErrorCode getValueH1 (MatrixDouble &pts)
 
MoFEMErrorCode getValueH1AinsworthBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueH1BernsteinBezierBase (MatrixDouble &pts)
 
MoFEMErrorCode getValueL2 (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::vector< MatrixDoubleN_volume_edge
 
ublas::vector< MatrixDoubleN_volume_face
 
MatrixDouble N_volume_bubble
 
ublas::matrix< MatrixDoublediffN_face_edge
 
ublas::vector< MatrixDoublediffN_face_bubble
 
ublas::vector< MatrixDoublediffN_volume_edge
 
ublas::vector< MatrixDoublediffN_volume_face
 
MatrixDouble diffN_volume_bubble
 
MatrixInt senseFaceAlpha
 

Detailed Description

Calculate base functions on tetrahedral.

Definition at line 31 of file TetPolynomialBase.hpp.

Constructor & Destructor Documentation

◆ TetPolynomialBase()

TetPolynomialBase::TetPolynomialBase ( )

Definition at line 42 of file TetPolynomialBase.cpp.

42 {}

◆ ~TetPolynomialBase()

TetPolynomialBase::~TetPolynomialBase ( )

Definition at line 41 of file TetPolynomialBase.cpp.

41 {}

Member Function Documentation

◆ getValue()

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

Reimplemented from MoFEM::BaseFunction.

Definition at line 1178 of file TetPolynomialBase.cpp.

1179  {
1181 
1183  CHKERR ctx_ptr->query_interface(IDD_TET_BASE_FUNCTION, &iface);
1184  cTx = reinterpret_cast<EntPolynomialBaseCtx *>(iface);
1185 
1186  int nb_gauss_pts = pts.size2();
1187  if (!nb_gauss_pts)
1189 
1190  if (pts.size1() < 3)
1191  SETERRQ(
1192  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1193  "Wrong dimension of pts, should be at least 3 rows with coordinates");
1194 
1196  const FieldApproximationBase base = cTx->bAse;
1198  if (cTx->copyNodeBase == LASTBASE) {
1199  data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4,
1200  false);
1202  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1203  &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
1204  } else {
1205  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(base) =
1206  data.dataOnEntities[MBVERTEX][0].getNSharedPtr(cTx->copyNodeBase);
1207  }
1208  if (data.dataOnEntities[MBVERTEX][0].getN(base).size1() !=
1209  (unsigned int)nb_gauss_pts) {
1210  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1211  "Base functions or nodes has wrong number of integration points "
1212  "for base %s",
1213  ApproximationBaseNames[base]);
1214  }
1215  data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, false);
1216  std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
1217  data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
1218  }
1219 
1220  switch (cTx->sPace) {
1221  case H1:
1222  CHKERR getValueH1(pts);
1223  break;
1224  case HDIV:
1225  CHKERR getValueHdiv(pts);
1226  break;
1227  case HCURL:
1228  CHKERR getValueHcurl(pts);
1229  break;
1230  case L2:
1231  CHKERR getValueL2(pts);
1232  break;
1233  default:
1234  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Unknown space");
1235  }
1236 
1238 }
field with continuous normal traction
Definition: definitions.h:173
data structure for finite element entityIt keeps that about indices of degrees of freedom,...
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
const FieldApproximationBase bAse
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition: Tools.hpp:143
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
MoFEMErrorCode getValueL2(MatrixDouble &pts)
static const MOFEMuuid IDD_TET_BASE_FUNCTION
const FieldApproximationBase copyNodeBase
static const char *const ApproximationBaseNames[]
Definition: definitions.h:158
FieldApproximationBase
approximation base
Definition: definitions.h:144
DataForcesAndSourcesCore & dAta
MoFEMErrorCode getValueHdiv(MatrixDouble &pts)
field with continuous tangents
Definition: definitions.h:172
MoFEMErrorCode getValueH1(MatrixDouble &pts)
MoFEMErrorCode getValueHcurl(MatrixDouble &pts)
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
continuous field
Definition: definitions.h:171
EntPolynomialBaseCtx * cTx
static MoFEMErrorCode shapeFunMBTET(double *shape, const double *ksi, const double *eta, const double *zeta, const double nb)
Calculate shape functions on tetrahedron.
Definition: Tools.hpp:448
field with C-1 continuity
Definition: definitions.h:174

◆ getValueH1()

MoFEMErrorCode TetPolynomialBase::getValueH1 ( MatrixDouble pts)
private

Definition at line 44 of file TetPolynomialBase.cpp.

44  {
46 
47  switch (cTx->bAse) {
51  break;
54  break;
55  default:
56  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
57  }
58 
60 }
const FieldApproximationBase bAse
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
MoFEMErrorCode getValueH1AinsworthBase(MatrixDouble &pts)
MoFEMErrorCode getValueH1BernsteinBezierBase(MatrixDouble &pts)
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:146
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
EntPolynomialBaseCtx * cTx

◆ getValueH1AinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueH1AinsworthBase ( MatrixDouble pts)
private

Definition at line 62 of file TetPolynomialBase.cpp.

62  {
64 
66  const FieldApproximationBase base = cTx->bAse;
67  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
68  double *diffL, const int dim) =
70 
71  int nb_gauss_pts = pts.size2();
72 
73  int sense[6], order[6];
74  if (data.spacesOnEntities[MBEDGE].test(H1)) {
75  // edges
76  if (data.dataOnEntities[MBEDGE].size() != 6) {
77  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
78  }
79  double *h1_edge_n[6], *diff_h1_egde_n[6];
80  for (int ee = 0; ee != 6; ++ee) {
81  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
82  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
83  "data inconsistency");
84  }
85  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
86  order[ee] = data.dataOnEntities[MBEDGE][ee].getDataOrder();
87  int nb_dofs = NBEDGE_H1(data.dataOnEntities[MBEDGE][ee].getDataOrder());
88  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, nb_dofs,
89  false);
90  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
91  3 * nb_dofs, false);
92  h1_edge_n[ee] =
93  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
94  diff_h1_egde_n[ee] =
95  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
96  }
98  sense, order,
99  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
100  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
101  h1_edge_n, diff_h1_egde_n, nb_gauss_pts, base_polynomials);
102  } else {
103  for (int ee = 0; ee != 6; ++ee) {
104  data.dataOnEntities[MBEDGE][ee].getN(base).resize(0, 0, false);
105  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(0, 0, false);
106  }
107  }
108 
109  if (data.spacesOnEntities[MBTRI].test(H1)) {
110  // faces
111  if (data.dataOnEntities[MBTRI].size() != 4) {
112  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
113  }
114  double *h1_face_n[4], *diff_h1_face_n[4];
115  for (int ff = 0; ff != 4; ++ff) {
116  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
117  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
118  "data inconsistency");
119  }
120  int nb_dofs = NBFACETRI_H1(data.dataOnEntities[MBTRI][ff].getDataOrder());
121  order[ff] = data.dataOnEntities[MBTRI][ff].getDataOrder();
122  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, nb_dofs,
123  false);
124  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
125  3 * nb_dofs, false);
126  h1_face_n[ff] =
127  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
128  diff_h1_face_n[ff] =
129  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
130  }
131  if (data.facesNodes.size1() != 4) {
132  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
133  }
134  if (data.facesNodes.size2() != 3) {
135  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
136  }
138  &*data.facesNodes.data().begin(), order,
139  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
140  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
141  h1_face_n, diff_h1_face_n, nb_gauss_pts, base_polynomials);
142 
143  } else {
144  for (int ff = 0; ff != 4; ++ff) {
145  data.dataOnEntities[MBTRI][ff].getN(base).resize(0, false);
146  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(0, 0, false);
147  }
148  }
149 
150  if (data.spacesOnEntities[MBTET].test(H1)) {
151  // volume
152  int order = data.dataOnEntities[MBTET][0].getDataOrder();
153  int nb_vol_dofs = NBVOLUMETET_H1(order);
154  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, nb_vol_dofs,
155  false);
156  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
157  3 * nb_vol_dofs, false);
159  data.dataOnEntities[MBTET][0].getDataOrder(),
160  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
161  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
162  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
163  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
164  nb_gauss_pts, base_polynomials);
165  } else {
166  data.dataOnEntities[MBTET][0].getN(base).resize(0, 0, false);
167  data.dataOnEntities[MBTET][0].getDiffN(base).resize(0, 0, false);
168  }
169 
171 }
PetscErrorCode H1_FaceShapeFunctions_MBTET(int *faces_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:317
PetscErrorCode H1_VolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *volumeN, double *diff_volumeN, int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: h1.c:419
#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_MBTET(int *sense, int *p, double *N, double *diffN, double *edgeN[], double *diff_edgeN[], int GDIM, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: h1.c:245
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
const FieldApproximationBase bAse
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
MatrixInt facesNodes
nodes on finite element faces
FieldApproximationBase
approximation base
Definition: definitions.h:144
DataForcesAndSourcesCore & dAta
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
#define CHKERR
Inline error check.
Definition: definitions.h:596
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
continuous field
Definition: definitions.h:171
EntPolynomialBaseCtx * cTx

◆ getValueH1BernsteinBezierBase()

MoFEMErrorCode TetPolynomialBase::getValueH1BernsteinBezierBase ( MatrixDouble pts)
private

Definition at line 174 of file TetPolynomialBase.cpp.

174  {
176 
178  const std::string field_name = cTx->fieldName;
179  const int nb_gauss_pts = pts.size2();
180 
181  if (data.dataOnEntities[MBVERTEX][0].getN(NOBASE).size1() !=
182  (unsigned int)nb_gauss_pts)
183  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
184  "Base functions or nodes has wrong number of integration points "
185  "for base %s",
187  auto &lambda = data.dataOnEntities[MBVERTEX][0].getN(NOBASE);
188 
189  auto get_alpha = [field_name](auto &data) -> MatrixInt & {
190  auto &ptr = data.getBBAlphaIndicesSharedPtr(field_name);
191  if (!ptr)
192  ptr.reset(new MatrixInt());
193  return *ptr;
194  };
195 
196  auto get_base = [field_name](auto &data) -> MatrixDouble & {
197  auto &ptr = data.getBBNSharedPtr(field_name);
198  if (!ptr)
199  ptr.reset(new MatrixDouble());
200  return *ptr;
201  };
202 
203  auto get_diff_base = [field_name](auto &data) -> MatrixDouble & {
204  auto &ptr = data.getBBDiffNSharedPtr(field_name);
205  if (!ptr)
206  ptr.reset(new MatrixDouble());
207  return *ptr;
208  };
209 
210  auto get_alpha_by_name_ptr =
211  [](auto &data,
212  const std::string &field_name) -> boost::shared_ptr<MatrixInt> & {
213  return data.getBBAlphaIndicesSharedPtr(field_name);
214  };
215 
216  auto get_base_by_name_ptr =
217  [](auto &data,
218  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
219  return data.getBBNSharedPtr(field_name);
220  };
221 
222  auto get_diff_base_by_name_ptr =
223  [](auto &data,
224  const std::string &field_name) -> boost::shared_ptr<MatrixDouble> & {
225  return data.getBBDiffNSharedPtr(field_name);
226  };
227 
228  auto get_alpha_by_order_ptr =
229  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixInt> & {
230  return data.getBBAlphaIndicesByOrderSharedPtr(o);
231  };
232 
233  auto get_base_by_order_ptr =
234  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
235  return data.getBBNByOrderSharedPtr(o);
236  };
237 
238  auto get_diff_base_by_order_ptr =
239  [](auto &data, const size_t o) -> boost::shared_ptr<MatrixDouble> & {
240  return data.getBBDiffNByOrderSharedPtr(o);
241  };
242 
243  auto &vert_ent_data = data.dataOnEntities[MBVERTEX][0];
244  auto &vertex_alpha = get_alpha(vert_ent_data);
245  vertex_alpha.resize(4, 4, false);
246  vertex_alpha.clear();
247  for (int n = 0; n != 4; ++n)
248  vertex_alpha(n, n) = data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n];
249 
250  auto &vert_get_n = get_base(vert_ent_data);
251  auto &vert_get_diff_n = get_diff_base(vert_ent_data);
252  vert_get_n.resize(nb_gauss_pts, 4, false);
253  vert_get_diff_n.resize(nb_gauss_pts, 12, false);
255  1, lambda.size1(), vertex_alpha.size1(), &vertex_alpha(0, 0),
256  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &vert_get_n(0, 0),
257  &vert_get_diff_n(0, 0));
258  for (int n = 0; n != 4; ++n) {
259  const double f = boost::math::factorial<double>(
260  data.dataOnEntities[MBVERTEX][0].getBBNodeOrder()[n]);
261  for (int g = 0; g != nb_gauss_pts; ++g) {
262  vert_get_n(g, n) *= f;
263  for (int d = 0; d != 3; ++d)
264  vert_get_diff_n(g, 3 * n + d) *= f;
265  }
266  }
267 
268  // edges
269  if (data.spacesOnEntities[MBEDGE].test(H1)) {
270  if (data.dataOnEntities[MBEDGE].size() != 6)
271  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
272  "Wrong size of ent data");
273 
274  constexpr int edges_nodes[6][2] = {{0, 1}, {1, 2}, {2, 0},
275  {0, 3}, {1, 3}, {2, 3}};
276  for (int ee = 0; ee != 6; ++ee) {
277  auto &ent_data = data.dataOnEntities[MBEDGE][ee];
278  const int sense = ent_data.getSense();
279  if (sense == 0)
280  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
281  "Sense of the edge unknown");
282  const int order = ent_data.getDataOrder();
283  const int nb_dofs = NBEDGE_H1(order);
284 
285  if (nb_dofs) {
286  if (get_alpha_by_order_ptr(ent_data, order)) {
287  get_alpha_by_name_ptr(ent_data, field_name) =
288  get_alpha_by_order_ptr(ent_data, order);
289  get_base_by_name_ptr(ent_data, field_name) =
290  get_base_by_order_ptr(ent_data, order);
291  get_diff_base_by_name_ptr(ent_data, field_name) =
292  get_diff_base_by_order_ptr(ent_data, order);
293  } else {
294  auto &get_n = get_base(ent_data);
295  auto &get_diff_n = get_diff_base(ent_data);
296  get_n.resize(nb_gauss_pts, nb_dofs, false);
297  get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
298 
299  auto &edge_alpha = get_alpha(data.dataOnEntities[MBEDGE][ee]);
300  edge_alpha.resize(nb_dofs, 4, false);
302  &edge_alpha(0, 0));
303  if (sense == -1) {
304  for (int i = 0; i != edge_alpha.size1(); ++i) {
305  int a = edge_alpha(i, edges_nodes[ee][0]);
306  edge_alpha(i, edges_nodes[ee][0]) =
307  edge_alpha(i, edges_nodes[ee][1]);
308  edge_alpha(i, edges_nodes[ee][1]) = a;
309  }
310  }
312  order, lambda.size1(), edge_alpha.size1(), &edge_alpha(0, 0),
313  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
314  &get_diff_n(0, 0));
315 
316  get_alpha_by_order_ptr(ent_data, order) =
317  get_alpha_by_name_ptr(ent_data, field_name);
318  get_base_by_order_ptr(ent_data, order) =
319  get_base_by_name_ptr(ent_data, field_name);
320  get_diff_base_by_order_ptr(ent_data, order) =
321  get_diff_base_by_name_ptr(ent_data, field_name);
322  }
323  }
324  }
325  } else {
326  for (int ee = 0; ee != 6; ++ee) {
327  auto &ent_data = data.dataOnEntities[MBEDGE][ee];
328  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
329  auto &get_n = get_base(ent_data);
330  auto &get_diff_n = get_diff_base(ent_data);
331  get_n.resize(nb_gauss_pts, 0, false);
332  get_diff_n.resize(nb_gauss_pts, 0, false);
333  }
334  }
335 
336  // face
337  if (data.spacesOnEntities[MBTRI].test(H1)) {
338  if (data.dataOnEntities[MBTRI].size() != 4)
339  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
340  "Wrong size of ent data");
341  if (data.facesNodes.size1() != 4)
342  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
343  if (data.facesNodes.size2() != 3)
344  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
345 
346  for (int ff = 0; ff != 4; ++ff) {
347  auto &ent_data = data.dataOnEntities[MBTRI][ff];
348  const int order = ent_data.getDataOrder();
349  const int nb_dofs = NBFACETRI_H1(order);
350 
351  if (nb_dofs) {
352  if (get_alpha_by_order_ptr(ent_data, order)) {
353  get_alpha_by_name_ptr(ent_data, field_name) =
354  get_alpha_by_order_ptr(ent_data, order);
355  get_base_by_name_ptr(ent_data, field_name) =
356  get_base_by_order_ptr(ent_data, order);
357  get_diff_base_by_name_ptr(ent_data, field_name) =
358  get_diff_base_by_order_ptr(ent_data, order);
359  } else {
360 
361  auto &get_n = get_base(ent_data);
362  auto &get_diff_n = get_diff_base(ent_data);
363  get_n.resize(nb_gauss_pts, nb_dofs, false);
364  get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
365 
366  auto &face_alpha = get_alpha(ent_data );
367  face_alpha.resize(nb_dofs, 4, false);
368 
370  &face_alpha(0, 0));
371  senseFaceAlpha.resize(face_alpha.size1(), face_alpha.size2(), false);
372  senseFaceAlpha.clear();
373  constexpr int tri_nodes[4][3] = {
374  {0, 1, 3}, {1, 2, 3}, {0, 2, 3}, {0, 1, 2}};
375  for (int d = 0; d != nb_dofs; ++d)
376  for (int n = 0; n != 3; ++n)
377  senseFaceAlpha(d, data.facesNodes(ff, n)) =
378  face_alpha(d, tri_nodes[ff][n]);
379  face_alpha.swap(senseFaceAlpha);
381  order, lambda.size1(), face_alpha.size1(), &face_alpha(0, 0),
382  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
383  &get_diff_n(0, 0));
384 
385  get_alpha_by_order_ptr(ent_data, order) =
386  get_alpha_by_name_ptr(ent_data, field_name);
387  get_base_by_order_ptr(ent_data, order) =
388  get_base_by_name_ptr(ent_data, field_name);
389  get_diff_base_by_order_ptr(ent_data, order) =
390  get_diff_base_by_name_ptr(ent_data, field_name);
391  }
392  }
393  }
394  } else {
395  for (int ff = 0; ff != 4; ++ff) {
396  auto &ent_data = data.dataOnEntities[MBTRI][ff];
397  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
398  auto &get_n = get_base(ent_data);
399  auto &get_diff_n = get_diff_base(ent_data);
400  get_n.resize(nb_gauss_pts, 0, false);
401  get_diff_n.resize(nb_gauss_pts, 0, false);
402  }
403  }
404 
405  if (data.spacesOnEntities[MBTET].test(H1)) {
406  if (data.dataOnEntities[MBTET].size() != 1)
407  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
408  "Wrong size ent of ent data");
409 
410  auto &ent_data = data.dataOnEntities[MBTET][0];
411  const int order = ent_data.getDataOrder();
412  const int nb_dofs = NBVOLUMETET_H1(order);
413  if (get_alpha_by_order_ptr(ent_data, order)) {
414  get_alpha_by_name_ptr(ent_data, field_name) =
415  get_alpha_by_order_ptr(ent_data, order);
416  get_base_by_name_ptr(ent_data, field_name) =
417  get_base_by_order_ptr(ent_data, order);
418  get_diff_base_by_name_ptr(ent_data, field_name) =
419  get_diff_base_by_order_ptr(ent_data, order);
420  } else {
421 
422  auto &get_n = get_base(ent_data);
423  auto &get_diff_n = get_diff_base(ent_data);
424  get_n.resize(nb_gauss_pts, nb_dofs, false);
425  get_diff_n.resize(nb_gauss_pts, 3 * nb_dofs, false);
426  if (nb_dofs) {
427  auto &tet_alpha = get_alpha(ent_data);
428  tet_alpha.resize(nb_dofs, 4, false);
429 
430  CHKERR BernsteinBezier::generateIndicesTetTet(order, &tet_alpha(0, 0));
432  order, lambda.size1(), tet_alpha.size1(), &tet_alpha(0, 0),
433  &lambda(0, 0), Tools::diffShapeFunMBTET.data(), &get_n(0, 0),
434  &get_diff_n(0, 0));
435 
436  get_alpha_by_order_ptr(ent_data, order) =
437  get_alpha_by_name_ptr(ent_data, field_name);
438  get_base_by_order_ptr(ent_data, order) =
439  get_base_by_name_ptr(ent_data, field_name);
440  get_diff_base_by_order_ptr(ent_data, order) =
441  get_diff_base_by_name_ptr(ent_data, field_name);
442  }
443  }
444  } else {
445  auto &ent_data = data.dataOnEntities[MBTET][0];
446  ent_data.getBBAlphaIndicesSharedPtr(field_name).reset();
447  auto &get_n = get_base(ent_data);
448  auto &get_diff_n = get_diff_base(ent_data);
449  get_n.resize(nb_gauss_pts, 0, false);
450  get_diff_n.resize(nb_gauss_pts, 0, false);
451  }
452 
454 }
#define NBEDGE_H1(P)
Numer of base function on edge for H1 space.
static MoFEMErrorCode generateIndicesTetTet(const int N, int *alpha)
data structure for finite element entityIt keeps that about indices of degrees of freedom,...
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
#define NBVOLUMETET_H1(P)
Number of base functions on tetrahedron for H1 space.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition: Tools.hpp:143
static MoFEMErrorCode generateIndicesTriTet(const int N[], int *alpha[])
static MoFEMErrorCode baseFunctionsTet(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< int, ublas::row_major, IntAllocator > MatrixInt
Definition: Types.hpp:73
static const char *const ApproximationBaseNames[]
Definition: definitions.h:158
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
MatrixInt facesNodes
nodes on finite element faces
static MoFEMErrorCode generateIndicesEdgeTet(const int N[], int *alpha[])
DataForcesAndSourcesCore & dAta
#define NBFACETRI_H1(P)
Number of base function on triangle for H1 space.
#define CHKERR
Inline error check.
Definition: definitions.h:596
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
continuous field
Definition: definitions.h:171
EntPolynomialBaseCtx * cTx

◆ getValueHcurl()

MoFEMErrorCode TetPolynomialBase::getValueHcurl ( MatrixDouble pts)
private

Definition at line 1159 of file TetPolynomialBase.cpp.

1159  {
1161 
1162  switch (cTx->bAse) {
1166  break;
1167  case DEMKOWICZ_JACOBI_BASE:
1169  break;
1170  default:
1171  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
1172  }
1173 
1175 }
MoFEMErrorCode getValueHcurlAinsworthBase(MatrixDouble &pts)
const FieldApproximationBase bAse
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:146
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
EntPolynomialBaseCtx * cTx
MoFEMErrorCode getValueHcurlDemkowiczBase(MatrixDouble &pts)

◆ getValueHcurlAinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueHcurlAinsworthBase ( MatrixDouble pts)
private

Definition at line 909 of file TetPolynomialBase.cpp.

909  {
911 
913  const FieldApproximationBase base = cTx->bAse;
914  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
915  double *diffL, const int dim) =
917 
918  int nb_gauss_pts = pts.size2();
919 
920  // edges
921  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
922  int sense[6], order[6];
923  if (data.dataOnEntities[MBEDGE].size() != 6) {
924  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
925  }
926  double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
927  for (int ee = 0; ee != 6; ee++) {
928  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
929  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
930  "data inconsistency");
931  }
932  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
933  order[ee] = data.dataOnEntities[MBEDGE][ee].getDataOrder();
934  int nb_dofs = NBEDGE_AINSWORTH_HCURL(
935  data.dataOnEntities[MBEDGE][ee].getDataOrder());
936  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
937  3 * nb_dofs, false);
938  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
939  9 * nb_dofs, false);
940  hcurl_edge_n[ee] =
941  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
942  diff_hcurl_edge_n[ee] =
943  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
944  }
946  sense, order,
947  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
948  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
949  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts, base_polynomials);
950  } else {
951  for (int ee = 0; ee != 6; ee++) {
952  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
953  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
954  false);
955  }
956  }
957 
958  // triangles
959  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
960  int order[4];
961  // faces
962  if (data.dataOnEntities[MBTRI].size() != 4) {
963  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
964  }
965  double *hcurl_base_n[4], *diff_hcurl_base_n[4];
966  for (int ff = 0; ff != 4; ff++) {
967  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
968  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
969  "data inconsistency");
970  }
971  order[ff] = data.dataOnEntities[MBTRI][ff].getDataOrder();
972  int nb_dofs = NBFACETRI_AINSWORTH_HCURL(order[ff]);
973  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
974  3 * nb_dofs, false);
975  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
976  9 * nb_dofs, false);
977  hcurl_base_n[ff] =
978  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
979  diff_hcurl_base_n[ff] =
980  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
981  }
982  if (data.facesNodes.size1() != 4) {
983  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
984  }
985  if (data.facesNodes.size2() != 3) {
986  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
987  }
989  &*data.facesNodes.data().begin(), order,
990  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
991  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
992  hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts, base_polynomials);
993  } else {
994  for (int ff = 0; ff != 4; ff++) {
995  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
996  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
997  false);
998  }
999  }
1000 
1001  if (data.spacesOnEntities[MBTET].test(HCURL)) {
1002 
1003  // volume
1004  int order = data.dataOnEntities[MBTET][0].getDataOrder();
1005  int nb_vol_dofs = NBVOLUMETET_AINSWORTH_HCURL(order);
1006  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
1007  3 * nb_vol_dofs, false);
1008  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
1009  9 * nb_vol_dofs, false);
1011  data.dataOnEntities[MBTET][0].getDataOrder(),
1012  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1013  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1014  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
1015  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
1016  nb_gauss_pts, base_polynomials);
1017 
1018  } else {
1019  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
1020  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
1021  }
1022 
1024 }
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,...
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
const FieldApproximationBase bAse
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
MatrixInt facesNodes
nodes on finite element faces
FieldApproximationBase
approximation base
Definition: definitions.h:144
DataForcesAndSourcesCore & dAta
#define NBEDGE_AINSWORTH_HCURL(P)
MoFEMErrorCode Hcurl_Ainsworth_FaceFunctions_MBTET(int *face_nodes, int *p, double *N, double *diffN, double *phi_f[4], double *diff_phi_f[4], 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:1052
field with continuous tangents
Definition: definitions.h:172
#define CHKERR
Inline error check.
Definition: definitions.h:596
MoFEMErrorCode Hcurl_Ainsworth_VolumeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v, double *diff_phi_v, int nb_integration_pts, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
H-curl volume base functions.
Definition: Hcurl.cpp:1403
MoFEMErrorCode Hcurl_Ainsworth_EdgeBaseFunctions_MBTET(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 tetrahedral.
Definition: Hcurl.cpp:16
#define NBVOLUMETET_AINSWORTH_HCURL(P)
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
#define NBFACETRI_AINSWORTH_HCURL(P)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
EntPolynomialBaseCtx * cTx

◆ getValueHcurlDemkowiczBase()

MoFEMErrorCode TetPolynomialBase::getValueHcurlDemkowiczBase ( MatrixDouble pts)
private

Definition at line 1027 of file TetPolynomialBase.cpp.

1027  {
1029 
1031  const FieldApproximationBase base = cTx->bAse;
1032  if (base != DEMKOWICZ_JACOBI_BASE) {
1033  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1034  "This should be used only with DEMKOWICZ_JACOBI_BASE "
1035  "but base is %s",
1036  ApproximationBaseNames[base]);
1037  }
1038 
1039  int nb_gauss_pts = pts.size2();
1040 
1041  // edges
1042  if (data.spacesOnEntities[MBEDGE].test(HCURL)) {
1043  int sense[6], order[6];
1044  if (data.dataOnEntities[MBEDGE].size() != 6) {
1045  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1046  "wrong size of data structure, expected space for six edges "
1047  "but is %d",
1048  data.dataOnEntities[MBEDGE].size());
1049  }
1050  double *hcurl_edge_n[6], *diff_hcurl_edge_n[6];
1051  for (int ee = 0; ee != 6; ee++) {
1052  if (data.dataOnEntities[MBEDGE][ee].getSense() == 0) {
1053  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1054  "orintation of edges is not set");
1055  }
1056  sense[ee] = data.dataOnEntities[MBEDGE][ee].getSense();
1057  order[ee] = data.dataOnEntities[MBEDGE][ee].getDataOrder();
1058  int nb_dofs = NBEDGE_DEMKOWICZ_HCURL(
1059  data.dataOnEntities[MBEDGE][ee].getDataOrder());
1060  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts,
1061  3 * nb_dofs, false);
1062  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts,
1063  9 * nb_dofs, false);
1064  hcurl_edge_n[ee] =
1065  &*data.dataOnEntities[MBEDGE][ee].getN(base).data().begin();
1066  diff_hcurl_edge_n[ee] =
1067  &*data.dataOnEntities[MBEDGE][ee].getDiffN(base).data().begin();
1068  }
1070  sense, order,
1071  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1072  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1073  hcurl_edge_n, diff_hcurl_edge_n, nb_gauss_pts);
1074  } else {
1075  // No DOFs on edges, resize base function matrices, indicating that no
1076  // dofs on them.
1077  for (int ee = 0; ee != 6; ee++) {
1078  data.dataOnEntities[MBEDGE][ee].getN(base).resize(nb_gauss_pts, 0, false);
1079  data.dataOnEntities[MBEDGE][ee].getDiffN(base).resize(nb_gauss_pts, 0,
1080  false);
1081  }
1082  }
1083 
1084  // triangles
1085  if (data.spacesOnEntities[MBTRI].test(HCURL)) {
1086  int order[4];
1087  // faces
1088  if (data.dataOnEntities[MBTRI].size() != 4) {
1089  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1090  "data structure for storing face h-curl base have wrong size "
1091  "should be four but is %d",
1092  data.dataOnEntities[MBTRI].size());
1093  }
1094  double *hcurl_base_n[4], *diff_hcurl_base_n[4];
1095  for (int ff = 0; ff != 4; ff++) {
1096  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
1097  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1098  "orintation of face is not set");
1099  }
1100  order[ff] = data.dataOnEntities[MBTRI][ff].getDataOrder();
1101  int nb_dofs = NBFACETRI_DEMKOWICZ_HCURL(order[ff]);
1102  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts,
1103  3 * nb_dofs, false);
1104  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts,
1105  9 * nb_dofs, false);
1106  hcurl_base_n[ff] =
1107  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
1108  diff_hcurl_base_n[ff] =
1109  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
1110  }
1111  if (data.facesNodes.size1() != 4) {
1112  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1113  "data inconsistency, should be four faces");
1114  }
1115  if (data.facesNodes.size2() != 3) {
1116  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1117  "data inconsistency, should be three nodes on face");
1118  }
1120  &*data.facesNodes.data().begin(), order,
1121  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1122  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1123  hcurl_base_n, diff_hcurl_base_n, nb_gauss_pts);
1124  } else {
1125  // No DOFs on faces, resize base function matrices, indicating that no
1126  // dofs on them.
1127  for (int ff = 0; ff != 4; ff++) {
1128  data.dataOnEntities[MBTRI][ff].getN(base).resize(nb_gauss_pts, 0, false);
1129  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(nb_gauss_pts, 0,
1130  false);
1131  }
1132  }
1133 
1134  if (data.spacesOnEntities[MBTET].test(HCURL)) {
1135  // volume
1136  int order = data.dataOnEntities[MBTET][0].getDataOrder();
1137  int nb_vol_dofs = NBVOLUMETET_DEMKOWICZ_HCURL(order);
1138  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts,
1139  3 * nb_vol_dofs, false);
1140  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts,
1141  9 * nb_vol_dofs, false);
1143  data.dataOnEntities[MBTET][0].getDataOrder(),
1144  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1145  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
1146  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
1147  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
1148  nb_gauss_pts);
1149  } else {
1150  // No DOFs on faces, resize base function matrices, indicating that no
1151  // dofs on them.
1152  data.dataOnEntities[MBTET][0].getN(base).resize(nb_gauss_pts, 0, false);
1153  data.dataOnEntities[MBTET][0].getDiffN(base).resize(nb_gauss_pts, 0, false);
1154  }
1155 
1157 }
data structure for finite element entityIt keeps that about indices of degrees of freedom,...
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
MoFEMErrorCode Hcurl_Demkowicz_VolumeBaseFunctions_MBTET(int p, double *n, double *diff_n, double *phi, double *diff_phi, int nb_integration_pts)
Volume base interior function.
Definition: Hcurl.cpp:2465
const FieldApproximationBase bAse
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
static const char *const ApproximationBaseNames[]
Definition: definitions.h:158
MoFEMErrorCode Hcurl_Demkowicz_FaceBaseFunctions_MBTET(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:2392
MatrixInt facesNodes
nodes on finite element faces
FieldApproximationBase
approximation base
Definition: definitions.h:144
MoFEMErrorCode Hcurl_Demkowicz_EdgeBaseFunctions_MBTET(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 tetrahedral.
Definition: Hcurl.cpp:2068
DataForcesAndSourcesCore & dAta
field with continuous tangents
Definition: definitions.h:172
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define NBFACETRI_DEMKOWICZ_HCURL(P)
#define NBEDGE_DEMKOWICZ_HCURL(P)
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
EntPolynomialBaseCtx * cTx
#define NBVOLUMETET_DEMKOWICZ_HCURL(P)

◆ getValueHdiv()

MoFEMErrorCode TetPolynomialBase::getValueHdiv ( MatrixDouble pts)
private

Definition at line 892 of file TetPolynomialBase.cpp.

892  {
894 
895  switch (cTx->bAse) {
898  return getValueHdivAinsworthBase(pts);
900  return getValueHdivDemkowiczBase(pts);
901  default:
902  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not implemented");
903  }
904 
906 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
MoFEMErrorCode getValueHdivDemkowiczBase(MatrixDouble &pts)
const FieldApproximationBase bAse
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
MoFEMErrorCode getValueHdivAinsworthBase(MatrixDouble &pts)
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:146
EntPolynomialBaseCtx * cTx

◆ getValueHdivAinsworthBase()

MoFEMErrorCode TetPolynomialBase::getValueHdivAinsworthBase ( MatrixDouble pts)
private

Definition at line 485 of file TetPolynomialBase.cpp.

485  {
487 
489  const FieldApproximationBase base = cTx->bAse;
490  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
491  double *diffL, const int dim) =
493 
494  int nb_gauss_pts = pts.size2();
495 
496  // face shape functions
497 
498  double *phi_f_e[4][3];
499  double *phi_f[4];
500  double *diff_phi_f_e[4][3];
501  double *diff_phi_f[4];
502 
503  N_face_edge.resize(4, 3, false);
504  N_face_bubble.resize(4, false);
505  diffN_face_edge.resize(4, 3, false);
506  diffN_face_bubble.resize(4, false);
507 
508  int faces_order[4];
509  for (int ff = 0; ff != 4; ++ff) {
510  if (data.dataOnEntities[MBTRI][ff].getSense() == 0) {
511  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
512  }
513  faces_order[ff] = data.dataOnEntities[MBTRI][ff].getDataOrder();
514  // three edges on face
515  for (int ee = 0; ee < 3; ee++) {
516  N_face_edge(ff, ee).resize(
517  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_EDGE_HDIV(faces_order[ff]),
518  false);
519  diffN_face_edge(ff, ee).resize(
520  nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_EDGE_HDIV(faces_order[ff]),
521  false);
522  phi_f_e[ff][ee] = &((N_face_edge(ff, ee))(0, 0));
523  diff_phi_f_e[ff][ee] = &((diffN_face_edge(ff, ee))(0, 0));
524  }
525  N_face_bubble[ff].resize(nb_gauss_pts,
526  3 * NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]),
527  false);
528  diffN_face_bubble[ff].resize(
529  nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]),
530  false);
531  phi_f[ff] = &*(N_face_bubble[ff].data().begin());
532  diff_phi_f[ff] = &*(diffN_face_bubble[ff].data().begin());
533  }
534 
536  &data.facesNodes(0, 0), faces_order,
537  &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
538  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_f_e,
539  diff_phi_f_e, nb_gauss_pts, base_polynomials);
540 
542  &data.facesNodes(0, 0), faces_order,
543  &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
544  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_f, diff_phi_f,
545  nb_gauss_pts, base_polynomials);
546 
547  // volume shape functions
548 
549  double *phi_v_e[6];
550  double *phi_v_f[4];
551  double *phi_v;
552  double *diff_phi_v_e[6];
553  double *diff_phi_v_f[4];
554  double *diff_phi_v;
555 
556  int volume_order = data.dataOnEntities[MBTET][0].getDataOrder();
557 
558  N_volume_edge.resize(6, false);
559  diffN_volume_edge.resize(6, false);
560  for (int ee = 0; ee != 6; ++ee) {
561  N_volume_edge[ee].resize(
562  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_EDGE_HDIV(volume_order), false);
563  diffN_volume_edge[ee].resize(
564  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_EDGE_HDIV(volume_order), false);
565  phi_v_e[ee] = &*(N_volume_edge[ee].data().begin());
566  diff_phi_v_e[ee] = &*(diffN_volume_edge[ee].data().begin());
567  }
569  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
570  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_v_e,
571  diff_phi_v_e, nb_gauss_pts, base_polynomials);
572 
573  N_volume_face.resize(4, false);
574  diffN_volume_face.resize(4, false);
575  for (int ff = 0; ff != 4; ++ff) {
576  N_volume_face[ff].resize(
577  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order), false);
578  diffN_volume_face[ff].resize(
579  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order), false);
580  phi_v_f[ff] = &*(N_volume_face[ff].data().begin());
581  diff_phi_v_f[ff] = &*(diffN_volume_face[ff].data().begin());
582  }
584  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
585  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_v_f,
586  diff_phi_v_f, nb_gauss_pts, base_polynomials);
587 
588  N_volume_bubble.resize(
589  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order), false);
590  diffN_volume_bubble.resize(
591  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order), false);
592  phi_v = &*(N_volume_bubble.data().begin());
593  diff_phi_v = &*(diffN_volume_bubble.data().begin());
595  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
596  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), phi_v, diff_phi_v,
597  nb_gauss_pts, base_polynomials);
598 
599  // Set shape functions into data structure Shape functions hast to be put
600  // in arrays in order which guarantee hierarchical series of degrees of
601  // freedom, i.e. in other words dofs form sub-entities has to be group
602  // by order.
603 
606 
607  // faces
608  if (data.dataOnEntities[MBTRI].size() != 4) {
609  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "data inconsistency");
610  }
611  for (int ff = 0; ff != 4; ff++) {
612  data.dataOnEntities[MBTRI][ff].getN(base).resize(
613  nb_gauss_pts, 3 * NBFACETRI_AINSWORTH_HDIV(faces_order[ff]), false);
614  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
615  nb_gauss_pts, 9 * NBFACETRI_AINSWORTH_HDIV(faces_order[ff]), false);
616  if (NBFACETRI_AINSWORTH_HDIV(faces_order[ff]) == 0)
617  continue;
618  // face
619  double *base_ptr =
620  &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
621  FTensor::Tensor1<double *, 3> t_base(base_ptr, &base_ptr[HVEC1],
622  &base_ptr[HVEC2], 3);
623  double *diff_base_ptr =
624  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
626  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
627  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
628  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
629  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
630  &diff_base_ptr[HVEC2_2], 9);
631  // face-face
632  boost::shared_ptr<FTensor::Tensor1<double *, 3>> t_base_f;
633  boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>> t_diff_base_f;
634  if (NBFACETRI_AINSWORTH_FACE_HDIV(faces_order[ff]) > 0) {
635  base_ptr = phi_f[ff];
636  t_base_f = boost::shared_ptr<FTensor::Tensor1<double *, 3>>(
637  new FTensor::Tensor1<double *, 3>(base_ptr, &base_ptr[HVEC1],
638  &base_ptr[HVEC2], 3));
639  diff_base_ptr = diff_phi_f[ff];
640  t_diff_base_f = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>>(
642  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
643  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
644  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
645  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
646  &diff_base_ptr[HVEC2_2], 9));
647  }
648  // edge-face
649  base_ptr = phi_f_e[ff][0];
650  FTensor::Tensor1<double *, 3> t_base_f_e0(base_ptr, &base_ptr[HVEC1],
651  &base_ptr[HVEC2], 3);
652  diff_base_ptr = diff_phi_f_e[ff][0];
653  FTensor::Tensor2<double *, 3, 3> t_diff_base_f_e0(
654  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
655  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
656  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
657  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
658  &diff_base_ptr[HVEC2_2], 9);
659  base_ptr = phi_f_e[ff][1];
660  FTensor::Tensor1<double *, 3> t_base_f_e1(base_ptr, &base_ptr[HVEC1],
661  &base_ptr[HVEC2], 3);
662  diff_base_ptr = diff_phi_f_e[ff][1];
663  FTensor::Tensor2<double *, 3, 3> t_diff_base_f_e1(
664  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
665  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
666  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
667  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
668  &diff_base_ptr[HVEC2_2], 9);
669  base_ptr = phi_f_e[ff][2];
670  FTensor::Tensor1<double *, 3> t_base_f_e2(base_ptr, &base_ptr[HVEC1],
671  &base_ptr[HVEC2], 3);
672  diff_base_ptr = diff_phi_f_e[ff][2];
673  FTensor::Tensor2<double *, 3, 3> t_diff_base_f_e2(
674  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
675  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
676  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
677  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
678  &diff_base_ptr[HVEC2_2], 9);
679  for (int gg = 0; gg != nb_gauss_pts; gg++) {
680  for (int oo = 0; oo != faces_order[ff]; oo++) {
681  for (int dd = NBFACETRI_AINSWORTH_EDGE_HDIV(oo);
682  dd != NBFACETRI_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
683  t_base(i) = t_base_f_e0(i);
684  ++t_base;
685  ++t_base_f_e0;
686  t_diff_base(i, j) = t_diff_base_f_e0(i, j);
687  ++t_diff_base;
688  ++t_diff_base_f_e0;
689  t_base(i) = t_base_f_e1(i);
690  ++t_base;
691  ++t_base_f_e1;
692  t_diff_base(i, j) = t_diff_base_f_e1(i, j);
693  ++t_diff_base;
694  ++t_diff_base_f_e1;
695  t_base(i) = t_base_f_e2(i);
696  ++t_base;
697  ++t_base_f_e2;
698  t_diff_base(i, j) = t_diff_base_f_e2(i, j);
699  ++t_diff_base;
700  ++t_diff_base_f_e2;
701  }
702  for (int dd = NBFACETRI_AINSWORTH_FACE_HDIV(oo);
703  dd != NBFACETRI_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
704  t_base(i) = (*t_base_f)(i);
705  ++t_base;
706  ++(*t_base_f);
707  t_diff_base(i, j) = (*t_diff_base_f)(i, j);
708  ++t_diff_base;
709  ++(*t_diff_base_f);
710  }
711  }
712  }
713  }
714 
715  // volume
716  data.dataOnEntities[MBTET][0].getN(base).resize(
717  nb_gauss_pts, 3 * NBVOLUMETET_AINSWORTH_HDIV(volume_order), false);
718  data.dataOnEntities[MBTET][0].getDiffN(base).resize(
719  nb_gauss_pts, 9 * NBVOLUMETET_AINSWORTH_HDIV(volume_order), false);
720  if (NBVOLUMETET_AINSWORTH_HDIV(volume_order) > 0) {
721  double *base_ptr =
722  &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
723  FTensor::Tensor1<double *, 3> t_base(base_ptr, &base_ptr[HVEC1],
724  &base_ptr[HVEC2], 3);
725  double *diff_base_ptr =
726  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
728  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
729  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
730  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
731  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
732  &diff_base_ptr[HVEC2_2], 9);
733  // edges
734  std::vector<FTensor::Tensor1<double *, 3>> t_base_v_e;
735  t_base_v_e.reserve(6);
736  std::vector<FTensor::Tensor2<double *, 3, 3>> t_diff_base_v_e;
737  t_diff_base_v_e.reserve(6);
738  for (int ee = 0; ee != 6; ee++) {
739  base_ptr = phi_v_e[ee];
740  diff_base_ptr = diff_phi_v_e[ee];
741  t_base_v_e.push_back(FTensor::Tensor1<double *, 3>(
742  base_ptr, &base_ptr[HVEC1], &base_ptr[HVEC2], 3));
743  t_diff_base_v_e.push_back(FTensor::Tensor2<double *, 3, 3>(
744  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
745  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
746  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
747  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
748  &diff_base_ptr[HVEC2_2], 9));
749  }
750  // faces
751  std::vector<FTensor::Tensor1<double *, 3>> t_base_v_f;
752  t_base_v_f.reserve(4);
753  std::vector<FTensor::Tensor2<double *, 3, 3>> t_diff_base_v_f;
754  t_diff_base_v_f.reserve(4);
755  if (NBVOLUMETET_AINSWORTH_FACE_HDIV(volume_order) > 0) {
756  for (int ff = 0; ff != 4; ff++) {
757  base_ptr = phi_v_f[ff];
758  diff_base_ptr = diff_phi_v_f[ff];
759  t_base_v_f.push_back(FTensor::Tensor1<double *, 3>(
760  base_ptr, &base_ptr[HVEC1], &base_ptr[HVEC2], 3));
761  t_diff_base_v_f.push_back(FTensor::Tensor2<double *, 3, 3>(
762  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
763  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
764  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
765  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
766  &diff_base_ptr[HVEC2_2], 9));
767  }
768  }
769  boost::shared_ptr<FTensor::Tensor1<double *, 3>> t_base_v;
770  boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>> t_diff_base_v;
771  if (NBVOLUMETET_AINSWORTH_VOLUME_HDIV(volume_order) > 0) {
772  base_ptr = phi_v;
773  t_base_v = boost::shared_ptr<FTensor::Tensor1<double *, 3>>(
774  new FTensor::Tensor1<double *, 3>(base_ptr, &base_ptr[HVEC1],
775  &base_ptr[HVEC2], 3));
776  diff_base_ptr = diff_phi_v;
777  t_diff_base_v = boost::shared_ptr<FTensor::Tensor2<double *, 3, 3>>(
779  &diff_base_ptr[HVEC0_0], &diff_base_ptr[HVEC0_1],
780  &diff_base_ptr[HVEC0_2], &diff_base_ptr[HVEC1_0],
781  &diff_base_ptr[HVEC1_1], &diff_base_ptr[HVEC1_2],
782  &diff_base_ptr[HVEC2_0], &diff_base_ptr[HVEC2_1],
783  &diff_base_ptr[HVEC2_2], 9));
784  }
785  for (int gg = 0; gg != nb_gauss_pts; gg++) {
786  for (int oo = 0; oo < volume_order; oo++) {
787  for (int dd = NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo);
788  dd < NBVOLUMETET_AINSWORTH_EDGE_HDIV(oo + 1); dd++) {
789  for (int ee = 0; ee < 6; ee++) {
790  t_base(i) = t_base_v_e[ee](i);
791  ++t_base;
792  ++t_base_v_e[ee];
793  t_diff_base(i, j) = t_diff_base_v_e[ee](i, j);
794  ++t_diff_base;
795  ++t_diff_base_v_e[ee];
796  }
797  }
798  for (int dd = NBVOLUMETET_AINSWORTH_FACE_HDIV(oo);
799  dd < NBVOLUMETET_AINSWORTH_FACE_HDIV(oo + 1); dd++) {
800  for (int ff = 0; ff < 4; ff++) {
801  t_base(i) = t_base_v_f[ff](i);
802  ++t_base;
803  ++t_base_v_f[ff];
804  t_diff_base(i, j) = t_diff_base_v_f[ff](i, j);
805  ++t_diff_base;
806  ++t_diff_base_v_f[ff];
807  }
808  }
809  for (int dd = NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo);
810  dd < NBVOLUMETET_AINSWORTH_VOLUME_HDIV(oo + 1); dd++) {
811  t_base(i) = (*t_base_v)(i);
812  ++t_base;
813  ++(*t_base_v);
814  t_diff_base(i, j) = (*t_diff_base_v)(i, j);
815  ++t_diff_base;
816  ++(*t_diff_base_v);
817  }
818  }
819  }
820  }
821 
823 }
ublas::matrix< MatrixDouble > diffN_face_edge
#define NBFACETRI_AINSWORTH_EDGE_HDIV(P)
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,...
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
const FieldApproximationBase bAse
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define NBVOLUMETET_AINSWORTH_FACE_HDIV(P)
#define NBVOLUMETET_AINSWORTH_VOLUME_HDIV(P)
ublas::vector< MatrixDouble > N_volume_face
MoFEMErrorCode Hdiv_Ainsworth_VolumeBubbleShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v, double *diff_phi_v, int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Interior bubble functions by Ainsworth .
Definition: Hdiv.cpp:484
ublas::matrix< MatrixDouble > N_face_edge
MoFEMErrorCode Hdiv_Ainsworth_FaceBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_f[], double *diff_phi_v_f[], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Definition: Hdiv.cpp:368
MatrixInt facesNodes
nodes on finite element faces
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
FieldApproximationBase
approximation base
Definition: definitions.h:144
DataForcesAndSourcesCore & dAta
ublas::vector< MatrixDouble > diffN_volume_edge
ublas::vector< MatrixDouble > diffN_face_bubble
ublas::vector< MatrixDouble > N_volume_edge
#define CHKERR
Inline error check.
Definition: definitions.h:596
ublas::vector< MatrixDouble > diffN_volume_face
MoFEMErrorCode Hdiv_Ainsworth_EdgeBasedVolumeShapeFunctions_MBTET(int p, double *N, double *diffN, double *phi_v_e[6], double *diff_phi_v_e[6], int gdim, PetscErrorCode(*base_polynomials)(int p, double s, double *diff_s, double *L, double *diffL, const int dim))
Hdiv base function, Edge-based interior (volume) functions by Ainsworth .
Definition: Hdiv.cpp:280
#define NBVOLUMETET_AINSWORTH_HDIV(P)
MoFEMErrorCode Hdiv_Ainsworth_EdgeFaceShapeFunctions_MBTET(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f_e[4][3], double *diff_phi_f_e[4][3], int gdim, 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:13
#define NBFACETRI_AINSWORTH_FACE_HDIV(P)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
#define NBFACETRI_AINSWORTH_HDIV(P)
MoFEMErrorCode Hdiv_Ainsworth_FaceBubbleShapeFunctions(int *faces_nodes, int *p, double *N, double *diffN, double *phi_f[], double *diff_phi_f[], int gdim, 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:137
EntPolynomialBaseCtx * cTx
ublas::vector< MatrixDouble > N_face_bubble
#define NBVOLUMETET_AINSWORTH_EDGE_HDIV(P)

◆ getValueHdivDemkowiczBase()

MoFEMErrorCode TetPolynomialBase::getValueHdivDemkowiczBase ( MatrixDouble pts)
private

Definition at line 825 of file TetPolynomialBase.cpp.

825  {
827 
829  const FieldApproximationBase base = cTx->bAse;
830  if (base != DEMKOWICZ_JACOBI_BASE) {
831  SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
832  "This should be used only with DEMKOWICZ_JACOBI_BASE "
833  "but base is %s",
834  ApproximationBaseNames[base]);
835  }
836  int nb_gauss_pts = pts.size2();
837 
838  int volume_order = data.dataOnEntities[MBTET][0].getDataOrder();
839 
840  int p_f[4];
841  double *phi_f[4];
842  double *diff_phi_f[4];
843 
844  // Calculate base function on tet faces
845  for (int ff = 0; ff != 4; ff++) {
846  int face_order = data.dataOnEntities[MBTRI][ff].getDataOrder();
847  int order = volume_order > face_order ? volume_order : face_order;
848  data.dataOnEntities[MBTRI][ff].getN(base).resize(
849  nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
850  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
851  nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(order), false);
852  p_f[ff] = order;
853  phi_f[ff] = &*data.dataOnEntities[MBTRI][ff].getN(base).data().begin();
854  diff_phi_f[ff] =
855  &*data.dataOnEntities[MBTRI][ff].getDiffN(base).data().begin();
856  if (NBFACETRI_DEMKOWICZ_HDIV(order) == 0)
857  continue;
859  &data.facesNodes(ff, 0), order,
860  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
861  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
862  phi_f[ff], diff_phi_f[ff], nb_gauss_pts, 4);
863  }
864 
865  // Calculate base functions in tet interior
866  if (NBVOLUMETET_DEMKOWICZ_HDIV(volume_order) > 0) {
867  data.dataOnEntities[MBTET][0].getN(base).resize(
868  nb_gauss_pts, 3 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
869  data.dataOnEntities[MBTET][0].getDiffN(base).resize(
870  nb_gauss_pts, 9 * NBVOLUMETET_DEMKOWICZ_HDIV(volume_order), false);
871  double *phi_v = &*data.dataOnEntities[MBTET][0].getN(base).data().begin();
872  double *diff_phi_v =
873  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin();
875  volume_order, &data.dataOnEntities[MBVERTEX][0].getN(base)(0, 0),
876  &data.dataOnEntities[MBVERTEX][0].getDiffN(base)(0, 0), p_f, phi_f,
877  diff_phi_f, phi_v, diff_phi_v, nb_gauss_pts);
878  }
879 
880  // Set size of face base correctly
881  for (int ff = 0; ff != 4; ff++) {
882  int face_order = data.dataOnEntities[MBTRI][ff].getDataOrder();
883  data.dataOnEntities[MBTRI][ff].getN(base).resize(
884  nb_gauss_pts, 3 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
885  data.dataOnEntities[MBTRI][ff].getDiffN(base).resize(
886  nb_gauss_pts, 9 * NBFACETRI_DEMKOWICZ_HDIV(face_order), true);
887  }
888 
890 }
data structure for finite element entityIt keeps that about indices of degrees of freedom,...
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
const FieldApproximationBase bAse
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
static const char *const ApproximationBaseNames[]
Definition: definitions.h:158
MoFEMErrorCode Hdiv_Demkowicz_Interior_MBTET(int p, double *N, double *diffN, int p_face[], double *phi_f[4], double *diff_phi_f[4], double *phi_v, double *diff_phi_v, int gdim)
Definition: Hdiv.cpp:762
MatrixInt facesNodes
nodes on finite element faces
#define NBFACETRI_DEMKOWICZ_HDIV(P)
FieldApproximationBase
approximation base
Definition: definitions.h:144
DataForcesAndSourcesCore & dAta
#define NBVOLUMETET_DEMKOWICZ_HDIV(P)
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
EntPolynomialBaseCtx * cTx

◆ getValueL2()

MoFEMErrorCode TetPolynomialBase::getValueL2 ( MatrixDouble pts)
private

Definition at line 456 of file TetPolynomialBase.cpp.

456  {
458 
460  const FieldApproximationBase base = cTx->bAse;
461  PetscErrorCode (*base_polynomials)(int p, double s, double *diff_s, double *L,
462  double *diffL, const int dim) =
464 
465  int nb_gauss_pts = pts.size2();
466 
467  data.dataOnEntities[MBTET][0].getN(base).resize(
468  nb_gauss_pts,
469  NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getDataOrder()), false);
470  data.dataOnEntities[MBTET][0].getDiffN(base).resize(
471  nb_gauss_pts,
472  3 * NBVOLUMETET_L2(data.dataOnEntities[MBTET][0].getDataOrder()), false);
473 
475  data.dataOnEntities[MBTET][0].getDataOrder(),
476  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
477  &*data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin(),
478  &*data.dataOnEntities[MBTET][0].getN(base).data().begin(),
479  &*data.dataOnEntities[MBTET][0].getDiffN(base).data().begin(),
480  nb_gauss_pts, base_polynomials);
481 
483 }
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 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,...
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
const FieldApproximationBase bAse
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
FieldApproximationBase
approximation base
Definition: definitions.h:144
DataForcesAndSourcesCore & dAta
#define CHKERR
Inline error check.
Definition: definitions.h:596
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
EntPolynomialBaseCtx * cTx
PetscErrorCode L2_Ainsworth_ShapeFunctions_MBTET(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 tetrahedron for L2 space.
Definition: l2.c:84

◆ query_interface()

MoFEMErrorCode TetPolynomialBase::query_interface ( const MOFEMuuid uuid,
BaseFunctionUnknownInterface **  iface 
) const
virtual

Reimplemented from MoFEM::BaseFunction.

Definition at line 25 of file TetPolynomialBase.cpp.

26  {
27 
29  *iface = NULL;
30  if (uuid == IDD_TET_BASE_FUNCTION) {
31  *iface = const_cast<TetPolynomialBase *>(this);
33  } else {
34  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "wrong interference");
35  }
36  ierr = BaseFunction::query_interface(uuid, iface);
37  CHKERRG(ierr);
39 }
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, MoFEM::BaseFunctionUnknownInterface **iface) const
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:544
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
static const MOFEMuuid IDD_TET_BASE_FUNCTION

Member Data Documentation

◆ cTx

EntPolynomialBaseCtx* MoFEM::TetPolynomialBase::cTx
private

Definition at line 43 of file TetPolynomialBase.hpp.

◆ diffN_face_bubble

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

Definition at line 58 of file TetPolynomialBase.hpp.

◆ diffN_face_edge

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

Definition at line 57 of file TetPolynomialBase.hpp.

◆ diffN_volume_bubble

MatrixDouble MoFEM::TetPolynomialBase::diffN_volume_bubble
private

Definition at line 61 of file TetPolynomialBase.hpp.

◆ diffN_volume_edge

ublas::vector<MatrixDouble> MoFEM::TetPolynomialBase::diffN_volume_edge
private

Definition at line 59 of file TetPolynomialBase.hpp.

◆ diffN_volume_face

ublas::vector<MatrixDouble> MoFEM::TetPolynomialBase::diffN_volume_face
private

Definition at line 60 of file TetPolynomialBase.hpp.

◆ N_face_bubble

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

Definition at line 52 of file TetPolynomialBase.hpp.

◆ N_face_edge

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

Definition at line 51 of file TetPolynomialBase.hpp.

◆ N_volume_bubble

MatrixDouble MoFEM::TetPolynomialBase::N_volume_bubble
private

Definition at line 55 of file TetPolynomialBase.hpp.

◆ N_volume_edge

ublas::vector<MatrixDouble> MoFEM::TetPolynomialBase::N_volume_edge
private

Definition at line 53 of file TetPolynomialBase.hpp.

◆ N_volume_face

ublas::vector<MatrixDouble> MoFEM::TetPolynomialBase::N_volume_face
private

Definition at line 54 of file TetPolynomialBase.hpp.

◆ senseFaceAlpha

MatrixInt MoFEM::TetPolynomialBase::senseFaceAlpha
private

Definition at line 74 of file TetPolynomialBase.hpp.


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