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, const EntityType type)
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 
24 
25  int order_data = getMaxDataOrder();
26  int order_row = getMaxRowOrder();
27  int order_col = getMaxColOrder();
28  const auto type = numeredEntFiniteElementPtr->getEntType();
29 
30  auto get_rule_by_type = [&]() {
31  switch (type) {
32  case MBHEX:
33  return getRule(order_row + 1, order_col + 1, order_data + 1);
34  default:
35  return getRule(order_row, order_col, order_data);
36  }
37  };
38 
39  const int rule = get_rule_by_type();
40 
41  auto calc_base_for_tet = [&]() {
43  const size_t nb_gauss_pts = gaussPts.size2();
44  auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
45  auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
46  base.resize(nb_gauss_pts, 4, false);
47  diff_base.resize(nb_gauss_pts, 12, false);
48  CHKERR Tools::shapeFunMBTET(&*base.data().begin(), &gaussPts(0, 0),
49  &gaussPts(1, 0), &gaussPts(2, 0), nb_gauss_pts);
50  double *diff_shape_ptr = &*diff_base.data().begin();
51  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
52  for (int nn = 0; nn != 4; ++nn) {
53  for (int dd = 0; dd != 3; ++dd, ++diff_shape_ptr) {
54  *diff_shape_ptr = Tools::diffShapeFunMBTET[3 * nn + dd];
55  }
56  }
57  }
59  };
60 
61  auto calc_base_for_hex = [&]() {
63  const size_t nb_gauss_pts = gaussPts.size2();
64  auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
65  auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
66  base.resize(nb_gauss_pts, 8, false);
67  diff_base.resize(nb_gauss_pts, 24, false);
68  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
69  const double ksi = gaussPts(0, gg);
70  const double zeta = gaussPts(1, gg);
71  const double eta = gaussPts(2, gg);
72  base(gg, 0) = N_MBHEX0(ksi, zeta, eta);
73  base(gg, 1) = N_MBHEX1(ksi, zeta, eta);
74  base(gg, 2) = N_MBHEX2(ksi, zeta, eta);
75  base(gg, 3) = N_MBHEX3(ksi, zeta, eta);
76  base(gg, 4) = N_MBHEX4(ksi, zeta, eta);
77  base(gg, 5) = N_MBHEX5(ksi, zeta, eta);
78  base(gg, 6) = N_MBHEX6(ksi, zeta, eta);
79  base(gg, 7) = N_MBHEX7(ksi, zeta, eta);
80  diff_base(gg, 0 * 3 + 0) = diffN_MBHEX0x(zeta, eta);
81  diff_base(gg, 1 * 3 + 0) = diffN_MBHEX1x(zeta, eta);
82  diff_base(gg, 2 * 3 + 0) = diffN_MBHEX2x(zeta, eta);
83  diff_base(gg, 3 * 3 + 0) = diffN_MBHEX3x(zeta, eta);
84  diff_base(gg, 4 * 3 + 0) = diffN_MBHEX4x(zeta, eta);
85  diff_base(gg, 5 * 3 + 0) = diffN_MBHEX5x(zeta, eta);
86  diff_base(gg, 6 * 3 + 0) = diffN_MBHEX6x(zeta, eta);
87  diff_base(gg, 7 * 3 + 0) = diffN_MBHEX7x(zeta, eta);
88  diff_base(gg, 0 * 3 + 1) = diffN_MBHEX0y(ksi, eta);
89  diff_base(gg, 1 * 3 + 1) = diffN_MBHEX1y(ksi, eta);
90  diff_base(gg, 2 * 3 + 1) = diffN_MBHEX2y(ksi, eta);
91  diff_base(gg, 3 * 3 + 1) = diffN_MBHEX3y(ksi, eta);
92  diff_base(gg, 4 * 3 + 1) = diffN_MBHEX4y(ksi, eta);
93  diff_base(gg, 5 * 3 + 1) = diffN_MBHEX5y(ksi, eta);
94  diff_base(gg, 6 * 3 + 1) = diffN_MBHEX6y(ksi, eta);
95  diff_base(gg, 7 * 3 + 1) = diffN_MBHEX7y(ksi, eta);
96  diff_base(gg, 0 * 3 + 2) = diffN_MBHEX0z(ksi, zeta);
97  diff_base(gg, 1 * 3 + 2) = diffN_MBHEX1z(ksi, zeta);
98  diff_base(gg, 2 * 3 + 2) = diffN_MBHEX2z(ksi, zeta);
99  diff_base(gg, 3 * 3 + 2) = diffN_MBHEX3z(ksi, zeta);
100  diff_base(gg, 4 * 3 + 2) = diffN_MBHEX4z(ksi, zeta);
101  diff_base(gg, 5 * 3 + 2) = diffN_MBHEX5z(ksi, zeta);
102  diff_base(gg, 6 * 3 + 2) = diffN_MBHEX6z(ksi, zeta);
103  diff_base(gg, 7 * 3 + 2) = diffN_MBHEX7z(ksi, zeta);
104  }
106  };
107 
108  auto set_integration_pts_for_hex = [&]() {
111  rule);
113  };
114 
115  auto set_integration_pts_for_tet = [&]() {
117  if (rule < QUAD_3D_TABLE_SIZE) {
118  if (QUAD_3D_TABLE[rule]->dim != 3) {
119  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
120  }
121  if (QUAD_3D_TABLE[rule]->order < rule) {
123  "wrong order %d != %d", QUAD_3D_TABLE[rule]->order, rule);
124  }
125  size_t nb_gauss_pts = QUAD_3D_TABLE[rule]->npoints;
126  gaussPts.resize(4, nb_gauss_pts, false);
127  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[1], 4,
128  &gaussPts(0, 0), 1);
129  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[2], 4,
130  &gaussPts(1, 0), 1);
131  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[3], 4,
132  &gaussPts(2, 0), 1);
133  cblas_dcopy(nb_gauss_pts, QUAD_3D_TABLE[rule]->weights, 1,
134  &gaussPts(3, 0), 1);
135 
136  // CHKERR calc_base_for_tet();
137 
138  auto &base = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
139  auto &diff_base = dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE);
140  base.resize(nb_gauss_pts, 4, false);
141  diff_base.resize(nb_gauss_pts, 12, false);
142  double *shape_ptr = &*base.data().begin();
143  cblas_dcopy(4 * nb_gauss_pts, QUAD_3D_TABLE[rule]->points, 1, shape_ptr,
144  1);
145  double *diff_shape_ptr = &*diff_base.data().begin();
146  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
147  for (int nn = 0; nn != 4; ++nn) {
148  for (int dd = 0; dd != 3; ++dd, ++diff_shape_ptr) {
149  *diff_shape_ptr = Tools::diffShapeFunMBTET[3 * nn + dd];
150  }
151  }
152  }
153 
154  } else {
155  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
156  "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
157  }
159  };
160 
161  if (rule >= 0) {
162  switch (type) {
163  case MBTET:
164  CHKERR set_integration_pts_for_tet();
165  break;
166  case MBHEX:
167  CHKERR set_integration_pts_for_hex();
168  CHKERR calc_base_for_hex();
169  break;
170  default:
171  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
172  "Element type not implemented: %d", type);
173  }
174  } else {
175  // If rule is negative, set user defined integration points
176  CHKERR setGaussPts(order_row, order_col, order_data);
177  const size_t nb_gauss_pts = gaussPts.size2();
178  if (nb_gauss_pts) {
179  switch (type) {
180  case MBTET: {
181  CHKERR calc_base_for_tet();
182  } break;
183  case MBHEX:
184  CHKERR calc_base_for_hex();
185  break;
186  default:
187  SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
188  "Element type not implemented: %d", type);
189  }
190  }
191  }
192 
194 }
195 
198  const auto ent = numeredEntFiniteElementPtr->getEnt();
199  const auto type = numeredEntFiniteElementPtr->getEntType();
200  CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
201  coords.resize(3 * num_nodes, false);
202  CHKERR mField.get_moab().get_coords(conn, num_nodes, &*coords.data().begin());
203 
204  auto get_tet_t_diff_n = [&]() {
208  };
209 
210  auto get_hex_t_diff_n = [&]() {
215  };
216 
217  auto get_t_diff_n = [&]() {
218  if (type == MBTET)
219  return get_tet_t_diff_n();
220  return get_hex_t_diff_n();
221  };
222 
223  auto t_diff_n = get_t_diff_n();
224 
226  &coords[0], &coords[1], &coords[2]);
229  jAc.clear();
230 
231  for (size_t n = 0; n != num_nodes; ++n) {
232  tJac(i, j) += t_coords(i) * t_diff_n(j);
233  ++t_coords;
234  ++t_diff_n;
235  }
238 
239  if (type == MBTET)
240  vOlume *= G_TET_W1[0] / 6.;
241 
242 #ifndef NDEBUG
243  if (vOlume <= std::numeric_limits<double>::epsilon() || vOlume != vOlume) {
244  MOFEM_LOG("SELF", Sev::warning) << "Volume is zero " << vOlume;
245  MOFEM_LOG("SELF", Sev::warning) << "Coords " << coords << endl;
246  }
247 #endif
248 
250 }
251 
255  // Get coords at Gauss points
257 
258  auto &shape_functions = dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE);
259  double *shape_functions_ptr = &*shape_functions.data().begin();
260  const size_t nb_base_functions = shape_functions.size2();
261  const size_t nb_gauss_pts = gaussPts.size2();
262  coordsAtGaussPts.resize(nb_gauss_pts, 3, false);
263  coordsAtGaussPts.clear();
264  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_coords_at_gauss_ptr(
265  &coordsAtGaussPts(0, 0), &coordsAtGaussPts(0, 1),
266  &coordsAtGaussPts(0, 2));
268  shape_functions_ptr);
269  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
271  &coords[0], &coords[1], &coords[2]);
272  for (int bb = 0; bb != nb_base_functions; ++bb) {
273  t_coords_at_gauss_ptr(i) += t_coords(i) * t_shape_functions;
274  ++t_coords;
275  ++t_shape_functions;
276  };
277  ++t_coords_at_gauss_ptr;
278  }
280 }
281 
285 
288 
289  auto type = numeredEntFiniteElementPtr->getEntType();
290  auto dim_type = CN::Dimension(type);
291 
292  auto get_data_on_ents = [&](auto lower_dim, auto space) {
294  auto data = dataOnElement[space];
295  data->facesNodes = dataH1.facesNodes;
296  data->facesNodesOrder = dataH1.facesNodesOrder;
297  for (auto dd = dim_type; dd >= lower_dim; --dd) {
298  int nb_ents = moab::CN::NumSubEntities(type, dd);
299  for (int ii = 0; ii != nb_ents; ++ii) {
300  auto sub_ent_type = moab::CN::SubEntityType(type, dd, ii);
301  if ((dataH1.spacesOnEntities[sub_ent_type]).test(space)) {
302  auto &data_on_ent = data->dataOnEntities[sub_ent_type];
303  CHKERR getEntitySense(sub_ent_type, data_on_ent);
304  CHKERR getEntityDataOrder(sub_ent_type, space, data_on_ent);
305  data->spacesOnEntities[sub_ent_type].set(space);
306  }
307  }
308  }
310  };
311 
312  CHKERR get_data_on_ents(1, H1); // H1
313  CHKERR get_data_on_ents(1, HCURL); // Hcurl
314  CHKERR get_data_on_ents(2, HDIV); // Hdiv
315  CHKERR get_data_on_ents(3, L2); // L2
316 
318 }
319 
322 
323  if (numeredEntFiniteElementPtr->getEntType() == MBTET) {
324  // OK that code is not nice.
325  MatrixDouble new_diff_n;
326  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
328  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
329  EntitiesFieldData::EntData &data = dataH1.dataOnEntities[MBVERTEX][0];
330  if ((data.getDiffN(base).size1() == 4) &&
331  (data.getDiffN(base).size2() == 3)) {
332  const size_t nb_gauss_pts = gaussPts.size2();
333  const size_t nb_base_functions = 4;
334  new_diff_n.resize(nb_gauss_pts, 3 * nb_base_functions, false);
335  double *new_diff_n_ptr = &*new_diff_n.data().begin();
337  new_diff_n_ptr, &new_diff_n_ptr[1], &new_diff_n_ptr[2]);
338  double *t_diff_n_ptr = &*data.getDiffN(base).data().begin();
339  for (unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
341  t_diff_n_ptr, &t_diff_n_ptr[1], &t_diff_n_ptr[2]);
342  for (unsigned int bb = 0; bb != nb_base_functions; bb++) {
343  t_new_diff_n(i) = t_diff_n(i);
344  ++t_new_diff_n;
345  ++t_diff_n;
346  }
347  }
348  data.getDiffN(base).resize(new_diff_n.size1(), new_diff_n.size2(),
349  false);
350  data.getDiffN(base).swap(new_diff_n);
351  }
352  }
353  }
354 
355  if (dataH1.spacesOnEntities[MBVERTEX].test(H1))
357 
358  std::array<std::bitset<LASTSPACE>, 4> spaces_by_dim{
359 
360  std::bitset<LASTSPACE>(0), //
361  std::bitset<LASTSPACE>(0), //
362  std::bitset<LASTSPACE>(0), //
363  std::bitset<LASTSPACE>(0)
364 
365  };
366  for (auto type = MBEDGE; type != MBENTITYSET; ++type) {
367  spaces_by_dim[CN::Dimension(type)] |= dataH1.spacesOnEntities[type];
368  }
369 
370  if (
371 
372  spaces_by_dim[1].test(HCURL) || //
373  spaces_by_dim[2].test(HCURL) || //
374  spaces_by_dim[3].test(HCURL)
375 
376  ) {
379  }
380 
381  if (spaces_by_dim[2].test(HDIV) || spaces_by_dim[3].test(HDIV)) {
383  // Fix for tetrahedrons
384  if (numeredEntFiniteElementPtr->getEntType() == MBTET) {
385  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
386  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
387  for (auto t : {MBTRI, MBTET}) {
388  for (auto &d : dataHdiv.dataOnEntities[t]) {
389  d.getN(base) /= 6;
390  d.getDiffN(base) /= 6;
391  }
392  }
393  }
394  }
396  }
397 
398  if (spaces_by_dim[3].test(L2)) {
400  }
401 
403 }
404 
406  ForcesAndSourcesCore *ptr) {
408  if (!(ptrFE = dynamic_cast<VolumeElementForcesAndSourcesCore *>(ptr)))
409  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
410  "User operator and finite element do not work together");
412 }
413 
416 
417  const auto type = numeredEntFiniteElementPtr->getEntType();
419  switch (type) {
420  case MBTET:
422  boost::shared_ptr<BaseFunction>(new TetPolynomialBase(this));
423  break;
424  case MBHEX:
426  boost::shared_ptr<BaseFunction>(new HexPolynomialBase());
427  break;
428  default:
430  }
433  }
434 
438  if (gaussPts.size2() == 0)
443 
445 
446  // Iterate over operators
448 
450 }
451 
452 } // 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::VolumeElementForcesAndSourcesCore
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)
Definition: VolumeElementForcesAndSourcesCore.cpp:9
MoFEM::VolumeElementForcesAndSourcesCore::calculateVolumeAndJacobian
virtual MoFEMErrorCode calculateVolumeAndJacobian()
Calculate element volume and Jacobian.
Definition: VolumeElementForcesAndSourcesCore.cpp:196
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:95
MoFEM::ForcesAndSourcesCore::setGaussPts
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
Definition: ForcesAndSourcesCore.cpp:1961
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:96
MoFEM::Tools::diffShapeFunMBTET
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition: Tools.hpp:271
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
MoFEM::VolumeElementForcesAndSourcesCore::calculateCoordinatesAtGaussPts
virtual MoFEMErrorCode calculateCoordinatesAtGaussPts()
Calculate coordinate at integration points.
Definition: VolumeElementForcesAndSourcesCore.cpp:253
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:126
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:89
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:100
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::VolumeElementForcesAndSourcesCore::setIntegrationPts
virtual MoFEMErrorCode setIntegrationPts()
Set integration points.
Definition: VolumeElementForcesAndSourcesCore.cpp:22
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
diffN_MBHEX2x
#define diffN_MBHEX2x(y, z)
Definition: fem_tools.h:81
MoFEM::VolumeElementForcesAndSourcesCore::transformBaseFunctions
virtual MoFEMErrorCode transformBaseFunctions()
Transform base functions based on geometric element Jacobian.
Definition: VolumeElementForcesAndSourcesCore.cpp:320
convert.type
type
Definition: convert.py:64
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator::setPtrFE
MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr)
Definition: VolumeElementForcesAndSourcesCore.cpp:405
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:101
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:170
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:1955
diffN_MBHEX7x
#define diffN_MBHEX7x(y, z)
Definition: fem_tools.h:86
MoFEM::VolumeElementForcesAndSourcesCore::opSetInvJacH1
OpSetInvJacH1 opSetInvJacH1
Definition: VolumeElementForcesAndSourcesCore.hpp:93
MoFEM::VolumeElementForcesAndSourcesCore::tJac
FTensor::Tensor2< double *, 3, 3 > tJac
Definition: VolumeElementForcesAndSourcesCore.hpp:102
diffN_MBHEX7y
#define diffN_MBHEX7y(x, z)
Definition: fem_tools.h:94
MoFEM::VolumeElementForcesAndSourcesCore::operator()
MoFEMErrorCode operator()()
function is run for every finite element
Definition: VolumeElementForcesAndSourcesCore.cpp:414
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:98
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:738
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:90
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:94
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::VolumeElementForcesAndSourcesCore::getSpaceBaseAndOrderOnElement
virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement()
Determine approximation space and order of base functions.
Definition: VolumeElementForcesAndSourcesCore.cpp:283
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:383
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
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:103
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:985
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::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