v0.8.23
VolumeElementForcesAndSourcesCore.cpp
Go to the documentation of this file.
1 /** \file VolumeElementForcesAndSourcesCore.cpp
2 
3 \brief Implementation of volume element
4 
5 */
6 
7 /* This file is part of MoFEM.
8  * MoFEM is free software: you can redistribute it and/or modify it under
9  * the terms of the GNU Lesser General Public License as published by the
10  * Free Software Foundation, either version 3 of the License, or (at your
11  * option) any later version.
12  *
13  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  * License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20 
21 #ifdef __cplusplus
22 extern "C" {
23 #endif
24 #include <cblas.h>
25 #include <lapack_wrap.h>
26 // #include <gm_rule.h>
27 #include <quad.h>
28 #ifdef __cplusplus
29 }
30 #endif
31 
32 namespace MoFEM {
33 
35  Interface &m_field, const EntityType type)
36  : ForcesAndSourcesCore(m_field), coords(12), jAc(3, 3), invJac(3, 3),
37  opSetInvJacH1(invJac), opContravariantPiolaTransform(vOlume, jAc),
38  opCovariantPiolaTransform(invJac), opSetInvJacHdivAndHcurl(invJac),
39  meshPositionsFieldName("MESH_NODE_POSITIONS"),
40  opHOatGaussPoints(hoCoordsAtGaussPts, hoGaussPtsJac),
41  opSetHoInvJacH1(hoGaussPtsInvJac),
42  opHoContravariantTransform(hoGaussPtsDetJac, hoGaussPtsJac),
43  opHoCovariantTransform(hoGaussPtsInvJac),
44  opSetHoInvJacHdivAndHcurl(hoGaussPtsInvJac),
45  tJac(&jAc(0, 0), &jAc(0, 1), &jAc(0, 2), &jAc(1, 0), &jAc(1, 1),
46  &jAc(1, 2), &jAc(2, 0), &jAc(2, 1), &jAc(2, 2)),
47  tInvJac(&invJac(0, 0), &invJac(0, 1), &invJac(0, 2), &invJac(1, 0),
48  &invJac(1, 1), &invJac(1, 2), &invJac(2, 0), &invJac(2, 1),
49  &invJac(2, 2)) {
51  boost::shared_ptr<BaseFunction>(new TetPolynomialBase());
52 }
53 
56  int order_data = getMaxDataOrder();
57  int order_row = getMaxRowOrder();
58  int order_col = getMaxColOrder();
59  int rule = getRule(order_row, order_col, order_data);
60 
61  if (rule >= 0) {
62  if (rule < QUAD_3D_TABLE_SIZE) {
63  if (QUAD_3D_TABLE[rule]->dim != 3) {
64  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "wrong dimension");
65  }
66  if (QUAD_3D_TABLE[rule]->order < rule) {
68  "wrong order %d != %d", QUAD_3D_TABLE[rule]->order, rule);
69  }
71  gaussPts.resize(4, nbGaussPts, false);
72  cblas_dcopy(nbGaussPts, &QUAD_3D_TABLE[rule]->points[1], 4,
73  &gaussPts(0, 0), 1);
74  cblas_dcopy(nbGaussPts, &QUAD_3D_TABLE[rule]->points[2], 4,
75  &gaussPts(1, 0), 1);
76  cblas_dcopy(nbGaussPts, &QUAD_3D_TABLE[rule]->points[3], 4,
77  &gaussPts(2, 0), 1);
78  cblas_dcopy(nbGaussPts, QUAD_3D_TABLE[rule]->weights, 1, &gaussPts(3, 0),
79  1);
80  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nbGaussPts, 4,
81  false);
82  double *shape_ptr =
83  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
84  cblas_dcopy(4 * nbGaussPts, QUAD_3D_TABLE[rule]->points, 1, shape_ptr, 1);
85  } else {
87  "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
88  nbGaussPts = 0;
89  }
90  } else {
91  CHKERR setGaussPts(order_row, order_col, order_data);
92  nbGaussPts = gaussPts.size2();
93  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nbGaussPts, 4,
94  false);
95  if (nbGaussPts > 0) {
97  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin(),
98  &gaussPts(0, 0), &gaussPts(1, 0), &gaussPts(2, 0), nbGaussPts);
99  }
100  }
102 }
103 
107  CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
108  CHKERR mField.get_moab().get_coords(conn, num_nodes, &*coords.data().begin());
113  &coords[0], &coords[1], &coords[2]);
116  jAc.clear();
117  for (auto n : {0, 1, 2, 3}) {
118  tJac(i, j) += t_coords(i) * t_diff_n(j);
119  ++t_coords;
120  ++t_diff_n;
121  }
124  vOlume *= G_TET_W1[0] / 6.;
126 }
127 
131  // Get coords at Gauss points
133 
134  double *shape_functions_ptr =
135  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
136  coordsAtGaussPts.resize(nbGaussPts, 3, false);
137  coordsAtGaussPts.clear();
138  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_coords_at_gauss_ptr(
139  &coordsAtGaussPts(0, 0), &coordsAtGaussPts(0, 1),
140  &coordsAtGaussPts(0, 2));
142  shape_functions_ptr);
143  for (unsigned int gg = 0; gg < nbGaussPts; gg++) {
145  &coords[0], &coords[1], &coords[2]);
146  for (int bb = 0; bb < 4; bb++) {
147  t_coords_at_gauss_ptr(i) += t_coords(i) * t_shape_functions;
148  ++t_coords;
149  ++t_shape_functions;
150  };
151  ++t_coords_at_gauss_ptr;
152  }
154 }
155 
159 
162  // H1
163  if ((dataH1.spacesOnEntities[MBEDGE]).test(H1)) {
164  CHKERR getEntitySense<MBEDGE>(dataH1);
165  CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
166  }
167  if ((dataH1.spacesOnEntities[MBTRI]).test(H1)) {
168  CHKERR getEntitySense<MBTRI>(dataH1);
169  CHKERR getEntityDataOrder<MBTRI>(dataH1, H1);
170  }
171  if ((dataH1.spacesOnEntities[MBTET]).test(H1)) {
172  CHKERR getEntityDataOrder<MBTET>(dataH1, H1);
173  }
174  // Hcurl
175  if ((dataH1.spacesOnEntities[MBEDGE]).test(HCURL)) {
176  CHKERR getEntitySense<MBEDGE>(dataHcurl);
177  CHKERR getEntityDataOrder<MBEDGE>(dataHcurl, HCURL);
178  dataHcurl.spacesOnEntities[MBEDGE].set(HCURL);
179  }
180  if ((dataH1.spacesOnEntities[MBTRI]).test(HCURL)) {
181  CHKERR getEntitySense<MBTRI>(dataHcurl);
183  CHKERR getEntityDataOrder<MBTRI>(dataHcurl, HCURL);
184  dataHcurl.spacesOnEntities[MBTRI].set(HCURL);
185  }
186  if ((dataH1.spacesOnEntities[MBTET]).test(HCURL)) {
187  CHKERR getEntityDataOrder<MBTET>(dataHcurl, HCURL);
188  dataHcurl.spacesOnEntities[MBTET].set(HCURL);
189  }
190  // Hdiv
191  if ((dataH1.spacesOnEntities[MBTRI]).test(HDIV)) {
192  CHKERR getEntitySense<MBTRI>(dataHdiv);
194  CHKERR getEntityDataOrder<MBTRI>(dataHdiv, HDIV);
195  dataHdiv.spacesOnEntities[MBTRI].set(HDIV);
196  }
197  if ((dataH1.spacesOnEntities[MBTET]).test(HDIV)) {
198  CHKERR getEntityDataOrder<MBTET>(dataHdiv, HDIV);
199  dataHdiv.spacesOnEntities[MBTET].set(HDIV);
200  }
201  // L2
202  if ((dataH1.spacesOnEntities[MBTET]).test(L2)) {
203  CHKERR getEntityDataOrder<MBTET>(dataL2, L2);
204  dataL2.spacesOnEntities[MBTET].set(L2);
205  }
207 }
208 
211 
213  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
216  }
217  if (dataH1.spacesOnEntities[MBTRI].test(HDIV)) {
220  }
221  if (dataH1.spacesOnEntities[MBTET].test(L2)) {
223  }
225 }
226 
229  if (dataPtr->get<FieldName_mi_tag>().find(meshPositionsFieldName) !=
230  dataPtr->get<FieldName_mi_tag>().end()) {
231  const Field *field_struture =
233  BitFieldId id = field_struture->getId();
234  if ((numeredEntFiniteElementPtr->getBitFieldIdData() & id).none()) {
235  SETERRQ(mField.get_comm(), MOFEM_NOT_FOUND,
236  "no MESH_NODE_POSITIONS in element data");
237  }
238 
240  if (dataH1.dataOnEntities[MBVERTEX][0].getFieldData().size() != 12) {
241  SETERRQ(mField.get_comm(), MOFEM_NOT_FOUND,
242  "no MESH_NODE_POSITIONS in element data or field has wrong "
243  "number of coefficients");
244  }
247  hoGaussPtsInvJac.resize(hoGaussPtsJac.size1(), hoGaussPtsJac.size2(),
248  false);
249  // Express Jacobian as tensor
251  &hoGaussPtsJac(0, 0), &hoGaussPtsJac(0, 1), &hoGaussPtsJac(0, 2),
252  &hoGaussPtsJac(0, 3), &hoGaussPtsJac(0, 4), &hoGaussPtsJac(0, 5),
253  &hoGaussPtsJac(0, 6), &hoGaussPtsJac(0, 7), &hoGaussPtsJac(0, 8));
255  &hoGaussPtsInvJac(0, 0), &hoGaussPtsInvJac(0, 1),
256  &hoGaussPtsInvJac(0, 2), &hoGaussPtsInvJac(0, 3),
257  &hoGaussPtsInvJac(0, 4), &hoGaussPtsInvJac(0, 5),
258  &hoGaussPtsInvJac(0, 6), &hoGaussPtsInvJac(0, 7),
259  &hoGaussPtsInvJac(0, 8));
260  hoGaussPtsDetJac.resize(nbGaussPts, false);
262  // Calculate inverse and determinant
263  for (unsigned int gg = 0; gg != nbGaussPts; gg++) {
264  CHKERR determinantTensor3by3(jac, det);
265  // if(det<0) {
266  // SETERRQ(mField.get_comm(),MOFEM_DATA_INCONSISTENCY,"Negative
267  // volume");
268  // }
269  CHKERR invertTensor3by3(jac, det, inv_jac);
270  ++jac;
271  ++inv_jac;
272  ++det;
273  }
274  } else {
275  hoCoordsAtGaussPts.resize(0, 0, false);
276  hoGaussPtsInvJac.resize(0, 0, false);
277  hoGaussPtsDetJac.resize(0, false);
278  }
280 }
281 
284  if (hoCoordsAtGaussPts.size1() > 0) {
285  // Transform derivatives of base functions and apply Piola transformation
286  // if needed.
287 
289  if (dataH1.spacesOnEntities[MBTET].test(L2)) {
291  }
292  if (dataH1.spacesOnEntities[MBTRI].test(HDIV)) {
295  }
296  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
299  }
300  }
302 }
303 
306 
307  if (numeredEntFiniteElementPtr->getEntType() != MBTET)
310 
314  if (nbGaussPts == 0)
319 
320  try {
321  MatrixDouble new_diff_n;
322  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
324  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
326  dataH1.dataOnEntities[MBVERTEX][0];
327  if ((data.getDiffN(base).size1() == 4) &&
328  (data.getDiffN(base).size2() == 3)) {
329  const unsigned int nb_base_functions = 4;
330  new_diff_n.resize(nbGaussPts, 3 * nb_base_functions, false);
331  double *new_diff_n_ptr = &*new_diff_n.data().begin();
333  new_diff_n_ptr, &new_diff_n_ptr[1], &new_diff_n_ptr[2]);
334  double *t_diff_n_ptr = &*data.getDiffN(base).data().begin();
335  for (unsigned int gg = 0; gg != nbGaussPts; gg++) {
337  t_diff_n_ptr, &t_diff_n_ptr[1], &t_diff_n_ptr[2]);
338  for (unsigned int bb = 0; bb != nb_base_functions; bb++) {
339  t_new_diff_n(i) = t_diff_n(i);
340  ++t_new_diff_n;
341  ++t_diff_n;
342  }
343  }
344  data.getDiffN(base).resize(new_diff_n.size1(), new_diff_n.size2(),
345  false);
346  data.getDiffN(base).data().swap(new_diff_n.data());
347  }
348  }
349  }
350  CATCH_ERRORS;
351 
354 
355  // Iterate over operators
357 
359 }
360 
362  getDivergenceOfHDivBaseFunctions(int side, EntityType type,
364  int gg, VectorDouble &div) {
366 
367  int nb_dofs = data.getFieldData().size();
368  if (nb_dofs == 0)
370 
371  if (data.getSpace() != HDIV && data.getSpace() != HCURL) {
373  "This function should be used for HDIV used but is used with %s",
374  FieldSpaceNames[data.getSpace()]);
375  }
376 
377  if ((unsigned int)nb_dofs != data.getDiffN().size2() / 9) {
379  "Data inositency, wrong number of dofs = %s "
380  "%d != %d/9",
381  FieldSpaceNames[data.getSpace()], nb_dofs,
382  data.getDiffN().size2());
383  }
384 
385  div.resize(nb_dofs, false);
386 
387  FTensor::Tensor0<double *> t_div(&*div.data().begin());
388  const double *grad_ptr = &data.getDiffN()(gg, 0);
390  grad_ptr, &grad_ptr[HVEC1_1], &grad_ptr[HVEC2_2]);
391  for (int dd = 0; dd < nb_dofs; dd++) {
392  t_div = t_grad_base(0) + t_grad_base(1) + t_grad_base(2);
393  ++t_div;
394  ++t_grad_base;
395  }
396 
398 }
399 
401  getCurlOfHCurlBaseFunctions(int side, EntityType type,
402  DataForcesAndSourcesCore::EntData &data, int gg,
403  MatrixDouble &curl) {
405 
406  int nb_dofs = data.getFieldData().size();
407  if (nb_dofs == 0)
409 
410  if (data.getSpace() != HDIV && data.getSpace() != HCURL) {
411  SETERRQ1(getVolumeFE()->mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
412  "This function should be used for primarily for HCURL"
413  " but will work with HDIV used but is used with %s",
414  FieldSpaceNames[data.getSpace()]);
415  }
416 
417  if ((unsigned int)nb_dofs != data.getDiffN().size2() / 9) {
418  SETERRQ3(getVolumeFE()->mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
419  "Data insistency, wrong number of dofs = %s "
420  "%d != %d/9",
421  FieldSpaceNames[data.getSpace()], nb_dofs,
422  data.getDiffN().size2());
423  }
424 
425  curl.resize(nb_dofs, 3, false);
427  &curl(0, 0), &curl(0, 1), &curl(0, 2));
428  const double *grad_ptr = &data.getDiffN()(gg, 0);
429 
431  grad_ptr, &grad_ptr[HVEC0_1], &grad_ptr[HVEC0_2], &grad_ptr[HVEC1_0],
432  &grad_ptr[HVEC1_1], &grad_ptr[HVEC1_2], &grad_ptr[HVEC2_0],
433  &grad_ptr[HVEC2_1], &grad_ptr[HVEC2_2]);
434  for (int dd = 0; dd != nb_dofs; ++dd) {
435  t_curl(0) = t_grad_base(2, 1) - t_grad_base(1, 2);
436  t_curl(1) = t_grad_base(0, 2) - t_grad_base(2, 0);
437  t_curl(2) = t_grad_base(1, 0) - t_grad_base(0, 1);
438  ++t_curl;
439  ++t_grad_base;
440  }
441 
443 }
444 
447  if (faceFEPtr == NULL) {
448  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
449  "Pointer to face element is not set");
450  }
451  const EntityHandle face_entity =
452  faceFEPtr->numeredEntFiniteElementPtr->getEnt();
453  SideNumber_multiIndex &side_table = const_cast<SideNumber_multiIndex &>(
454  numeredEntFiniteElementPtr->getSideNumberTable());
455  SideNumber_multiIndex::nth_index<0>::type::iterator sit =
456  side_table.get<0>().find(face_entity);
457  if (sit == side_table.get<0>().end()) {
458  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
459  "Face can not be found on volume element");
460  }
461  faceSense = (*sit)->sense;
462  faceSideNumber = (*sit)->side_number;
463  fill(tetConnMap, &tetConnMap[4], -1);
464  for (int nn = 0; nn != 3; nn++) {
465  faceConnMap[nn] =
466  std::distance(conn, find(conn, &conn[4], faceFEPtr->conn[nn]));
467  tetConnMap[faceConnMap[nn]] = nn;
468  if (faceConnMap[nn] > 3) {
469  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
470  "No common node on face and element can not be found");
471  }
472  }
473  oppositeNode =
474  std::distance(tetConnMap, find(tetConnMap, &tetConnMap[4], -1));
475  const int nb_gauss_pts = faceFEPtr->gaussPts.size2();
476  gaussPts.resize(4, nb_gauss_pts, false);
477  gaussPts.clear();
478  DataForcesAndSourcesCore &dataH1_on_face = *faceFEPtr->dataOnElement[H1];
479  const MatrixDouble &face_shape_funtions =
480  dataH1_on_face.dataOnEntities[MBVERTEX][0].getN(NOBASE);
481  const double tet_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
482  for (int gg = 0; gg != nb_gauss_pts; gg++) {
483  gaussPts(0, gg) =
484  face_shape_funtions(gg, 0) * tet_coords[3 * faceConnMap[0] + 0] +
485  face_shape_funtions(gg, 1) * tet_coords[3 * faceConnMap[1] + 0] +
486  face_shape_funtions(gg, 2) * tet_coords[3 * faceConnMap[2] + 0];
487  gaussPts(1, gg) =
488  face_shape_funtions(gg, 0) * tet_coords[3 * faceConnMap[0] + 1] +
489  face_shape_funtions(gg, 1) * tet_coords[3 * faceConnMap[1] + 1] +
490  face_shape_funtions(gg, 2) * tet_coords[3 * faceConnMap[2] + 1];
491  gaussPts(2, gg) =
492  face_shape_funtions(gg, 0) * tet_coords[3 * faceConnMap[0] + 2] +
493  face_shape_funtions(gg, 1) * tet_coords[3 * faceConnMap[1] + 2] +
494  face_shape_funtions(gg, 2) * tet_coords[3 * faceConnMap[2] + 2];
495  gaussPts(3, gg) = faceFEPtr->gaussPts(2, gg);
496  }
498 }
499 
500 VectorDouble &
502  return getFaceFEPtr()->nOrmal;
503 }
504 
507  return getFaceFEPtr()->normalsAtGaussPts;
508 }
509 
512  return getFaceFEPtr()->coordsAtGaussPts;
513 }
514 
515 /** \brief if higher order geometry return normals at Gauss pts.
516  *
517  * \param gg gauss point number
518  */
519 ublas::matrix_row<MatrixDouble>
521  const int gg) {
522  return ublas::matrix_row<MatrixDouble>(getNormalsAtGaussPts(), gg);
523 }
524 
525 } // namespace MoFEM
DataForcesAndSourcesCore & dataL2
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:337
int getMaxRowOrder() const
Get max order of approximation for field in rows.
std::bitset< LASTSPACE > spacesOnEntities[MBMAXTYPE]
spaces on entity types
MoFEMErrorCode calculateBaseFunctionsOnElement()
Calculate base functions.
#define QUAD_3D_TABLE_SIZE
Definition: quad.h:186
virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement()
Determine approximation space and order of base functions.
MatrixDouble gaussPts
Matrix of integration points.
field with continuous normal traction
Definition: definitions.h:172
virtual moab::Interface & get_moab()=0
OpGetDataAndGradient< 3, 3 > opHOatGaussPoints
higher order geometry data at Gauss pts
data structure for finite element entityIt keeps that about indices of degrees of freedom,...
Calculate base functions on tetrahedral.
MoFEMErrorCode createDataOnElement()
Create a entity data on element object.
Provide data structure for (tensor) field approximation.The Field is intended to provide support for ...
DataForcesAndSourcesCore & dataH1
int npoints
Definition: quad.h:29
void cblas_dcopy(const int N, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_dcopy.c:11
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:77
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant.
Definition: Templates.hpp:357
MoFEMErrorCode getNodesFieldData(const boost::string_ref field_name, FEDofEntity_multiIndex &dofs, VectorDouble &nodes_data, VectorDofs &nodes_dofs, FieldSpace &space, FieldApproximationBase &base) const
Get field data on nodes.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition: Tools.hpp:116
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)
multi_index_container< boost::shared_ptr< SideNumber >, indexed_by< hashed_unique< member< SideNumber, EntityHandle, &SideNumber::ent > >, ordered_non_unique< composite_key< SideNumber, const_mem_fun< SideNumber, EntityType, &SideNumber::getEntType >, member< SideNumber, char, &SideNumber::side_number > > >, ordered_non_unique< const_mem_fun< SideNumber, EntityType, &SideNumber::getEntType > > > > SideNumber_multiIndex
SideNumber_multiIndex for SideNumber.
MoFEMErrorCode getFaceTriNodes(DataForcesAndSourcesCore &data) const
Get nodes on triangles.
MoFEMErrorCode operator()()
function is run for every finite element
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507
virtual MoFEMErrorCode calculateCoordinatesAtGaussPts()
Calculate coordinate at integration points.
boost::shared_ptr< const FEDofEntity_multiIndex > dataPtr
Pointer to finite element data dofs.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MatrixDouble & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
MoFEMErrorCode getCurlOfHCurlBaseFunctions(int side, EntityType type, DataForcesAndSourcesCore::EntData &data, int gg, MatrixDouble &curl)
Get curl of base functions at integration point.
boost::ptr_vector< EntData > dataOnEntities[MBMAXTYPE]
OpSetHoContravariantPiolaTransform opHoContravariantTransform
MoFEMErrorCode loopOverOperators()
Iterate user data operators.
virtual const MatrixDouble & getDiffN(const FieldApproximationBase base) const
get derivatives of base functions
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
DataForcesAndSourcesCore & dataHcurl
unsigned int nbGaussPts
Number of integration points.
MoFEMErrorCode getDivergenceOfHDivBaseFunctions(int side, EntityType type, DataForcesAndSourcesCore::EntData &data, int gg, VectorDouble &div)
Get divergence of base functions at integration point.
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
virtual MoFEMErrorCode calculateHoJacobian()
Calculate Jacobian for HO geometry.
FieldApproximationBase
approximation base
Definition: definitions.h:142
virtual MoFEMErrorCode calculateVolumeAndJacobian()
Calculate element volume and Jacobian.
static const double G_TET_W1[]
Definition: fem_tools.h:495
virtual MoFEMErrorCode transformBaseFunctions()
Transform base functions based on geometric element Jacobian.
const VectorDouble & getFieldData() const
get dofs values
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
const BitFieldId & getId() const
Get unique field id.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:144
static const char *const FieldSpaceNames[]
Definition: definitions.h:177
field with continuous tangents
Definition: definitions.h:171
FieldSpace & getSpace()
Get field space.
virtual MoFEMErrorCode setIntegrationPts()
Set integration points.
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:53
#define CHKERR
Inline error check.
Definition: definitions.h:595
virtual MoFEMErrorCode opRhs(DataForcesAndSourcesCore &data, const bool do_vertices, const bool do_edges, const bool do_quads, const bool do_tris, const bool do_tets, const bool do_prisms, const bool error_if_no_base=true)
MoFEMErrorCode getEntityFieldData(DataForcesAndSourcesCore &data, const std::string &field_name, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
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:176
int order
Definition: quad.h:28
const VolumeElementForcesAndSourcesCore * getVolumeFE()
return pointer to Generic Volume Finite Element object
MultiIndex Tag for field name.
DataForcesAndSourcesCore & dataHdiv
Data on single entity (This is passed as argument to DataOperator::doWork)
OpSetContravariantPiolaTransform opContravariantPiolaTransform
structure to get information form mofem into DataForcesAndSourcesCore
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:432
virtual MPI_Comm & get_comm() const =0
continuous field
Definition: definitions.h:170
virtual const Field * get_field_structure(const std::string &name)=0
get field structure
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:76
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
virtual MoFEMErrorCode transformHoBaseFunctions()
Transform base functions based on ho-geometry element Jacobian.
MoFEMErrorCode getSpacesAndBaseOnEntities(DataForcesAndSourcesCore &data) const
Get field approximation space and base on entities.
int getMaxColOrder() const
Get max order of approximation for field in columns.
static QUAD *const QUAD_3D_TABLE[]
Definition: quad.h:187
int getMaxDataOrder() const
Get max order of approximation for data fields.
const int order
approximation order
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
field with C-1 continuity
Definition: definitions.h:173