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 
9 VolumeElementForcesAndSourcesCore::VolumeElementForcesAndSourcesCore(
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 
22 MoFEMErrorCode VolumeElementForcesAndSourcesCore::setIntegrationPts() {
23  MoFEMFunctionBegin;
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 = [&]() {
42  MoFEMFunctionBegin;
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  }
58  MoFEMFunctionReturn(0);
59  };
60 
61  auto calc_base_for_hex = [&]() {
62  MoFEMFunctionBegin;
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  }
105  MoFEMFunctionReturn(0);
106  };
107 
108  auto set_integration_pts_for_hex = [&]() {
109  MoFEMFunctionBegin;
110  CHKERR Tools::outerProductOfEdgeIntegrationPtsForHex(gaussPts, rule, rule,
111  rule);
112  MoFEMFunctionReturn(0);
113  };
114 
115  auto set_integration_pts_for_tet = [&]() {
116  MoFEMFunctionBegin;
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) {
122  SETERRQ2(mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
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  }
158  MoFEMFunctionReturn(0);
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 
193  MoFEMFunctionReturn(0);
194 }
195 
196 MoFEMErrorCode VolumeElementForcesAndSourcesCore::calculateVolumeAndJacobian() {
197  MoFEMFunctionBegin;
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 = [&]() {
205  return FTensor::Tensor1<FTensor::PackPtr<const double *, 3>, 3>(
206  &Tools::diffShapeFunMBTET[0], &Tools::diffShapeFunMBTET[1],
207  &Tools::diffShapeFunMBTET[2]);
208  };
209 
210  auto get_hex_t_diff_n = [&]() {
211  return FTensor::Tensor1<FTensor::PackPtr<const double *, 3>, 3>(
212  &Tools::diffShapeFunMBHEXAtCenter[0],
213  &Tools::diffShapeFunMBHEXAtCenter[1],
214  &Tools::diffShapeFunMBHEXAtCenter[2]);
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 
225  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_coords(
226  &coords[0], &coords[1], &coords[2]);
227  FTensor::Index<'i', 3> i;
228  FTensor::Index<'j', 3> j;
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  }
236  CHKERR determinantTensor3by3(tJac, vOlume);
237  CHKERR invertTensor3by3(tJac, vOlume, tInvJac);
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 
249  MoFEMFunctionReturn(0);
250 }
251 
252 MoFEMErrorCode
253 VolumeElementForcesAndSourcesCore::calculateCoordinatesAtGaussPts() {
254  MoFEMFunctionBegin;
255  // Get coords at Gauss points
256  FTensor::Index<'i', 3> i;
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));
267  FTensor::Tensor0<FTensor::PackPtr<double *, 1>> t_shape_functions(
268  shape_functions_ptr);
269  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
270  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_coords(
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  }
279  MoFEMFunctionReturn(0);
280 }
281 
282 MoFEMErrorCode
283 VolumeElementForcesAndSourcesCore::getSpaceBaseAndOrderOnElement() {
284  MoFEMFunctionBegin;
285 
286  CHKERR getSpacesAndBaseOnEntities(dataH1);
287  CHKERR getFaceNodes(dataH1);
288  // H1
289  if ((dataH1.spacesOnEntities[MBEDGE]).test(H1)) {
290  CHKERR getEntitySense<MBEDGE>(dataH1);
291  CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
292  }
293  if ((dataH1.spacesOnEntities[MBTRI]).test(H1)) {
294  CHKERR getEntitySense<MBTRI>(dataH1);
295  CHKERR getEntityDataOrder<MBTRI>(dataH1, H1);
296  }
297  if ((dataH1.spacesOnEntities[MBTET]).test(H1)) {
298  CHKERR getEntityDataOrder<MBTET>(dataH1, H1);
299  }
300  if ((dataH1.spacesOnEntities[MBQUAD]).test(H1)) {
301  CHKERR getEntitySense<MBQUAD>(dataH1);
302  CHKERR getEntityDataOrder<MBQUAD>(dataH1, H1);
303  }
304  if ((dataH1.spacesOnEntities[MBHEX]).test(H1)) {
305  CHKERR getEntityDataOrder<MBHEX>(dataH1, H1);
306  }
307  // Hcurl
308  if ((dataH1.spacesOnEntities[MBEDGE]).test(HCURL)) {
309  CHKERR getEntitySense<MBEDGE>(dataHcurl);
310  CHKERR getEntityDataOrder<MBEDGE>(dataHcurl, HCURL);
311  dataHcurl.spacesOnEntities[MBEDGE].set(HCURL);
312  }
313  if ((dataH1.spacesOnEntities[MBTRI]).test(HCURL)) {
314  dataHcurl.facesNodes = dataH1.facesNodes;
315  CHKERR getEntitySense<MBTRI>(dataHcurl);
316  CHKERR getEntityDataOrder<MBTRI>(dataHcurl, HCURL);
317  dataHcurl.spacesOnEntities[MBTRI].set(HCURL);
318  }
319  if ((dataH1.spacesOnEntities[MBTET]).test(HCURL)) {
320  CHKERR getEntityDataOrder<MBTET>(dataHcurl, HCURL);
321  dataHcurl.spacesOnEntities[MBTET].set(HCURL);
322  }
323  if ((dataH1.spacesOnEntities[MBQUAD]).test(HCURL)) {
324  dataHcurl.facesNodes = dataH1.facesNodes;
325  dataHcurl.facesNodesOrder = dataH1.facesNodesOrder;
326  CHKERR getEntitySense<MBQUAD>(dataHcurl);
327  CHKERR getEntityDataOrder<MBQUAD>(dataHcurl, HCURL);
328  dataHcurl.spacesOnEntities[MBQUAD].set(HCURL);
329  }
330  if ((dataH1.spacesOnEntities[MBHEX]).test(HCURL)) {
331  CHKERR getEntityDataOrder<MBHEX>(dataHcurl, HCURL);
332  dataHcurl.spacesOnEntities[MBHEX].set(HCURL);
333  }
334  // Hdiv
335  if ((dataH1.spacesOnEntities[MBTRI]).test(HDIV)) {
336  dataHdiv.facesNodes = dataH1.facesNodes;
337  CHKERR getEntitySense<MBTRI>(dataHdiv);
338  CHKERR getEntityDataOrder<MBTRI>(dataHdiv, HDIV);
339  dataHdiv.spacesOnEntities[MBTRI].set(HDIV);
340  }
341  if ((dataH1.spacesOnEntities[MBTET]).test(HDIV)) {
342  CHKERR getEntityDataOrder<MBTET>(dataHdiv, HDIV);
343  dataHdiv.spacesOnEntities[MBTET].set(HDIV);
344  }
345  if ((dataH1.spacesOnEntities[MBQUAD]).test(HDIV)) {
346  dataHdiv.facesNodes = dataH1.facesNodes;
347  dataHdiv.facesNodesOrder = dataH1.facesNodesOrder;
348  CHKERR getEntitySense<MBQUAD>(dataHdiv);
349  CHKERR getEntityDataOrder<MBQUAD>(dataHdiv, HDIV);
350  dataHdiv.spacesOnEntities[MBQUAD].set(HDIV);
351  }
352  if ((dataH1.spacesOnEntities[MBHEX]).test(HDIV)) {
353  CHKERR getEntityDataOrder<MBHEX>(dataHdiv, HDIV);
354  dataHdiv.spacesOnEntities[MBHEX].set(HDIV);
355  }
356  // L2
357  if ((dataH1.spacesOnEntities[MBTET]).test(L2)) {
358  CHKERR getEntityDataOrder<MBTET>(dataL2, L2);
359  dataL2.spacesOnEntities[MBTET].set(L2);
360  }
361  if ((dataH1.spacesOnEntities[MBHEX]).test(L2)) {
362  CHKERR getEntityDataOrder<MBHEX>(dataL2, L2);
363  dataL2.spacesOnEntities[MBHEX].set(L2);
364  }
365  MoFEMFunctionReturn(0);
366 }
367 
368 MoFEMErrorCode VolumeElementForcesAndSourcesCore::transformBaseFunctions() {
369  MoFEMFunctionBegin;
370 
371  if (numeredEntFiniteElementPtr->getEntType() == MBTET) {
372  // OK that code is not nice.
373  MatrixDouble new_diff_n;
374  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
375  FTensor::Index<'i', 3> i;
376  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
377  EntitiesFieldData::EntData &data = dataH1.dataOnEntities[MBVERTEX][0];
378  if ((data.getDiffN(base).size1() == 4) &&
379  (data.getDiffN(base).size2() == 3)) {
380  const size_t nb_gauss_pts = gaussPts.size2();
381  const size_t nb_base_functions = 4;
382  new_diff_n.resize(nb_gauss_pts, 3 * nb_base_functions, false);
383  double *new_diff_n_ptr = &*new_diff_n.data().begin();
384  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_new_diff_n(
385  new_diff_n_ptr, &new_diff_n_ptr[1], &new_diff_n_ptr[2]);
386  double *t_diff_n_ptr = &*data.getDiffN(base).data().begin();
387  for (unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
388  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_diff_n(
389  t_diff_n_ptr, &t_diff_n_ptr[1], &t_diff_n_ptr[2]);
390  for (unsigned int bb = 0; bb != nb_base_functions; bb++) {
391  t_new_diff_n(i) = t_diff_n(i);
392  ++t_new_diff_n;
393  ++t_diff_n;
394  }
395  }
396  data.getDiffN(base).resize(new_diff_n.size1(), new_diff_n.size2(),
397  false);
398  data.getDiffN(base).swap(new_diff_n);
399  }
400  }
401  }
402 
403  if (dataH1.spacesOnEntities[MBVERTEX].test(H1))
404  CHKERR opSetInvJacH1.opRhs(dataH1);
405  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
406  CHKERR opCovariantPiolaTransform.opRhs(dataHcurl);
407  CHKERR opSetInvJacHdivAndHcurl.opRhs(dataHcurl);
408  }
409  for (EntityType t = CN::TypeDimensionMap[2].first;
410  t <= CN::TypeDimensionMap[2].second; ++t) {
411  if (dataH1.spacesOnEntities[t].test(HDIV)) {
412  CHKERR opContravariantPiolaTransform.opRhs(dataHdiv);
413  // Fix for tetrahedrons
414  if (numeredEntFiniteElementPtr->getEntType() == MBTET) {
415  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
416  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
417  for (auto t : {MBTRI, MBTET}) {
418  for (auto &d : dataHdiv.dataOnEntities[t]) {
419  d.getN(base) /= 6;
420  d.getDiffN(base) /= 6;
421  }
422  }
423  }
424  }
425  CHKERR opSetInvJacHdivAndHcurl.opRhs(dataHdiv);
426  break;
427  }
428  }
429  for (EntityType t = CN::TypeDimensionMap[3].first;
430  t <= CN::TypeDimensionMap[3].second; ++t) {
431  if (dataH1.spacesOnEntities[t].test(L2)) {
432  CHKERR opSetInvJacH1.opRhs(dataL2);
433  break;
434  }
435  }
436 
437  MoFEMFunctionReturn(0);
438 }
439 
440 MoFEMErrorCode VolumeElementForcesAndSourcesCore::UserDataOperator::setPtrFE(
441  ForcesAndSourcesCore *ptr) {
442  MoFEMFunctionBeginHot;
443  if (!(ptrFE = dynamic_cast<VolumeElementForcesAndSourcesCore *>(ptr)))
444  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
445  "User operator and finite element do not work together");
446  MoFEMFunctionReturnHot(0);
447 }
448 
449 MoFEMErrorCode VolumeElementForcesAndSourcesCore::operator()() {
450  MoFEMFunctionBegin;
451 
452  const auto type = numeredEntFiniteElementPtr->getEntType();
453  if (type != lastEvaluatedElementEntityType) {
454  switch (type) {
455  case MBTET:
456  getElementPolynomialBase() =
457  boost::shared_ptr<BaseFunction>(new TetPolynomialBase());
458  break;
459  case MBHEX:
460  getElementPolynomialBase() =
461  boost::shared_ptr<BaseFunction>(new HexPolynomialBase());
462  break;
463  default:
464  MoFEMFunctionReturnHot(0);
465  }
466  CHKERR createDataOnElement(type);
467  lastEvaluatedElementEntityType = type;
468  }
469 
470  CHKERR calculateVolumeAndJacobian();
471  CHKERR getSpaceBaseAndOrderOnElement();
472  CHKERR setIntegrationPts();
473  if (gaussPts.size2() == 0)
474  MoFEMFunctionReturnHot(0);
475  CHKERR calculateCoordinatesAtGaussPts();
476  CHKERR calHierarchicalBaseFunctionsOnElement();
477  CHKERR calBernsteinBezierBaseFunctionsOnElement();
478 
479  CHKERR transformBaseFunctions();
480 
481  // Iterate over operators
482  CHKERR loopOverOperators();
483 
484  MoFEMFunctionReturn(0);
485 }
486 
487 } // namespace MoFEM
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
MoFEM::VolumeElementForcesAndSourcesCore::setIntegrationPts
virtual MoFEMErrorCode setIntegrationPts()
Set integration points.
Definition: VolumeElementForcesAndSourcesCore.cpp:22
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:47
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::ForcesAndSourcesCore::calBernsteinBezierBaseFunctionsOnElement
MoFEMErrorCode calBernsteinBezierBaseFunctionsOnElement()
Calculate Bernstein-Bezier base.
Definition: ForcesAndSourcesCore.cpp:1106
sdf_hertz.d
float d
Definition: sdf_hertz.py:5
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::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::ForcesAndSourcesCore::setGaussPts
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
Definition: ForcesAndSourcesCore.cpp:1926
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::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:1316
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:1086
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:50
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:46
zeta
double zeta
Viscous hardening.
Definition: plastic.cpp:177
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:1559
diffN_MBHEX5y
#define diffN_MBHEX5y(x, z)
Definition: fem_tools.h:92
MoFEM::VolumeElementForcesAndSourcesCore::VolumeElementForcesAndSourcesCore
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)
Definition: VolumeElementForcesAndSourcesCore.cpp:9
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:535
MoFEM::VolumeElementForcesAndSourcesCore::num_nodes
int num_nodes
Definition: VolumeElementForcesAndSourcesCore.hpp:100
MoFEM::VolumeElementForcesAndSourcesCore::calculateVolumeAndJacobian
virtual MoFEMErrorCode calculateVolumeAndJacobian()
Calculate element volume and Jacobian.
Definition: VolumeElementForcesAndSourcesCore.cpp:196
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:283
diffN_MBHEX2x
#define diffN_MBHEX2x(y, z)
Definition: fem_tools.h:81
convert.type
type
Definition: convert.py:64
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator::setPtrFE
MoFEMErrorCode setPtrFE(ForcesAndSourcesCore *ptr)
Definition: VolumeElementForcesAndSourcesCore.cpp:440
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:1305
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:1920
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
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:59
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:1511
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
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:56
MoFEM::ForcesAndSourcesCore::loopOverOperators
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
Definition: ForcesAndSourcesCore.cpp:1347
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:440
MoFEM::VolumeElementForcesAndSourcesCore::operator()
MoFEMErrorCode operator()()
function is run for every finite element
Definition: VolumeElementForcesAndSourcesCore.cpp:449
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:416
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:984
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
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:368
MoFEM::ForcesAndSourcesCore::mField
Interface & mField
Definition: ForcesAndSourcesCore.hpp:24
MoFEM::TetPolynomialBase
Calculate base functions on tetrahedral.
Definition: TetPolynomialBase.hpp:19
Generated by Doxygen 1.8.17 and hosted at University of Glasgow