v0.9.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 /* 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 namespace MoFEM {
22 
24  Interface &m_field, const EntityType type)
25  : ForcesAndSourcesCore(m_field),
26  meshPositionsFieldName("MESH_NODE_POSITIONS"), coords(12), jAc(3, 3),
27  invJac(3, 3), opSetInvJacH1(invJac),
28  opContravariantPiolaTransform(vOlume, jAc),
29  opCovariantPiolaTransform(invJac), opSetInvJacHdivAndHcurl(invJac),
30  opHOatGaussPoints(hoCoordsAtGaussPts, hoGaussPtsJac),
31  opSetHoInvJacH1(hoGaussPtsInvJac),
32  opHoContravariantTransform(hoGaussPtsDetJac, hoGaussPtsJac),
33  opHoCovariantTransform(hoGaussPtsInvJac),
34  opSetHoInvJacHdivAndHcurl(hoGaussPtsInvJac),
35  tJac(&jAc(0, 0), &jAc(0, 1), &jAc(0, 2), &jAc(1, 0), &jAc(1, 1),
36  &jAc(1, 2), &jAc(2, 0), &jAc(2, 1), &jAc(2, 2)),
37  tInvJac(&invJac(0, 0), &invJac(0, 1), &invJac(0, 2), &invJac(1, 0),
38  &invJac(1, 1), &invJac(1, 2), &invJac(2, 0), &invJac(2, 1),
39  &invJac(2, 2)) {
41  boost::shared_ptr<BaseFunction>(new TetPolynomialBase());
42 }
43 
46  int order_data = getMaxDataOrder();
47  int order_row = getMaxRowOrder();
48  int order_col = getMaxColOrder();
49  int rule = getRule(order_row, order_col, order_data);
50 
51  if (rule >= 0) {
52  if (rule < QUAD_3D_TABLE_SIZE) {
53  if (QUAD_3D_TABLE[rule]->dim != 3) {
54  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "wrong dimension");
55  }
56  if (QUAD_3D_TABLE[rule]->order < rule) {
58  "wrong order %d != %d", QUAD_3D_TABLE[rule]->order, rule);
59  }
60  size_t nb_gauss_pts = QUAD_3D_TABLE[rule]->npoints;
61  gaussPts.resize(4, nb_gauss_pts, false);
62  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[1], 4,
63  &gaussPts(0, 0), 1);
64  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[2], 4,
65  &gaussPts(1, 0), 1);
66  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[3], 4,
67  &gaussPts(2, 0), 1);
68  cblas_dcopy(nb_gauss_pts, QUAD_3D_TABLE[rule]->weights, 1, &gaussPts(3, 0),
69  1);
70  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4,
71  false);
72  double *shape_ptr =
73  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
74  cblas_dcopy(4 * nb_gauss_pts, QUAD_3D_TABLE[rule]->points, 1, shape_ptr, 1);
75  } else {
77  "rule > quadrature order %d < %d", rule, QUAD_3D_TABLE_SIZE);
78  }
79  } else {
80  CHKERR setGaussPts(order_row, order_col, order_data);
81  const size_t nb_gauss_pts = gaussPts.size2();
82  dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4,
83  false);
84  if (nb_gauss_pts > 0) {
86  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin(),
87  &gaussPts(0, 0), &gaussPts(1, 0), &gaussPts(2, 0), nb_gauss_pts);
88  }
89  }
91 }
92 
97  CHKERR mField.get_moab().get_connectivity(ent, conn, num_nodes, true);
98  CHKERR mField.get_moab().get_coords(conn, num_nodes, &*coords.data().begin());
103  &coords[0], &coords[1], &coords[2]);
106  jAc.clear();
107  for (auto n : {0, 1, 2, 3}) {
108  tJac(i, j) += t_coords(i) * t_diff_n(j);
109  ++t_coords;
110  ++t_diff_n;
111  }
114  vOlume *= G_TET_W1[0] / 6.;
116 }
117 
121  // Get coords at Gauss points
123 
124  double *shape_functions_ptr =
125  &*dataH1.dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
126  const size_t nb_gauss_pts = gaussPts.size2();
127  coordsAtGaussPts.resize(nb_gauss_pts, 3, false);
128  coordsAtGaussPts.clear();
129  FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_coords_at_gauss_ptr(
130  &coordsAtGaussPts(0, 0), &coordsAtGaussPts(0, 1),
131  &coordsAtGaussPts(0, 2));
133  shape_functions_ptr);
134  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
136  &coords[0], &coords[1], &coords[2]);
137  for (int bb = 0; bb != 4; ++bb) {
138  t_coords_at_gauss_ptr(i) += t_coords(i) * t_shape_functions;
139  ++t_coords;
140  ++t_shape_functions;
141  };
142  ++t_coords_at_gauss_ptr;
143  }
145 }
146 
150 
153  // H1
154  if ((dataH1.spacesOnEntities[MBEDGE]).test(H1)) {
155  CHKERR getEntitySense<MBEDGE>(dataH1);
156  CHKERR getEntityDataOrder<MBEDGE>(dataH1, H1);
157  }
158  if ((dataH1.spacesOnEntities[MBTRI]).test(H1)) {
159  CHKERR getEntitySense<MBTRI>(dataH1);
160  CHKERR getEntityDataOrder<MBTRI>(dataH1, H1);
161  }
162  if ((dataH1.spacesOnEntities[MBTET]).test(H1)) {
163  CHKERR getEntityDataOrder<MBTET>(dataH1, H1);
164  }
165  // Hcurl
166  if ((dataH1.spacesOnEntities[MBEDGE]).test(HCURL)) {
167  CHKERR getEntitySense<MBEDGE>(dataHcurl);
168  CHKERR getEntityDataOrder<MBEDGE>(dataHcurl, HCURL);
169  dataHcurl.spacesOnEntities[MBEDGE].set(HCURL);
170  }
171  if ((dataH1.spacesOnEntities[MBTRI]).test(HCURL)) {
172  CHKERR getEntitySense<MBTRI>(dataHcurl);
174  CHKERR getEntityDataOrder<MBTRI>(dataHcurl, HCURL);
175  dataHcurl.spacesOnEntities[MBTRI].set(HCURL);
176  }
177  if ((dataH1.spacesOnEntities[MBTET]).test(HCURL)) {
178  CHKERR getEntityDataOrder<MBTET>(dataHcurl, HCURL);
179  dataHcurl.spacesOnEntities[MBTET].set(HCURL);
180  }
181  // Hdiv
182  if ((dataH1.spacesOnEntities[MBTRI]).test(HDIV)) {
183  CHKERR getEntitySense<MBTRI>(dataHdiv);
185  CHKERR getEntityDataOrder<MBTRI>(dataHdiv, HDIV);
186  dataHdiv.spacesOnEntities[MBTRI].set(HDIV);
187  }
188  if ((dataH1.spacesOnEntities[MBTET]).test(HDIV)) {
189  CHKERR getEntityDataOrder<MBTET>(dataHdiv, HDIV);
190  dataHdiv.spacesOnEntities[MBTET].set(HDIV);
191  }
192  // L2
193  if ((dataH1.spacesOnEntities[MBTET]).test(L2)) {
194  CHKERR getEntityDataOrder<MBTET>(dataL2, L2);
195  dataL2.spacesOnEntities[MBTET].set(L2);
196  }
198 }
199 
202 
204  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
207  }
208  if (dataH1.spacesOnEntities[MBTRI].test(HDIV)) {
211  }
212  if (dataH1.spacesOnEntities[MBTET].test(L2)) {
214  }
215 
216  MatrixDouble new_diff_n;
217  for (int b = AINSWORTH_LEGENDRE_BASE; b != LASTBASE; b++) {
219  FieldApproximationBase base = static_cast<FieldApproximationBase>(b);
221  dataH1.dataOnEntities[MBVERTEX][0];
222  if ((data.getDiffN(base).size1() == 4) &&
223  (data.getDiffN(base).size2() == 3)) {
224  const size_t nb_gauss_pts = gaussPts.size2();
225  const size_t nb_base_functions = 4;
226  new_diff_n.resize(nb_gauss_pts, 3 * nb_base_functions, false);
227  double *new_diff_n_ptr = &*new_diff_n.data().begin();
229  new_diff_n_ptr, &new_diff_n_ptr[1], &new_diff_n_ptr[2]);
230  double *t_diff_n_ptr = &*data.getDiffN(base).data().begin();
231  for (unsigned int gg = 0; gg != nb_gauss_pts; gg++) {
233  t_diff_n_ptr, &t_diff_n_ptr[1], &t_diff_n_ptr[2]);
234  for (unsigned int bb = 0; bb != nb_base_functions; bb++) {
235  t_new_diff_n(i) = t_diff_n(i);
236  ++t_new_diff_n;
237  ++t_diff_n;
238  }
239  }
240  data.getDiffN(base).resize(new_diff_n.size1(), new_diff_n.size2(), false);
241  data.getDiffN(base).data().swap(new_diff_n.data());
242  }
243  }
244 
246 }
247 
250  if (dataPtr->get<FieldName_mi_tag>().find(meshPositionsFieldName) !=
251  dataPtr->get<FieldName_mi_tag>().end()) {
252  const Field *field_struture =
254  BitFieldId id = field_struture->getId();
255  if ((numeredEntFiniteElementPtr->getBitFieldIdData() & id).none()) {
256  SETERRQ(mField.get_comm(), MOFEM_NOT_FOUND,
257  "no MESH_NODE_POSITIONS in element data");
258  }
259 
261  if (dataH1.dataOnEntities[MBVERTEX][0].getFieldData().size() != 12) {
262  SETERRQ(mField.get_comm(), MOFEM_NOT_FOUND,
263  "no MESH_NODE_POSITIONS in element data or field has wrong "
264  "number of coefficients");
265  }
268  hoGaussPtsInvJac.resize(hoGaussPtsJac.size1(), hoGaussPtsJac.size2(),
269  false);
270  // Express Jacobian as tensor
272  &hoGaussPtsJac(0, 0), &hoGaussPtsJac(0, 1), &hoGaussPtsJac(0, 2),
273  &hoGaussPtsJac(0, 3), &hoGaussPtsJac(0, 4), &hoGaussPtsJac(0, 5),
274  &hoGaussPtsJac(0, 6), &hoGaussPtsJac(0, 7), &hoGaussPtsJac(0, 8));
276  &hoGaussPtsInvJac(0, 0), &hoGaussPtsInvJac(0, 1),
277  &hoGaussPtsInvJac(0, 2), &hoGaussPtsInvJac(0, 3),
278  &hoGaussPtsInvJac(0, 4), &hoGaussPtsInvJac(0, 5),
279  &hoGaussPtsInvJac(0, 6), &hoGaussPtsInvJac(0, 7),
280  &hoGaussPtsInvJac(0, 8));
281 
282  const size_t nb_gauss_pts = gaussPts.size2();
283  hoGaussPtsDetJac.resize(nb_gauss_pts, false);
285  // Calculate inverse and determinant
286  for (unsigned int gg = 0; gg != nb_gauss_pts; ++gg) {
287  CHKERR determinantTensor3by3(jac, det);
288  // if(det<0) {
289  // SETERRQ(mField.get_comm(),MOFEM_DATA_INCONSISTENCY,"Negative
290  // volume");
291  // }
292  CHKERR invertTensor3by3(jac, det, inv_jac);
293  ++jac;
294  ++inv_jac;
295  ++det;
296  }
297  } else {
298  hoCoordsAtGaussPts.resize(0, 0, false);
299  hoGaussPtsInvJac.resize(0, 0, false);
300  hoGaussPtsDetJac.resize(0, false);
301  }
303 }
304 
308  if (hoCoordsAtGaussPts.size1() > 0) {
309  // Transform derivatives of base functions and apply Piola transformation
310  // if needed.
311 
313  if (dataH1.spacesOnEntities[MBTET].test(L2)) {
315  }
316  if (dataH1.spacesOnEntities[MBTRI].test(HDIV)) {
319  }
320  if (dataH1.spacesOnEntities[MBEDGE].test(HCURL)) {
323  }
324  }
326 }
327 
329  getDivergenceOfHDivBaseFunctions(int side, EntityType type,
331  int gg, VectorDouble &div) {
333 
334  int nb_dofs = data.getFieldData().size();
335  if (nb_dofs == 0)
337 
338  if (data.getSpace() != HDIV && data.getSpace() != HCURL) {
340  "This function should be used for HDIV used but is used with %s",
341  FieldSpaceNames[data.getSpace()]);
342  }
343 
344  if ((unsigned int)nb_dofs != data.getDiffN().size2() / 9) {
346  "Data inositency, wrong number of dofs = %s "
347  "%d != %d/9",
348  FieldSpaceNames[data.getSpace()], nb_dofs,
349  data.getDiffN().size2());
350  }
351 
352  div.resize(nb_dofs, false);
353 
354  FTensor::Tensor0<double *> t_div(&*div.data().begin());
355  const double *grad_ptr = &data.getDiffN()(gg, 0);
357  grad_ptr, &grad_ptr[HVEC1_1], &grad_ptr[HVEC2_2]);
358  for (int dd = 0; dd < nb_dofs; dd++) {
359  t_div = t_grad_base(0) + t_grad_base(1) + t_grad_base(2);
360  ++t_div;
361  ++t_grad_base;
362  }
363 
365 }
366 
368  getCurlOfHCurlBaseFunctions(int side, EntityType type,
369  DataForcesAndSourcesCore::EntData &data, int gg,
370  MatrixDouble &curl) {
372 
373  int nb_dofs = data.getFieldData().size();
374  if (nb_dofs == 0)
376 
377  if (data.getSpace() != HDIV && data.getSpace() != HCURL) {
378  SETERRQ1(getVolumeFE()->mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
379  "This function should be used for primarily for HCURL"
380  " but will work with HDIV used but is used with %s",
381  FieldSpaceNames[data.getSpace()]);
382  }
383 
384  if ((unsigned int)nb_dofs != data.getDiffN().size2() / 9) {
385  SETERRQ3(getVolumeFE()->mField.get_comm(), MOFEM_DATA_INCONSISTENCY,
386  "Data insistency, wrong number of dofs = %s "
387  "%d != %d/9",
388  FieldSpaceNames[data.getSpace()], nb_dofs,
389  data.getDiffN().size2());
390  }
391 
392  curl.resize(nb_dofs, 3, false);
394  &curl(0, 0), &curl(0, 1), &curl(0, 2));
395  const double *grad_ptr = &data.getDiffN()(gg, 0);
396 
398  grad_ptr, &grad_ptr[HVEC0_1], &grad_ptr[HVEC0_2], &grad_ptr[HVEC1_0],
399  &grad_ptr[HVEC1_1], &grad_ptr[HVEC1_2], &grad_ptr[HVEC2_0],
400  &grad_ptr[HVEC2_1], &grad_ptr[HVEC2_2]);
401  for (int dd = 0; dd != nb_dofs; ++dd) {
402  t_curl(0) = t_grad_base(2, 1) - t_grad_base(1, 2);
403  t_curl(1) = t_grad_base(0, 2) - t_grad_base(2, 0);
404  t_curl(2) = t_grad_base(1, 0) - t_grad_base(0, 1);
405  ++t_curl;
406  ++t_grad_base;
407  }
408 
410 }
411 
412 } // namespace MoFEM
virtual MoFEMErrorCode calculateCoordinatesAtGaussPts()
Calculate coordinate at integration points.
DataForcesAndSourcesCore & dataL2
MoFEMErrorCode getCurlOfHCurlBaseFunctions(int side, EntityType type, DataForcesAndSourcesCore::EntData &data, int gg, MatrixDouble &curl)
Get curl of base functions at integration point.
int getMaxRowOrder() const
Get max order of approximation for field in rows.
virtual MoFEMErrorCode calculateVolumeAndJacobian()
Calculate element volume and Jacobian.
#define QUAD_3D_TABLE_SIZE
Definition: quad.h:186
MatrixDouble gaussPts
Matrix of integration points.
field with continuous normal traction
Definition: definitions.h:173
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:415
virtual MoFEMErrorCode getSpaceBaseAndOrderOnElement()
Determine approximation space and order of base functions.
virtual moab::Interface & get_moab()=0
virtual int getRule(int order_row, int order_col, int order_data)
another variant of getRule
virtual MoFEMErrorCode setIntegrationPts()
Set integration points.
Calculate base functions on tetrahedral.
Provide data structure for (tensor) field approximation.The Field is intended to provide support for ...
DataForcesAndSourcesCore & dataH1
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
FieldSpace & getSpace()
Get field space.
VolumeElementForcesAndSourcesCoreBase * getVolumeFE() const
return pointer to Generic Volume Finite Element object
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:74
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
static constexpr std::array< double, 12 > diffShapeFunMBTET
Definition: Tools.hpp:143
MoFEMErrorCode getFaceTriNodes(DataForcesAndSourcesCore &data) const
Get nodes on triangles.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
virtual MoFEMErrorCode transformBaseFunctions()
Transform base functions based on geometric element Jacobian.
boost::shared_ptr< const FEDofEntity_multiIndex > dataPtr
Pointer to finite element data dofs.
OpGetDataAndGradient< 3, 3 > opHOatGaussPoints
higher order geometry data at Gauss pts
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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:396
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
MoFEMErrorCode getNodesFieldData(DataForcesAndSourcesCore &data, const std::string &field_name) const
Get data on nodes.
virtual MoFEMErrorCode calculateHoJacobian()
Calculate Jacobian for HO geometry.
DataForcesAndSourcesCore & dataHcurl
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:144
virtual MoFEMErrorCode transformHoBaseFunctions()
Transform base functions based on ho-geometry element Jacobian.
static const double G_TET_W1[]
Definition: fem_tools.h:495
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:146
static const char *const FieldSpaceNames[]
Definition: definitions.h:178
field with continuous tangents
Definition: definitions.h:172
virtual MoFEMErrorCode setGaussPts(int order_row, int order_col, int order_data)
set user specific integration rule
std::bitset< BITFIELDID_SIZE > BitFieldId
Field Id.
Definition: Types.hpp:53
#define CHKERR
Inline error check.
Definition: definitions.h:596
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
int order
Definition: quad.h:28
MultiIndex Tag for field name.
VolumeElementForcesAndSourcesCoreBase(Interface &m_field, const EntityType type=MBTET)
DataForcesAndSourcesCore & dataHdiv
Data on single entity (This is passed as argument to DataOperator::doWork)
std::array< std::bitset< LASTSPACE >, MBMAXTYPE > spacesOnEntities
spaces on entity types
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:407
virtual MPI_Comm & get_comm() const =0
continuous field
Definition: definitions.h:171
virtual const Field * get_field_structure(const std::string &name)=0
get field structure
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
const VectorDouble & getFieldData() const
get dofs values
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
MoFEMErrorCode getDivergenceOfHDivBaseFunctions(int side, EntityType type, DataForcesAndSourcesCore::EntData &data, int gg, VectorDouble &div)
Get divergence of base functions at integration point.
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:448
int getMaxDataOrder() const
Get max order of approximation for data fields.
auto & getElementPolynomialBase()
Get the Entity Polynomial Base object.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
field with C-1 continuity
Definition: definitions.h:174