v0.14.0
VolumeElementForcesAndSourcesCore.cpp
Go to the documentation of this file.
1 /** \file VolumeElementForcesAndSourcesCore.cpp
2 
3 \brief Implementation of volume element
4 
5 */
6 
7 namespace MoFEM {
8 
10  Interface &m_field)
11  : ForcesAndSourcesCore(m_field), vOlume(elementMeasure),
12  meshPositionsFieldName("MESH_NODE_POSITIONS"), coords(24), jAc(3, 3),
13  invJac(3, 3), opSetInvJacH1(invJac),
14  opContravariantPiolaTransform(elementMeasure, jAc),
15  opCovariantPiolaTransform(invJac), opSetInvJacHdivAndHcurl(invJac),
16  tJac(&jAc(0, 0), &jAc(0, 1), &jAc(0, 2), &jAc(1, 0), &jAc(1, 1),
17  &jAc(1, 2), &jAc(2, 0), &jAc(2, 1), &jAc(2, 2)),
18  tInvJac(&invJac(0, 0), &invJac(0, 1), &invJac(0, 2), &invJac(1, 0),
19  &invJac(1, 1), &invJac(1, 2), &invJac(2, 0), &invJac(2, 1),
20  &invJac(2, 2)) {}
21 
23  Interface &m_field, const EntityType type)
25 
28 
29  int order_data = getMaxDataOrder();
30  int order_row = getMaxRowOrder();
31  int order_col = getMaxColOrder();
32  const auto type = numeredEntFiniteElementPtr->getEntType();
33 
34  auto get_rule_by_type = [&]() {
35  switch (type) {
36  case MBHEX:
37  return getRule(order_row + 1, order_col + 1, order_data + 1);
38  default:
39  return getRule(order_row, order_col, order_data);
40  }
41  };
42 
43  const int rule = get_rule_by_type();
44 
45  auto calc_base_for_tet = [&]() {
47  const size_t nb_gauss_pts = gaussPts.size2();
48  auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
49  auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
50  base.resize(nb_gauss_pts, 4, false);
51  diff_base.resize(nb_gauss_pts, 12, false);
52  CHKERR Tools::shapeFunMBTET(&*base.data().begin(), &gaussPts(0, 0),
53  &gaussPts(1, 0), &gaussPts(2, 0), nb_gauss_pts);
54  double *diff_shape_ptr = &*diff_base.data().begin();
55  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
56  for (int nn = 0; nn != 4; ++nn) {
57  for (int dd = 0; dd != 3; ++dd, ++diff_shape_ptr) {
58  *diff_shape_ptr = Tools::diffShapeFunMBTET[3 * nn + dd];
59  }
60  }
61  }
63  };
64 
65  auto calc_base_for_hex = [&]() {
67  const size_t nb_gauss_pts = gaussPts.size2();
68  auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
69  auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
70  base.resize(nb_gauss_pts, 8, false);
71  diff_base.resize(nb_gauss_pts, 24, false);
72  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
73  const double ksi = gaussPts(0, gg);
74  const double zeta = gaussPts(1, gg);
75  const double eta = gaussPts(2, gg);
76  base(gg, 0) = N_MBHEX0(ksi, zeta, eta);
77  base(gg, 1) = N_MBHEX1(ksi, zeta, eta);
78  base(gg, 2) = N_MBHEX2(ksi, zeta, eta);
79  base(gg, 3) = N_MBHEX3(ksi, zeta, eta);
80  base(gg, 4) = N_MBHEX4(ksi, zeta, eta);
81  base(gg, 5) = N_MBHEX5(ksi, zeta, eta);
82  base(gg, 6) = N_MBHEX6(ksi, zeta, eta);
83  base(gg, 7) = N_MBHEX7(ksi, zeta, eta);
84  diff_base(gg, 0 * 3 + 0) = diffN_MBHEX0x(zeta, eta);
85  diff_base(gg, 1 * 3 + 0) = diffN_MBHEX1x(zeta, eta);
86  diff_base(gg, 2 * 3 + 0) = diffN_MBHEX2x(zeta, eta);
87  diff_base(gg, 3 * 3 + 0) = diffN_MBHEX3x(zeta, eta);
88  diff_base(gg, 4 * 3 + 0) = diffN_MBHEX4x(zeta, eta);
89  diff_base(gg, 5 * 3 + 0) = diffN_MBHEX5x(zeta, eta);
90  diff_base(gg, 6 * 3 + 0) = diffN_MBHEX6x(zeta, eta);
91  diff_base(gg, 7 * 3 + 0) = diffN_MBHEX7x(zeta, eta);
92  diff_base(gg, 0 * 3 + 1) = diffN_MBHEX0y(ksi, eta);
93  diff_base(gg, 1 * 3 + 1) = diffN_MBHEX1y(ksi, eta);
94  diff_base(gg, 2 * 3 + 1) = diffN_MBHEX2y(ksi, eta);
95  diff_base(gg, 3 * 3 + 1) = diffN_MBHEX3y(ksi, eta);
96  diff_base(gg, 4 * 3 + 1) = diffN_MBHEX4y(ksi, eta);
97  diff_base(gg, 5 * 3 + 1) = diffN_MBHEX5y(ksi, eta);
98  diff_base(gg, 6 * 3 + 1) = diffN_MBHEX6y(ksi, eta);
99  diff_base(gg, 7 * 3 + 1) = diffN_MBHEX7y(ksi, eta);
100  diff_base(gg, 0 * 3 + 2) = diffN_MBHEX0z(ksi, zeta);
101  diff_base(gg, 1 * 3 + 2) = diffN_MBHEX1z(ksi, zeta);
102  diff_base(gg, 2 * 3 + 2) = diffN_MBHEX2z(ksi, zeta);
103  diff_base(gg, 3 * 3 + 2) = diffN_MBHEX3z(ksi, zeta);
104  diff_base(gg, 4 * 3 + 2) = diffN_MBHEX4z(ksi, zeta);
105  diff_base(gg, 5 * 3 + 2) = diffN_MBHEX5z(ksi, zeta);
106  diff_base(gg, 6 * 3 + 2) = diffN_MBHEX6z(ksi, zeta);
107  diff_base(gg, 7 * 3 + 2) = diffN_MBHEX7z(ksi, zeta);
108  }
110  };
111 
112  auto set_integration_pts_for_hex = [&]() {
115  rule);
117  };
118 
119  auto set_integration_pts_for_tet = [&]() {
121  if (rule < QUAD_3D_TABLE_SIZE) {
122  if (QUAD_3D_TABLE[rule]->dim != 3) {
123  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
124  }
125  if (QUAD_3D_TABLE[rule]->order < rule) {
127  "wrong order %d != %d", QUAD_3D_TABLE[rule]->order, rule);
128  }
129  size_t nb_gauss_pts = QUAD_3D_TABLE[rule]->npoints;
130  gaussPts.resize(4, nb_gauss_pts, false);
131  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[1], 4,
132  &gaussPts(0, 0), 1);
133  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[2], 4,
134  &gaussPts(1, 0), 1);
135  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[3], 4,
136  &gaussPts(2, 0), 1);
137  cblas_dcopy(nb_gauss_pts, QUAD_3D_TABLE[rule]->weights, 1,
138  &gaussPts(3, 0), 1);
139 
140  // CHKERR calc_base_for_tet();
141 
142  auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
143  auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
144  base.resize(nb_gauss_pts, 4, false);
145  diff_base.resize(nb_gauss_pts, 12, false);
146  double *shape_ptr = &*base.data().begin();
147  cblas_dcopy(4 * nb_gauss_pts, QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
148  1);
149  double *diff_shape_ptr = &*diff_base.data().begin();
150  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
151  for (int nn = 0; nn != 4; ++nn) {
152  for (int dd = 0; dd != 3; ++dd, ++diff_shape_ptr) {
153  *diff_shape_ptr = Tools::diffShapeFunMBTET[3 * nn + dd];
154  }
155  }
156  }
157 
158  } else {
159  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
160  "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
161  }
163  };
164 
165  if (rule >= 0) {
166  switch (type) {
167  case MBTET:
168  CHKERR set_integration_pts_for_tet();
169  break;
170  case MBHEX:
171  CHKERR set_integration_pts_for_hex();
172  CHKERR calc_base_for_hex();
173  break;
174  default:
175  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
176  "Element type not implemented: %d", type);
177  }
178  } else {
179  // If rule is negative, set user defined integration points
180  CHKERR setGaussPts(order_row, order_col, order_data);
181  const size_t nb_gauss_pts = gaussPts.size2();
182  if (nb_gauss_pts) {
183  switch (type) {
184  case MBTET: {
185  CHKERR calc_base_for_tet();
186  } break;
187  case MBHEX:
188  CHKERR calc_base_for_hex();
189  break;
190  default:
191  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
192  "Element type not implemented: %d", type);
193  }
194  }
195  }
196 
198 }
199 
202  const auto ent = numeredEntFiniteElementPtr->getEnt();
203  const auto type = numeredEntFiniteElementPtr->getEntType();
204  CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
205  coords.resize(3 * num_nodes, false);
206  CHKERR mField.get_moab().get_coords(conn, num_nodes, &*coords.data().begin());
207 
208  auto get_tet_t_diff_n = [&]() {
212  };
213 
214  auto get_hex_t_diff_n = [&]() {
219  };
220 
221  auto get_t_diff_n = [&]() {
222  if (type == MBTET)
223  return get_tet_t_diff_n();
224  return get_hex_t_diff_n();
225  };
226 
227  auto t_diff_n = get_t_diff_n();
228 
230  &coords[0], &coords[1], &coords[2]);
233  jAc.clear();
234 
235  for (size_t n = 0; n != num_nodes; ++n) {
236  tJac(i, j) += t_coords(i) * t_diff_n(j);
237  ++t_coords;
238  ++t_diff_n;
239  }
242 
243  if (type == MBTET)
244  vOlume *= G_TET_W1[0] / 6.;
245 
246 #ifndef NDEBUG
247  if (vOlume <= std::numeric_limits<double>::epsilon() || vOlume != vOlume) {
248  MOFEM_LOG("SELF", Sev::warning) << "Volume is zero " << vOlume;
249  MOFEM_LOG("SELF", Sev::warning) << "Coords " << coords << endl;
250  }
251 #endif
252 
254 }
255 
259  // Get coords at Gauss points
261 
262  auto &shape_functions = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
263  double *shape_functions_ptr = &*shape_functions.data().begin();
264  const size_t nb_base_functions = shape_functions.size2();
265  const size_t nb_gauss_pts = gaussPts.size2();
266  coordsAtGaussPts.resize(nb_gauss_pts, 3, false);
267  coordsAtGaussPts.clear();
268  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_coords_at_gauss_ptr(
269  &coordsAtGaussPts(0, 0), &coordsAtGaussPts(0, 1),
270  &coordsAtGaussPts(0, 2));
272  shape_functions_ptr);
273  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
275  &coords[0], &coords[1], &coords[2]);
276  for (int bb = 0; bb != nb_base_functions; ++bb) {
277  t_coords_at_gauss_ptr(i) += t_coords(i) * t_shape_functions;
278  ++t_coords;
279  ++t_shape_functions;
280  };
281  ++t_coords_at_gauss_ptr;
282  }
284 }
285 
289 
292 
293  auto type = numeredEntFiniteElementPtr->getEntType();
294  auto dim_type = CN::Dimension(type);
295 
296  auto get_data_on_ents = [&](auto lower_dim, auto space) {
298  auto data = dataOnElement[space];
299  data->facesNodes = dataH1.facesNodes;
300  data->facesNodesOrder = dataH1.facesNodesOrder;
301  for (auto dd = dim_type; dd >= lower_dim; --dd) {
302  int nb_ents = moab::CN::NumSubEntities(type, dd);
303  for (int ii = 0; ii != nb_ents; ++ii) {
304  auto sub_ent_type = moab::CN::SubEntityType(type, dd, ii);
305  if ((dataH1.spacesOnEntities[sub_ent_type]).test(space)) {
306  auto &data_on_ent = data->dataOnEntities[sub_ent_type];
307  CHKERR getEntitySense(sub_ent_type, data_on_ent);
308  CHKERR getEntityDataOrder(sub_ent_type, space, data_on_ent);
309  data->spacesOnEntities[sub_ent_type].set(space);
310  }
311  }
312  }
314  };
315 
316  CHKERR get_data_on_ents(1, H1); // H1
317  CHKERR get_data_on_ents(1, HCURL); // Hcurl
318  CHKERR get_data_on_ents(2, HDIV); // Hdiv
319  CHKERR get_data_on_ents(3, L2); // L2
320 
322 }
323 
326 
327  if (numeredEntFiniteElementPtr->getEntType() == MBTET) {
328  // OK that code is not nice.
329  MatrixDouble new_diff_n;
330  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
332  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
333  EntitiesFieldData::EntData &data = dataH1.dataOnEntities[MBVERTEX][0];
334  if ((data.getDiffN(base).size1() == 4) &&
335  (data.getDiffN(base).size2() == 3)) {
336  const size_t nb_gauss_pts = gaussPts.size2();
337  const size_t nb_base_functions = 4;
338  new_diff_n.resize(nb_gauss_pts, 3 * nb_base_functions, false);
339  double *new_diff_n_ptr = &*new_diff_n.data().begin();
341  new_diff_n_ptr, &new_diff_n_ptr[1], &new_diff_n_ptr[2]);
342  double *t_diff_n_ptr = &*data.getDiffN(base).data().begin();
343  for (unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
345  t_diff_n_ptr, &t_diff_n_ptr[1], &t_diff_n_ptr[2]);
346  for (unsigned int bb = 0; bb != nb_base_functions; bb++) {
347  t_new_diff_n(i) = t_diff_n(i);
348  ++t_new_diff_n;
349  ++t_diff_n;
350  }
351  }
352  data.getDiffN(base).resize(new_diff_n.size1(), new_diff_n.size2(),
353  false);
354  data.getDiffN(base).swap(new_diff_n);
355  }
356  }
357  }
358 
359  if (dataH1.spacesOnEntities[MBVERTEX].test(H1))
361 
362  std::array<std::bitset<LASTSPACE>, 4> spaces_by_dim{
363 
364  std::bitset<LASTSPACE>(0), //
365  std::bitset<LASTSPACE>(0), //
366  std::bitset<LASTSPACE>(0), //
367  std::bitset<LASTSPACE>(0)
368 
369  };
370  for (auto type = MBEDGE; type != MBENTITYSET; ++type) {
371  spaces_by_dim[CN::Dimension(type)] |= dataH1.spacesOnEntities[type];
372  }
373 
374  if (
375 
376  spaces_by_dim[1].test(HCURL) || //
377  spaces_by_dim[2].test(HCURL) || //
378  spaces_by_dim[3].test(HCURL)
379 
380  ) {
383  }
384 
385  if (spaces_by_dim[2].test(HDIV) || spaces_by_dim[3].test(HDIV)) {
387  // Fix for tetrahedrons
388  if (numeredEntFiniteElementPtr->getEntType() == MBTET) {
389  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
390  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
391  for (auto t : {MBTRI, MBTET}) {
392  for (auto &d : dataHdiv.dataOnEntities[t]) {
393  d.getN(base) /= 6;
394  d.getDiffN(base) /= 6;
395  }
396  }
397  }
398  }
400  }
401 
402  if (spaces_by_dim[3].test(L2)) {
404  }
405 
407 }
408 
410  ForcesAndSourcesCore *ptr) {
412  if (!(ptrFE = dynamic_cast<VolumeElementForcesAndSourcesCore *>(ptr)))
413  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
414  "User operator and finite element do not work together");
416 }
417 
420 
421  const auto type = numeredEntFiniteElementPtr->getEntType();
423  switch (type) {
424  case MBTET:
426  boost::shared_ptr<BaseFunction>(new TetPolynomialBase(this));
427  break;
428  case MBHEX:
430  boost::shared_ptr<BaseFunction>(new HexPolynomialBase());
431  break;
432  default:
434  }
437  }
438 
442  if (gaussPts.size2() == 0)
447 
449 
450  // Iterate over operators
452 
454 }
455 
456 } // namespace MoFEM
UBlasMatrix< double >
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
MoFEM::VolumeElementForcesAndSourcesCore::setIntegrationPts
virtual MoFEMErrorCode setIntegrationPts()
Set integration points.
Definition: VolumeElementForcesAndSourcesCore.cpp:26
MoFEM::HexPolynomialBase
Calculate base functions on tetrahedral.
Definition: HexPolynomialBase.hpp:19
diffN_MBHEX0y
#define diffN_MBHEX0y(x, z)
Definition: fem_tools.h:87
MoFEM::EntitiesFieldData::facesNodesOrder
MatrixInt facesNodesOrder
order of face nodes on element
Definition: EntitiesFieldData.hpp:46
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::ForcesAndSourcesCore::calBernsteinBezierBaseFunctionsOnElement
MoFEMErrorCode calBernsteinBezierBaseFunctionsOnElement()
Calculate Bernstein-Bezier base.
Definition: ForcesAndSourcesCore.cpp:1132
MoFEM::DataOperator::opRhs
virtual MoFEMErrorCode opRhs(EntitiesFieldData &data, const bool error_if_no_base=false)
Definition: DataOperators.cpp:140
LASTBASE
@ LASTBASE
Definition: definitions.h:69
MoFEM::ForcesAndSourcesCore::getMaxRowOrder
int getMaxRowOrder() const
Get max order of approximation for field in rows.
Definition: ForcesAndSourcesCore.cpp:131
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
diffN_MBHEX1z
#define diffN_MBHEX1z(x, y)
Definition: fem_tools.h:96
N_MBHEX5
#define N_MBHEX5(x, y, z)
Definition: fem_tools.h:76
MoFEM::ForcesAndSourcesCore::getFaceNodes
MoFEMErrorCode getFaceNodes(EntitiesFieldData &data) const
Get nodes on faces.
Definition: ForcesAndSourcesCore.cpp:921
MoFEM::ForcesAndSourcesCore::lastEvaluatedElementEntityType
EntityType lastEvaluatedElementEntityType
Last evaluated type of element entity.
Definition: ForcesAndSourcesCore.hpp:489
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
NOBASE
@ NOBASE
Definition: definitions.h:59
diffN_MBHEX3x
#define diffN_MBHEX3x(y, z)
Definition: fem_tools.h:82
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
N_MBHEX0
#define N_MBHEX0(x, y, z)
Definition: fem_tools.h:71
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::VolumeElementForcesAndSourcesCore::opCovariantPiolaTransform
OpSetCovariantPiolaTransform opCovariantPiolaTransform
Definition: VolumeElementForcesAndSourcesCore.hpp:96
MoFEM::ForcesAndSourcesCore::setGaussPts
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
Definition: ForcesAndSourcesCore.cpp:1979
MoFEM::ForcesAndSourcesCore::getMaxDataOrder
int getMaxDataOrder() const
Get max order of approximation for data fields.
Definition: ForcesAndSourcesCore.cpp:120
MoFEM::Tools::diffShapeFunMBHEXAtCenter
static constexpr std::array< double, 24 > diffShapeFunMBHEXAtCenter
Definition: Tools.hpp:218
QUAD_::order
int order
Definition: quad.h:28
MoFEM::ForcesAndSourcesCore::dataHdiv
EntitiesFieldData & dataHdiv
Definition: ForcesAndSourcesCore.hpp:473
MoFEM::VolumeElementForcesAndSourcesCore::opSetInvJacHdivAndHcurl
OpSetInvJacHdivAndHcurl opSetInvJacHdivAndHcurl
Definition: VolumeElementForcesAndSourcesCore.hpp:97
MoFEM::Tools::diffShapeFunMBTET
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition: Tools.hpp:271
MoFEM::VolumeElementForcesAndSourcesCore::calculateCoordinatesAtGaussPts
virtual MoFEMErrorCode calculateCoordinatesAtGaussPts()
Calculate coordinate at integration points.
Definition: VolumeElementForcesAndSourcesCore.cpp:257
MoFEM::ForcesAndSourcesCore::getEntityDataOrder
MoFEMErrorCode getEntityDataOrder(const EntityType type, const FieldSpace space, boost::ptr_vector< EntitiesFieldData::EntData > &data) const
Get the entity data order.
Definition: ForcesAndSourcesCore.cpp:139
N_MBHEX7
#define N_MBHEX7(x, y, z)
Definition: fem_tools.h:78
MoFEM::EntitiesFieldData::EntData::getDiffN
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
Definition: EntitiesFieldData.hpp:1329
MoFEM::ForcesAndSourcesCore::getElementPolynomialBase
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
Definition: ForcesAndSourcesCore.hpp:90
MoFEM::ForcesAndSourcesCore::getMaxColOrder
int getMaxColOrder() const
Get max order of approximation for field in columns.
Definition: ForcesAndSourcesCore.cpp:135
MoFEM::Tools::outerProductOfEdgeIntegrationPtsForHex
static MoFEMErrorCode outerProductOfEdgeIntegrationPtsForHex(MatrixDouble &pts, const int edge0, const int edge1, const int edge2)
Definition: Tools.cpp:659
MoFEM::ForcesAndSourcesCore::calHierarchicalBaseFunctionsOnElement
MoFEMErrorCode calHierarchicalBaseFunctionsOnElement()
Calculate base functions.
Definition: ForcesAndSourcesCore.cpp:1112
diffN_MBHEX5x
#define diffN_MBHEX5x(y, z)
Definition: fem_tools.h:84
G_TET_W1
static const double G_TET_W1[]
Definition: fem_tools.h:1115
MoFEM::EntitiesFieldData::spacesOnEntities
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
Definition: EntitiesFieldData.hpp:49
diffN_MBHEX2y
#define diffN_MBHEX2y(x, z)
Definition: fem_tools.h:89
MoFEM::EntitiesFieldData::facesNodes
MatrixInt facesNodes
nodes on finite element faces
Definition: EntitiesFieldData.hpp:45
zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:130
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::VolumeElementForcesAndSourcesCore::coords
VectorDouble coords
Definition: VolumeElementForcesAndSourcesCore.hpp:90
MoFEM::invertTensor3by3
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:1588
diffN_MBHEX5y
#define diffN_MBHEX5y(x, z)
Definition: fem_tools.h:92
N_MBHEX3
#define N_MBHEX3(x, y, z)
Definition: fem_tools.h:74
MoFEM::ForcesAndSourcesCore::dataL2
EntitiesFieldData & dataL2
Definition: ForcesAndSourcesCore.hpp:474
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::VolumeElementForcesAndSourcesCore::num_nodes
int num_nodes
Definition: VolumeElementForcesAndSourcesCore.hpp:101
MoFEM::VolumeElementForcesAndSourcesCore::calculateVolumeAndJacobian
virtual MoFEMErrorCode calculateVolumeAndJacobian()
Calculate element volume and Jacobian.
Definition: VolumeElementForcesAndSourcesCore.cpp:200
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
diffN_MBHEX1x
#define diffN_MBHEX1x(y, z)
Definition: fem_tools.h:80
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::ForcesAndSourcesCore::coordsAtGaussPts
MatrixDouble coordsAtGaussPts
coordinated at gauss points
Definition: ForcesAndSourcesCore.hpp:544
diffN_MBHEX4z
#define diffN_MBHEX4z(x, y)
Definition: fem_tools.h:99
MoFEM::VolumeElementForcesAndSourcesCore::getSpaceBaseAndOrderOnElement
virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement()
Determine approximation space and order of base functions.
Definition: VolumeElementForcesAndSourcesCore.cpp:287
diffN_MBHEX2x
#define diffN_MBHEX2x(y, z)
Definition: fem_tools.h:81
MoFEM::VolumeElementForcesAndSourcesCore::VolumeElementForcesAndSourcesCore
VolumeElementForcesAndSourcesCore(Interface &m_field)
Definition: VolumeElementForcesAndSourcesCore.cpp:9
convert.type
type
Definition: convert.py:64
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator::setPtrFE
MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr)
Definition: VolumeElementForcesAndSourcesCore.cpp:409
MoFEM::ForcesAndSourcesCore::dataOnElement
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
Definition: ForcesAndSourcesCore.hpp:461
MoFEM::VolumeElementForcesAndSourcesCore::conn
const EntityHandle * conn
Definition: VolumeElementForcesAndSourcesCore.hpp:102
MoFEM::ForcesAndSourcesCore::createDataOnElement
MoFEMErrorCode createDataOnElement(EntityType type)
Create a entity data on element object.
Definition: ForcesAndSourcesCore.cpp:1337
diffN_MBHEX3z
#define diffN_MBHEX3z(x, y)
Definition: fem_tools.h:98
eta
double eta
Definition: free_surface.cpp:172
N_MBHEX1
#define N_MBHEX1(x, y, z)
Definition: fem_tools.h:72
diffN_MBHEX2z
#define diffN_MBHEX2z(x, y)
Definition: fem_tools.h:97
QUAD_::npoints
int npoints
Definition: quad.h:29
MoFEM::ForcesAndSourcesCore::getRule
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
Definition: ForcesAndSourcesCore.cpp:1973
diffN_MBHEX7x
#define diffN_MBHEX7x(y, z)
Definition: fem_tools.h:86
MoFEM::VolumeElementForcesAndSourcesCore::opSetInvJacH1
OpSetInvJacH1 opSetInvJacH1
Definition: VolumeElementForcesAndSourcesCore.hpp:94
MoFEM::VolumeElementForcesAndSourcesCore::tJac
FTensor::Tensor2< double *, 3, 3 > tJac
Definition: VolumeElementForcesAndSourcesCore.hpp:103
diffN_MBHEX7y
#define diffN_MBHEX7y(x, z)
Definition: fem_tools.h:94
diffN_MBHEX5z
#define diffN_MBHEX5z(x, y)
Definition: fem_tools.h:100
diffN_MBHEX3y
#define diffN_MBHEX3y(x, z)
Definition: fem_tools.h:90
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
diffN_MBHEX1y
#define diffN_MBHEX1y(x, z)
Definition: fem_tools.h:88
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
QUAD_3D_TABLE
static QUAD *const QUAD_3D_TABLE[]
Definition: quad.h:187
MoFEM::ForcesAndSourcesCore::getSpacesAndBaseOnEntities
MoFEMErrorCode getSpacesAndBaseOnEntities(EntitiesFieldData &data) const
Get field approximation space and base on entities.
Definition: ForcesAndSourcesCore.cpp:990
FTensor::Index< 'i', 3 >
convert.n
n
Definition: convert.py:82
diffN_MBHEX4y
#define diffN_MBHEX4y(x, z)
Definition: fem_tools.h:91
diffN_MBHEX6x
#define diffN_MBHEX6x(y, z)
Definition: fem_tools.h:85
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into EntitiesFieldData
Definition: ForcesAndSourcesCore.hpp:22
MoFEM::VolumeElementForcesAndSourcesCore::vOlume
double & vOlume
Definition: VolumeElementForcesAndSourcesCore.hpp:99
QUAD_3D_TABLE_SIZE
#define QUAD_3D_TABLE_SIZE
Definition: quad.h:186
MoFEM::ForcesAndSourcesCore::dataH1
EntitiesFieldData & dataH1
Definition: ForcesAndSourcesCore.hpp:471
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1540
MoFEM::Tools::shapeFunMBTET
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:747
diffN_MBHEX4x
#define diffN_MBHEX4x(y, z)
Definition: fem_tools.h:83
FTensor::dd
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
diffN_MBHEX6y
#define diffN_MBHEX6y(x, z)
Definition: fem_tools.h:93
FTensor::Tensor0
Definition: Tensor0.hpp:16
N_MBHEX2
#define N_MBHEX2(x, y, z)
Definition: fem_tools.h:73
MoFEM::VolumeElementForcesAndSourcesCore::jAc
MatrixDouble3by3 jAc
Definition: VolumeElementForcesAndSourcesCore.hpp:91
N_MBHEX4
#define N_MBHEX4(x, y, z)
Definition: fem_tools.h:75
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
MoFEM::ForcesAndSourcesCore::dataHcurl
EntitiesFieldData & dataHcurl
Definition: ForcesAndSourcesCore.hpp:472
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::ForcesAndSourcesCore::gaussPts
MatrixDouble gaussPts
Matrix of integration points.
Definition: ForcesAndSourcesCore.hpp:109
MoFEM::VolumeElementForcesAndSourcesCore::opContravariantPiolaTransform
OpSetContravariantPiolaTransform opContravariantPiolaTransform
Definition: VolumeElementForcesAndSourcesCore.hpp:95
diffN_MBHEX0z
#define diffN_MBHEX0z(x, y)
Definition: fem_tools.h:95
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
diffN_MBHEX6z
#define diffN_MBHEX6z(x, y)
Definition: fem_tools.h:101
MoFEM::EntitiesFieldData::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition: EntitiesFieldData.hpp:57
MoFEM::ForcesAndSourcesCore::loopOverOperators
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
Definition: ForcesAndSourcesCore.cpp:1379
MoFEM::FEMethod::numeredEntFiniteElementPtr
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
Definition: LoopMethods.hpp:387
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MoFEM::VolumeElementForcesAndSourcesCore::operator()
MoFEMErrorCode operator()()
function is run for every finite element
Definition: VolumeElementForcesAndSourcesCore.cpp:418
N_MBHEX6
#define N_MBHEX6(x, y, z)
Definition: fem_tools.h:77
invJac
MatrixDouble invJac
Definition: HookeElement.hpp:683
MoFEM::VolumeElementForcesAndSourcesCore::tInvJac
FTensor::Tensor2< double *, 3, 3 > tInvJac
Definition: VolumeElementForcesAndSourcesCore.hpp:104
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
diffN_MBHEX7z
#define diffN_MBHEX7z(x, y)
Definition: fem_tools.h:102
MoFEM::ForcesAndSourcesCore::UserDataOperator::ptrFE
ForcesAndSourcesCore * ptrFE
Definition: ForcesAndSourcesCore.hpp:986
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
diffN_MBHEX0x
#define diffN_MBHEX0x(y, z)
Definition: fem_tools.h:79
MoFEM::VolumeElementForcesAndSourcesCore::transformBaseFunctions
virtual MoFEMErrorCode transformBaseFunctions()
Transform base functions based on geometric element Jacobian.
Definition: VolumeElementForcesAndSourcesCore.cpp:324
MoFEM::ForcesAndSourcesCore::mField
Interface & mField
Definition: ForcesAndSourcesCore.hpp:24
MoFEM::TetPolynomialBase
Calculate base functions on tetrahedral.
Definition: TetPolynomialBase.hpp:17
MoFEM::ForcesAndSourcesCore::getEntitySense
MoFEMErrorCode getEntitySense(const EntityType type, boost::ptr_vector< EntitiesFieldData::EntData > &data) const
get sense (orientation) of entity
Definition: ForcesAndSourcesCore.cpp:84