v0.8.23
FaceElementForcesAndSourcesCore.cpp
Go to the documentation of this file.
1 /** \file FaceElementForcesAndSourcesCore.cpp
2 
3 \brief Implementation of face 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)
36  : ForcesAndSourcesCore(m_field),
37  meshPositionsFieldName("MESH_NODE_POSITIONS"),
38  opHOCoordsAndNormals(hoCoordsAtGaussPts, normalsAtGaussPts,
39  tangentOneAtGaussPts, tangentTwoAtGaussPts),
40  opContravariantTransform(nOrmal, normalsAtGaussPts),
41  opCovariantTransform(nOrmal, normalsAtGaussPts, tangentOne,
42  tangentOneAtGaussPts, tangentTwo,
43  tangentTwoAtGaussPts) {
45  boost::shared_ptr<BaseFunction>(new TriPolynomialBase());
46 }
47 
50  const string &fe_name, VolumeElementForcesAndSourcesCoreOnSide &method) {
52 
53  const EntityHandle ent = getNumeredEntFiniteElementPtr()->getEnt();
54  const Problem *problem_ptr = getFEMethod()->problemPtr;
55  Range adjacent_volumes;
56  CHKERR getFaceFE()->mField.getInterface<BitRefManager>()->getAdjacenciesAny(
57  ent, 3, adjacent_volumes);
58  typedef NumeredEntFiniteElement_multiIndex::index<
59  Composite_Name_And_Ent_mi_tag>::type FEByComposite;
60  FEByComposite &numered_fe =
62 
63  method.feName = fe_name;
64 
65  CHKERR method.setFaceFEPtr(getFaceFE());
67  CHKERR method.copyKsp(*getFEMethod());
68  CHKERR method.copySnes(*getFEMethod());
69  CHKERR method.copyTs(*getFEMethod());
70 
71  CHKERR method.preProcess();
72 
73  int nn = 0;
74  method.loopSize = adjacent_volumes.size();
75  for (Range::iterator vit = adjacent_volumes.begin();
76  vit != adjacent_volumes.end(); vit++) {
77  FEByComposite::iterator miit =
78  numered_fe.find(boost::make_tuple(fe_name, *vit));
79  if (miit != numered_fe.end()) {
80  method.nInTheLoop = nn++;
81  method.numeredEntFiniteElementPtr = *miit;
82  method.dataFieldEntsPtr = (*miit)->sPtr->data_field_ents_view;
83  method.rowFieldEntsPtr = (*miit)->sPtr->row_field_ents_view;
84  method.colFieldEntsPtr = (*miit)->sPtr->col_field_ents_view;
85  method.dataPtr = (*miit)->sPtr->data_dofs;
86  method.rowPtr = (*miit)->rows_dofs;
87  method.colPtr = (*miit)->cols_dofs;
88  CHKERR method();
89  }
90  }
91 
92  CHKERR method.postProcess();
93 
95 }
96 
100  CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
101  coords.resize(num_nodes * 3, false);
102  CHKERR mField.get_moab().get_coords(conn, num_nodes, &*coords.data().begin());
103  double diff_n[6];
104  CHKERR ShapeDiffMBTRI(diff_n);
105 
106  nOrmal.resize(3, false);
107  tangentOne.resize(3, false);
108  tangentTwo.resize(3, false);
110  &coords[0], &coords[1], &coords[2]);
112  &nOrmal[0], &nOrmal[1], &nOrmal[2]);
114  &tangentOne[0], &tangentOne[1], &tangentOne[2]);
116  &tangentTwo[0], &tangentTwo[1], &tangentTwo[2]);
118  &diff_n[1]);
119 
125  t_t1(i) = 0;
126  t_t2(i) = 0;
127  for (int nn = 0; nn != 3; ++nn) {
128  t_t1(i) += t_coords(i) * t_diff(N0);
129  t_t2(i) += t_coords(i) * t_diff(N1);
130  ++t_coords;
131  ++t_diff;
132  }
133  t_normal(j) = FTensor::levi_civita(i, j, k) * t_t1(k) * t_t2(i);
134  aRea = sqrt(t_normal(i) * t_normal(i)) / 2.;
135 
137 }
138 
141  // Set integration points
142  int order_data = getMaxDataOrder();
143  int order_row = getMaxRowOrder();
144  int order_col = getMaxColOrder();
145  int rule = getRule(order_row, order_col, order_data);
146 
147  if (rule >= 0) {
148  if (rule < QUAD_2D_TABLE_SIZE) {
149  if (QUAD_2D_TABLE[rule]->dim != 2) {
150  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong dimension");
151  }
152  if (QUAD_2D_TABLE[rule]->order < rule) {
153  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
154  "wrong order %d != %d", QUAD_2D_TABLE[rule]->order, rule);
155  }
157  gaussPts.resize(3, nbGaussPts, false);
158  cblas_dcopy(nbGaussPts, &QUAD_2D_TABLE[rule]->points[1], 3,
159  &gaussPts(0, 0), 1);
160  cblas_dcopy(nbGaussPts, &QUAD_2D_TABLE[rule]->points[2], 3,
161  &gaussPts(1, 0), 1);
162  cblas_dcopy(nbGaussPts, QUAD_2D_TABLE[rule]->weights, 1, &gaussPts(2, 0),
163  1);
164  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nbGaussPts, 3,
165  false);
166  double *shape_ptr =
167  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
168  cblas_dcopy(3 * nbGaussPts, QUAD_2D_TABLE[rule]->points, 1, shape_ptr, 1);
169  } else {
170  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
171  "rule > quadrature order %d < %d", rule, QUAD_2D_TABLE_SIZE);
172  nbGaussPts = 0;
173  }
174  } else {
175  // If rule is negative, set user defined integration points
176  CHKERR setGaussPts(order_row, order_col, order_data);
177  nbGaussPts = gaussPts.size2();
178  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nbGaussPts, 3,
179  false);
180  if (nbGaussPts) {
182  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin(),
183  &gaussPts(0, 0), &gaussPts(1, 0), nbGaussPts);
184  }
185  }
187 }
188 
192  // Get spaces order/base and sense of entities.
193 
195 
196  // H1
197  if (dataH1.spacesOnEntities[MBEDGE].test(H1)) {
198  CHKERR getEntitySense<MBEDGE>(dataH1);
199  CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
200  }
201  if (dataH1.spacesOnEntities[MBTRI].test(H1)) {
202  CHKERR getEntitySense<MBTRI>(dataH1);
203  CHKERR getEntityDataOrder<MBTRI>(dataH1, H1);
204  }
205 
206  // Hcurl
207  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
208  CHKERR getEntitySense<MBEDGE>(dataHcurl);
209  CHKERR getEntityDataOrder<MBEDGE>(dataHcurl, HCURL);
210  dataHcurl.spacesOnEntities[MBEDGE].set(HCURL);
211  }
212  if (dataH1.spacesOnEntities[MBTRI].test(HCURL)) {
213  CHKERR getEntitySense<MBTRI>(dataHcurl);
214  CHKERR getEntityDataOrder<MBTRI>(dataHcurl, HCURL);
215  dataHcurl.spacesOnEntities[MBTRI].set(HCURL);
216  }
217 
218  // Hdiv
219  if (dataH1.spacesOnEntities[MBTRI].test(HDIV)) {
220  CHKERR getEntitySense<MBTRI>(dataHdiv);
221  CHKERR getEntityDataOrder<MBTRI>(dataHdiv, HDIV);
222  dataHdiv.spacesOnEntities[MBTRI].set(HDIV);
223  }
224 
225  // L2
226  if (dataH1.spacesOnEntities[MBTRI].test(L2)) {
227  CHKERR getEntitySense<MBTRI>(dataL2);
228  CHKERR getEntityDataOrder<MBTRI>(dataL2, L2);
229  dataL2.spacesOnEntities[MBTRI].set(L2);
230  }
231 
233 }
234 
238 
239  double *shape_functions =
240  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
241  coordsAtGaussPts.resize(nbGaussPts, 3, false);
242  for (int gg = 0; gg != nbGaussPts; ++gg) {
243  for (int dd = 0; dd != 3; ++dd) {
244  coordsAtGaussPts(gg, dd) =
245  cblas_ddot(3, &shape_functions[3 * gg], 1, &coords[dd], 3);
246  }
247  }
249 }
250 
253  // Check if field for high-order geometry is set and if it is set calculate
254  // higher-order normals and face tangent vectors.
255  if (dataPtr->get<FieldName_mi_tag>().find(meshPositionsFieldName) !=
256  dataPtr->get<FieldName_mi_tag>().end()) {
257 
258  const Field *field_struture =
260  BitFieldId id = field_struture->getId();
261 
262  if ((numeredEntFiniteElementPtr->getBitFieldIdData() & id).none()) {
263  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_FOUND,
264  "no MESH_NODE_POSITIONS in element data");
265  }
266 
267  // Calculate normal for high-order geometry
268 
274 
275  } else {
276  hoCoordsAtGaussPts.resize(0, 0, false);
277  normalsAtGaussPts.resize(0, 0, false);
278  tangentOneAtGaussPts.resize(0, 0, false);
279  tangentTwoAtGaussPts.resize(0, 0, false);
280  }
282 }
283 
286 
287  if (numeredEntFiniteElementPtr->getEntType() != MBTRI)
290 
291  // Calculate normal and tangent vectors for face geometry given by 3 nodes.
294 
296  if (nbGaussPts == 0)
298 
301 
302  dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).resize(3, 2, false);
303  std::copy(Tools::diffShapeFunMBTRI.begin(), Tools::diffShapeFunMBTRI.end(),
304  dataH1.dataOnEntities[MBVERTEX][0].getDiffN(NOBASE).data().begin());
305 
306  /// Use the some node base
310 
311  // Apply Piola transform to HDiv and HCurl spaces, uses previously calculated
312  // faces normal and tangent vectors.
313  if (dataH1.spacesOnEntities[MBTRI].test(HDIV)) {
315  }
316  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
317  CHKERR opCovariantTransform.opRhs(data_curl);
318  }
319 
320  // Iterate over operators
322 
324 }
325 
327 OpCalculateInvJacForFace::doWork(int side, EntityType type,
329 
331 
332  if (getNumeredEntFiniteElementPtr()->getEntType() != MBTRI) {
333  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
334  "This operator can be used only with element which is triangle");
335  }
336 
337  if (type == MBVERTEX) {
338  VectorDouble &coords = getCoords();
339  double *coords_ptr = &*coords.data().begin();
340  double j00 = 0, j01 = 0, j10 = 0, j11 = 0;
341 
342  // this is triangle, derivative of nodal shape functions is constant.
343  // So only need to do one node.
344  for (auto n : {0, 1, 2}) {
345  j00 += coords_ptr[3 * n + 0] * Tools::diffShapeFunMBTRI[2 * n + 0];
346  j01 += coords_ptr[3 * n + 0] * Tools::diffShapeFunMBTRI[2 * n + 1];
347  j10 += coords_ptr[3 * n + 1] * Tools::diffShapeFunMBTRI[2 * n + 0];
348  j11 += coords_ptr[3 * n + 1] * Tools::diffShapeFunMBTRI[2 * n + 1];
349  }
350  const double det = j00 * j11 - j01 * j10;
351 
352  invJac.resize(2, 2, false);
353  invJac(0, 0) = j11 / det;
354  invJac(0, 1) = -j01 / det;
355  invJac(1, 0) = -j10 / det;
356  invJac(1, 1) = j00 / det;
357  }
358 
359  doVertices = true;
360  doEdges = false;
361  doQuads = false;
362  doTris = false;
363  doTets = false;
364  doPrisms = false;
365 
367 }
368 
370 OpSetInvJacH1ForFace::doWork(int side, EntityType type,
373 
374  if (getNumeredEntFiniteElementPtr()->getEntType() != MBTRI &&
375  getNumeredEntFiniteElementPtr()->getEntType() != MBQUAD) {
376  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
377  "This operator can be used only with element which is triangle");
378  }
379 
380  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
381 
382  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
383 
384  unsigned int nb_functions = data.getN(base).size2();
385  if (nb_functions) {
386  unsigned int nb_gauss_pts = data.getN(base).size1();
387  diffNinvJac.resize(nb_gauss_pts, 2 * nb_functions, false);
388 
389  if (type != MBVERTEX) {
390  if (nb_functions != data.getDiffN(base).size2() / 2) {
391  SETERRQ2(
392  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
393  "data inconsistency nb_functions != data.diffN.size2()/2 ( %u != "
394  "%u/2 )",
395  nb_functions, data.getDiffN(base).size2());
396  }
397  }
398 
400  t_inv_jac(0, 0) = invJac(0, 0);
401  t_inv_jac(0, 1) = invJac(0, 1);
402  t_inv_jac(1, 0) = invJac(1, 0);
403  t_inv_jac(1, 1) = invJac(1, 1);
404 
405  switch (type) {
406  case MBVERTEX:
407  case MBEDGE:
408  case MBTRI:
409  case MBQUAD: {
414  &diffNinvJac(0, 0), &diffNinvJac(0, 1));
416  &data.getDiffN(base)(0, 0), &data.getDiffN(base)(0, 1));
417  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
418  for (unsigned int dd = 0; dd != nb_functions; ++dd) {
419  t_diff_n(i) = t_inv_jac(k, i) * t_diff_n_ref(k);
420  ++t_diff_n;
421  ++t_diff_n_ref;
422  }
423  }
424  data.getDiffN(base).data().swap(diffNinvJac.data());
425  } break;
426  default:
427  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "not implemented");
428  }
429  }
430  }
431 
433 }
434 
436 OpSetInvJacHcurlFace::doWork(int side, EntityType type,
439 
440  if (type != MBEDGE && type != MBTRI)
442 
443  if (getNumeredEntFiniteElementPtr()->getEntType() != MBTRI) {
444  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
445  "This operator can be used only with element which is triangle");
446  }
447 
449  &invJac(0, 0), &invJac(0, 1), &invJac(1, 0), &invJac(1, 1));
450 
454 
455  for (int b = AINSWORTH_LEGENDRE_BASE; b != USER_BASE; b++) {
456 
457  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
458 
459  const unsigned int nb_base_functions = data.getDiffN(base).size2() / 6;
460  if (nb_base_functions) {
461  const unsigned int nb_gauss_pts = data.getDiffN(base).size1();
462 
463  diffHcurlInvJac.resize(nb_gauss_pts, data.getDiffN(base).size2(), false);
464 
465  auto t_diff_n = data.getFTensor2DiffN<3, 2>(base);
466  double *t_inv_diff_n_ptr = &*diffHcurlInvJac.data().begin();
468  t_inv_diff_n_ptr, &t_inv_diff_n_ptr[HVEC0_1],
469  &t_inv_diff_n_ptr[HVEC1_0], &t_inv_diff_n_ptr[HVEC1_1],
470  &t_inv_diff_n_ptr[HVEC2_0], &t_inv_diff_n_ptr[HVEC2_1]);
471 
472  for (unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
473  for (unsigned int bb = 0; bb != nb_base_functions; bb++) {
474  t_inv_diff_n(i, j) = t_diff_n(i, k) * t_inv_jac(k, j);
475  ++t_diff_n;
476  ++t_inv_diff_n;
477  }
478  }
479 
480  data.getDiffN(base).data().swap(diffHcurlInvJac.data());
481  }
482  }
483 
485 }
486 
487 } // namespace MoFEM
MoFEMErrorCode copyBasicMethod(const BasicMethod &basic)
Copy data from other base method to this base method.
Definition: LoopMethods.cpp:99
DataForcesAndSourcesCore & dataL2
int getMaxRowOrder() const
Get max order of approximation for field in rows.
std::bitset< LASTSPACE > spacesOnEntities[MBMAXTYPE]
spaces on entity types
virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement()
Determine approximation space and order of base functions.
MoFEMErrorCode calculateBaseFunctionsOnElement()
Calculate base functions.
MatrixDouble gaussPts
Matrix of integration points.
field with continuous normal traction
Definition: definitions.h:172
user implemented approximation base
Definition: definitions.h:153
virtual moab::Interface & get_moab()=0
data structure for finite element entityIt keeps that about indices of degrees of freedom,...
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 loopSize
local number oe methods to process
virtual MoFEMErrorCode setIntegrationPts()
Set integration points.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:500
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
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
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.
boost::shared_ptr< const FieldEntity_vector_view > colFieldEntsPtr
Pointer to finite element field entities column view.
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:476
const boost::shared_ptr< DataForcesAndSourcesCore > dataOnElement[LASTSPACE]
Entity data on element entity rows fields.
std::string feName
Name of finite element.
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
boost::shared_ptr< const FieldEntity_vector_view > rowFieldEntsPtr
Pointer to finite element field entities row view.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:507
keeps basic data about problemThis is low level structure with information about problem,...
virtual const MatrixDouble & getN(const FieldApproximationBase base) const
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
boost::shared_ptr< const FEDofEntity_multiIndex > dataPtr
Pointer to finite element data dofs.
MoFEMErrorCode loopSideVolumes(const string &fe_name, VolumeElementForcesAndSourcesCoreOnSide &method)
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
boost::ptr_vector< EntData > dataOnEntities[MBMAXTYPE]
const Problem * problemPtr
raw pointer to problem
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: cblas_ddot.c:12
#define QUAD_2D_TABLE_SIZE
Definition: quad.h:174
static QUAD *const QUAD_2D_TABLE[]
Definition: quad.h:175
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
MoFEMErrorCode copyKsp(const KspMethod &ksp)
copy data form another method
Definition: LoopMethods.cpp:32
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode copyTs(const TSMethod &ts)
Copy TS solver data.
Definition: LoopMethods.cpp:75
OpSetCovariantPiolaTransformOnTriangle opCovariantTransform
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
DataForcesAndSourcesCore & dataHcurl
virtual MoFEMErrorCode preProcess()
function is run at the beginning of loop
PetscErrorCode ShapeDiffMBTRI(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:206
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
FieldApproximationBase
approximation base
Definition: definitions.h:142
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
const BitFieldId & getId() const
Get unique field id.
Managing BitRefLevels.
virtual MoFEMErrorCode calculateHoNormal()
Calculate normal on curved elements.
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:144
std::string meshPositionsFieldName
Name of the field with geometry.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
field with continuous tangents
Definition: definitions.h:171
static constexpr std::array< double, 6 > diffShapeFunMBTRI
Definition: Tools.hpp:83
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)
int nInTheLoop
number currently of processed method
MoFEMErrorCode getEntityFieldData(DataForcesAndSourcesCore &data, const std::string &field_name, const EntityType type_lo=MBVERTEX, const EntityType type_hi=MBPOLYHEDRON) const
PetscErrorCode ShapeMBTRI(double *N, const double *X, const double *Y, const int G_DIM)
calculate shape functions on triangle
Definition: fem_tools.c:194
Calculate base functions on triangle.
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
MultiIndex Tag for field name.
MoFEMErrorCode copySnes(const SnesMethod &snes)
Copy snes data.
Definition: LoopMethods.cpp:53
DataForcesAndSourcesCore & dataHdiv
OpSetContravariantPiolaTransformOnTriangle opContravariantTransform
Data on single entity (This is passed as argument to DataOperator::doWork)
boost::shared_ptr< const FieldEntity_multiIndex_spaceType_view > dataFieldEntsPtr
Pointer to finite element field entities data view.
structure to get information form mofem into DataForcesAndSourcesCore
virtual MoFEMErrorCode calculateCoordinatesAtGaussPts()
Calculate coordinate at integration points.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:406
virtual MoFEMErrorCode postProcess()
function is run at the end of loop
virtual MoFEMErrorCode calculateAreaAndNormal()
Calculate element area and normal of the face.
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
boost::shared_ptr< NumeredEntFiniteElement_multiIndex > numeredFiniteElements
store finite elements
const FaceElementForcesAndSourcesCore * getFaceFE()
return pointer to Generic Triangle Finite Element object
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode getSpacesAndBaseOnEntities(DataForcesAndSourcesCore &data) const
Get field approximation space and base on entities.
boost::shared_ptr< const FENumeredDofEntity_multiIndex > colPtr
Pointer to finite element columns dofs view.
int getMaxColOrder() const
Get max order of approximation for field in columns.
int getMaxDataOrder() const
Get max order of approximation for data fields.
const int order
approximation order
boost::shared_ptr< const FENumeredDofEntity_multiIndex > rowPtr
Pointer to finite element rows dofs view.
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
MoFEMErrorCode setFaceFEPtr(const FaceElementForcesAndSourcesCore *face_fe_ptr)
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use
field with C-1 continuity
Definition: definitions.h:173