v0.13.1
DataOperators.hpp
Go to the documentation of this file.
1/** \file DataOperators.hpp
2
3 Number of structures preforming operations on integration points.
4
5 For example:
6 - calculate Jacobian
7 - calculate Piola-Transform on shape functions
8 - calculate normals and tangent vectors at integration points on triangle
9
10 */
11
12/* This file is part of MoFEM.
13 * MoFEM is free software: you can redistribute it and/or modify it under
14 * the terms of the GNU Lesser General Public License as published by the
15 * Free Software Foundation, either version 3 of the License, or (at your
16 * option) any later version.
17 *
18 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
19 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
21 * License for more details.
22 *
23 * You should have received a copy of the GNU Lesser General Public
24 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
25
26#ifndef __DATAOPERATORS_HPP
27#define __DATAOPERATORS_HPP
28
29using namespace boost::numeric;
30
31namespace MoFEM {
32
33/** \brief base operator to do operations at Gauss Pt. level
34 * \ingroup mofem_forces_and_sources
35 */
37
38 DataOperator(const bool symm = true);
39
40 virtual ~DataOperator() = default;
41
42 using DoWorkLhsHookFunType = boost::function<MoFEMErrorCode(
43 DataOperator *op_ptr, int row_side, int col_side, EntityType row_type,
44 EntityType col_type, EntitiesFieldData::EntData &row_data,
46
48
49 /** \brief Operator for bi-linear form, usually to calculate values on
50 * left hand side
51 */
52 virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
53 EntityType col_type,
57 if (doWorkLhsHook) {
58 CHKERR doWorkLhsHook(this, row_side, col_side, row_type, col_type,
59 row_data, col_data);
60 } else {
61 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
62 "doWork function is not implemented for this operator");
63 }
65 }
66
67 virtual MoFEMErrorCode opLhs(EntitiesFieldData &row_data,
68 EntitiesFieldData &col_data);
69
70 using DoWorkRhsHookFunType = boost::function<MoFEMErrorCode(
71 DataOperator *op_ptr, int side, EntityType type,
73
75
76 /** \brief Operator for linear form, usually to calculate values on right hand
77 * side
78 */
82 if (doWorkRhsHook) {
83 CHKERR doWorkRhsHook(this, side, type, data);
84 } else {
85 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
86 "doWork function is not implemented for this operator");
87 }
89 }
90
92 const bool error_if_no_base = false);
93
94 bool sYmm; ///< If true assume that matrix is symmetric structure
95
96 std::array<bool, MBMAXTYPE>
97 doEntities; ///< If true operator is executed for entity.
98
99 // Deprecated variables. Use doEntities instead. I keep them for back
100 // compatibility with some older modules. It will be removed in some future.
101
102 bool &doVertices; ///< \deprectaed If false skip vertices
103 bool &doEdges; ///< \deprectaed If false skip edges
104 bool &doQuads; ///< \deprectaed
105 bool &doTris; ///< \deprectaed
106 bool &doTets; ///< \deprectaed
107 bool &doPrisms; ///< \deprectaed
108
109 /**
110 * \brief Get if operator uses symmetry of DOFs or not
111 *
112 * If symmetry is used, only not repeating combinations of entities are
113 * looped. For an example pair of (Vertex, Edge_0) and (Edge_0, Vertex) will
114 * calculate the same matrices only transposed. Implementing that this can be
115 * exploited by integrating only one pair.
116 *
117 * @return true if symmetry
118 */
119 inline bool getSymm() const { return sYmm; }
120
121 /// set if operator is executed taking in account symmetry
122 inline void setSymm() { sYmm = true; }
123
124 /// unset if operator is executed for non symmetric problem
125 inline void unSetSymm() { sYmm = false; }
126
127private:
128 template <bool Symm>
129 inline MoFEMErrorCode opLhs(EntitiesFieldData &row_data,
130 EntitiesFieldData &col_data);
131
132 template <bool ErrorIfNoBase>
134 const std::array<bool, MBMAXTYPE> &do_entities);
135};
136
137/**
138 * \brief Transform local reference derivatives of shape function to global
139 derivatives
140
141 * \ingroup mofem_forces_and_sources
142 */
144
148
150 : tInvJac(&inv_jac(0, 0), &inv_jac(0, 1), &inv_jac(0, 2), &inv_jac(1, 0),
151 &inv_jac(1, 1), &inv_jac(1, 2), &inv_jac(2, 0), &inv_jac(2, 1),
152 &inv_jac(2, 2)) {}
153
157};
158
159/**
160 * \brief brief Transform local reference derivatives of shape function to
161 global derivatives
162
163 * \ingroup mofem_forces_and_sources
164 */
166
171
173 : tInvJac(&inv_jac(0, 0), &inv_jac(0, 1), &inv_jac(0, 2), &inv_jac(1, 0),
174 &inv_jac(1, 1), &inv_jac(1, 2), &inv_jac(2, 0), &inv_jac(2, 1),
175 &inv_jac(2, 2)) {}
176
180};
181
182/** \brief apply contravariant (Piola) transfer to Hdiv space
183
184Contravariant Piola transformation
185\f[
186\psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
187\left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
188=
189\frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
190\f]
191
192* \ingroup mofem_forces_and_sources
193
194*/
196
197 double &vOlume;
198
203
205 : vOlume(volume),
206 // jAc(jac),
207 tJac(&jac(0, 0), &jac(0, 1), &jac(0, 2), &jac(1, 0), &jac(1, 1),
208 &jac(1, 2), &jac(2, 0), &jac(2, 1), &jac(2, 2)) {}
209
212
215};
216
217/** \brief apply covariant transfer to Hcurl space
218
219Contravariant Piola transformation
220\f[
221\psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
222\left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
223=
224\frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
225\f]
226
227
228* \ingroup mofem_forces_and_sources
229*/
231
236
238 : tInvJac(&inv_jac(0, 0), &inv_jac(0, 1), &inv_jac(0, 2), &inv_jac(1, 0),
239 &inv_jac(1, 1), &inv_jac(1, 2), &inv_jac(2, 0), &inv_jac(2, 1),
240 &inv_jac(2, 2)) {}
241
244
247};
248
249/** \brief Get field values and gradients at Gauss points
250 * \ingroup mofem_forces_and_sources
251 */
252template <int RANK, int DIM> struct OpGetDataAndGradient : public DataOperator {
253
256
258 MatrixDouble &data_grad_at_gauss_pt)
259 : dataAtGaussPts(data_at_gauss_pt),
260 dataGradAtGaussPts(data_grad_at_gauss_pt) {}
261
262 /**
263 * Return tensor associated with matrix storing values
264 */
265 template <int R>
267 THROW_MESSAGE("Not implemented");
268 }
269
270 /**
271 * Return tensor associated with matrix storing gradient values
272 */
273 template <int R, int D>
275 THROW_MESSAGE("Not implemented");
276 }
277
278 /**
279 * \brief Calculate gradient and values at integration points
280 * @param side side of entity on element
281 * @param type type of entity
282 * @param data data stored on entity (dofs values, dofs indices, etc.)
283 * @return error code
284 */
288 const int nb_base_functions = data.getN().size2();
289 bool constant_diff = false;
290 if (type == MBVERTEX && data.getDiffN().size1() * data.getDiffN().size2() ==
291 DIM * nb_base_functions) {
292 constant_diff = true;
293 }
294 const int nb_dofs = data.getFieldData().size();
295 for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
296 double *data_ptr, *n_ptr, *diff_n_ptr;
297 n_ptr = &data.getN()(gg, 0);
298 if (constant_diff) {
299 diff_n_ptr = &data.getDiffN()(0, 0);
300 } else {
301 diff_n_ptr = &data.getDiffN()(gg, 0);
302 }
303 data_ptr = &*data.getFieldData().data().begin();
304 for (int rr = 0; rr < RANK; rr++, data_ptr++) {
305 dataAtGaussPts(gg, rr) +=
306 cblas_ddot(nb_dofs / RANK, n_ptr, 1, data_ptr, RANK);
307 double *diff_n_ptr2 = diff_n_ptr;
308 for (unsigned int dd = 0; dd < DIM; dd++, diff_n_ptr2++) {
309 dataGradAtGaussPts(gg, DIM * rr + dd) +=
310 cblas_ddot(nb_dofs / RANK, diff_n_ptr2, DIM, data_ptr, RANK);
311 }
312 }
313 }
315 }
316
319
321
322 if (data.getFieldData().size() == 0) {
324 }
325
326 unsigned int nb_dofs = data.getFieldData().size();
327 if (nb_dofs == 0)
329
330 if (nb_dofs % RANK != 0) {
331 SETERRQ4(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
332 "data inconsistency, type %d, side %d, nb_dofs %d, rank %d",
333 type, side, nb_dofs, RANK);
334 }
335 if (nb_dofs / RANK > data.getN().size2()) {
336 std::cerr << side << " " << type << " "
337 << ApproximationBaseNames[data.getBase()] << std::endl;
338 std::cerr << data.getN() << std::endl;
339 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
340 "data inconsistency nb_dofs >= data.N.size2(), i.e. %u >= %u",
341 nb_dofs, data.getN().size2());
342 }
343
344 if (type == MBVERTEX) {
345 dataAtGaussPts.resize(data.getN().size1(), RANK, false);
346 dataGradAtGaussPts.resize(data.getN().size1(), RANK * DIM, false);
347 dataAtGaussPts.clear();
348 dataGradAtGaussPts.clear();
349 }
350
351 CHKERR calculateValAndGrad(side, type, data);
352
354 }
355};
356
357/**
358 * \brief Specialization for field with 3 coefficients in 3 dimension
359 */
360template <>
361template <>
363OpGetDataAndGradient<3, 3>::getValAtGaussPtsTensor<3>(MatrixDouble &data);
364
365/**
366 * \brief Specialization for field with 3 coefficients in 3 dimension
367 */
368template <>
369template <>
371OpGetDataAndGradient<3, 3>::getGradAtGaussPtsTensor<3, 3>(MatrixDouble &data);
372
373/**
374 * \brief Specialization for field with 3 coefficients in 3 dimension
375 */
376template <>
379
380/**
381 * \brief Specialization for field with for scalar field in 3 dimension
382 */
383template <>
386
387/** \brief calculate normals at Gauss points of triangle element
388 * \ingroup mofem_forces_and_sources
389 */
391
400
402 MatrixDouble &coords_at_gaussptf3, MatrixDouble &normals_at_gaussptf3,
403 MatrixDouble &tangent1_at_gaussptf3, MatrixDouble &tangent2_at_gaussptf3,
404 MatrixDouble &coords_at_gaussptf4, MatrixDouble &normals_at_gaussptf4,
405 MatrixDouble &tangent1_at_gaussptf4, MatrixDouble &tangent2_at_gaussptf4)
406 : cOords_at_GaussPtF3(coords_at_gaussptf3),
407 nOrmals_at_GaussPtF3(normals_at_gaussptf3),
408 tAngent1_at_GaussPtF3(tangent1_at_gaussptf3),
409 tAngent2_at_GaussPtF3(tangent2_at_gaussptf3),
410 cOords_at_GaussPtF4(coords_at_gaussptf4),
411 nOrmals_at_GaussPtF4(normals_at_gaussptf4),
412 tAngent1_at_GaussPtF4(tangent1_at_gaussptf4),
413 tAngent2_at_GaussPtF4(tangent2_at_gaussptf4) {}
414
418
420};
421
422/** \brief transform Hdiv base fluxes from reference element to physical
423 * triangle \ingroup mofem_forces_and_sources
424 *
425 * \deprecated It is used in contact elements. Contact elements should be
426 * minified to work as face element,
427 */
429
430 const VectorDouble
431 *normalRawPtr; //< Normal of the element for linear geometry
432 const MatrixDouble *normalsAtGaussPtsRawPtr; //< Normals at integration points
433 // for nonlinear geometry
434
435 /**
436 * @brief Shift in vector for linear geometry
437 *
438 * Normal can have size larger than three, for example normal for contact
439 * prism and flat prims element have six comonents, for top and bottom
440 * triangle of the prims.
441 *
442 */
444
446 const MatrixDouble &normals_at_pts,
447 const int normal_shift = 0)
448 : normalRawPtr(&normal), normalsAtGaussPtsRawPtr(&normals_at_pts),
449 normalShift(normal_shift) {}
450
452 const VectorDouble *normal_raw_ptr = nullptr,
453 const MatrixDouble *normals_at_pts_ptr = nullptr,
454 const int normal_shift = 0)
455 : normalRawPtr(normal_raw_ptr),
456 normalsAtGaussPtsRawPtr(normals_at_pts_ptr), normalShift(normal_shift) {
457 }
458
461};
462
463/** \brief transform Hcurl base fluxes from reference element to physical
464 * triangle \ingroup mofem_forces_and_sources
465 *
466 * \deprecated It is used in contact elements. Contact elements should be
467 * minified to work as face element,
468 */
470
477
479 const MatrixDouble &normals_at_pts,
480 const VectorDouble &tangent0,
481 const MatrixDouble &tangent0_at_pts,
482 const VectorDouble &tangent1,
483 const MatrixDouble &tangent1_at_pts)
484 : nOrmal(normal), normalsAtGaussPts(normals_at_pts), tAngent0(tangent0),
485 tangent0AtGaussPt(tangent0_at_pts), tAngent1(tangent1),
486 tangent1AtGaussPt(tangent1_at_pts) {}
487
490};
491
492/** \brief Calculate tangent vector on edge form HO geometry approximation
493 * \ingroup mofem_forces_and_sources
494 */
496
498
500
503};
504
505/** \brief transform Hcurl base fluxes from reference element to physical edge
506 * \ingroup mofem_forces_and_sources
507 */
509
512
514 const MatrixDouble &tangent_at_pts)
515 : tAngent(tangent), tangentAtGaussPt(tangent_at_pts) {}
516
519};
520
521} // namespace MoFEM
522
523#endif //__DATAOPERATORS_HPP
524
525/**
526 * \defgroup mofem_forces_and_sources Forces and sources
527 **/
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:45
static const char *const ApproximationBaseNames[]
Definition: definitions.h:85
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:574
FTensor::Tensor1< double, SPACE_DIM > normal(FTensor::Tensor1< T1, 3 > &t_coords, FTensor::Tensor1< T2, SPACE_DIM > &t_disp)
Definition: ContactOps.hpp:77
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
base operator to do operations at Gauss Pt. level
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
virtual MoFEMErrorCode opLhs(EntitiesFieldData &row_data, EntitiesFieldData &col_data)
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
bool & doTris
\deprectaed
bool & doPrisms
\deprectaed
bool sYmm
If true assume that matrix is symmetric structure.
boost::function< MoFEMErrorCode(DataOperator *op_ptr, int side, EntityType type, EntitiesFieldData::EntData &data)> DoWorkRhsHookFunType
boost::function< MoFEMErrorCode(DataOperator *op_ptr, int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)> DoWorkLhsHookFunType
void setSymm()
set if operator is executed taking in account symmetry
virtual MoFEMErrorCode opRhs(EntitiesFieldData &data, const bool error_if_no_base=false)
bool getSymm() const
Get if operator uses symmetry of DOFs or not.
virtual MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
bool & doVertices
\deprectaed If false skip vertices
bool & doTets
\deprectaed
DataOperator(const bool symm=true)
virtual ~DataOperator()=default
void unSetSymm()
unset if operator is executed for non symmetric problem
DoWorkLhsHookFunType doWorkLhsHook
DoWorkRhsHookFunType doWorkRhsHook
bool & doEdges
\deprectaed If false skip edges
bool & doQuads
\deprectaed
Data on single entity (This is passed as argument to DataOperator::doWork)
FieldApproximationBase & getBase()
Get approximation base.
MatrixDouble & getDiffN(const FieldApproximationBase base)
get derivatives of base functions
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
data structure for finite element entity
calculate normals at Gauss points of triangle element
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpGetCoordsAndNormalsOnPrism(MatrixDouble &coords_at_gaussptf3, MatrixDouble &normals_at_gaussptf3, MatrixDouble &tangent1_at_gaussptf3, MatrixDouble &tangent2_at_gaussptf3, MatrixDouble &coords_at_gaussptf4, MatrixDouble &normals_at_gaussptf4, MatrixDouble &tangent1_at_gaussptf4, MatrixDouble &tangent2_at_gaussptf4)
Get field values and gradients at Gauss points.
FTensor::Tensor2< double *, R, D > getGradAtGaussPtsTensor(MatrixDouble &data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode calculateValAndGrad(int side, EntityType type, EntitiesFieldData::EntData &data)
Calculate gradient and values at integration points.
OpGetDataAndGradient(MatrixDouble &data_at_gauss_pt, MatrixDouble &data_grad_at_gauss_pt)
FTensor::Tensor1< double *, R > getValAtGaussPtsTensor(MatrixDouble &data)
Calculate tangent vector on edge form HO geometry approximation.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpGetHOTangentOnEdge(MatrixDouble &tangent)
apply contravariant (Piola) transfer to Hdiv space
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
FTensor::Tensor2< double *, 3, 3 > tJac
OpSetContravariantPiolaTransform(double &volume, MatrixDouble3by3 &jac)
transform Hdiv base fluxes from reference element to physical triangle
OpSetContravariantPiolaTransformOnFace(const VectorDouble &normal, const MatrixDouble &normals_at_pts, const int normal_shift=0)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSetContravariantPiolaTransformOnFace(const VectorDouble *normal_raw_ptr=nullptr, const MatrixDouble *normals_at_pts_ptr=nullptr, const int normal_shift=0)
int normalShift
Shift in vector for linear geometry.
apply covariant transfer to Hcurl space
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSetCovariantPiolaTransform(MatrixDouble3by3 &inv_jac)
FTensor::Tensor2< double *, 3, 3 > tInvJac
transform Hcurl base fluxes from reference element to physical edge
OpSetCovariantPiolaTransformOnEdge(const VectorDouble &tangent, const MatrixDouble &tangent_at_pts)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
transform Hcurl base fluxes from reference element to physical triangle
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSetCovariantPiolaTransformOnFace(const VectorDouble &normal, const MatrixDouble &normals_at_pts, const VectorDouble &tangent0, const MatrixDouble &tangent0_at_pts, const VectorDouble &tangent1, const MatrixDouble &tangent1_at_pts)
Transform local reference derivatives of shape function to global derivatives.
FTensor::Index< 'j', 3 > j
MatrixDouble diffNinvJac
FTensor::Index< 'i', 3 > i
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSetInvJacH1(MatrixDouble3by3 &inv_jac)
FTensor::Tensor2< double *, 3, 3 > tInvJac
brief Transform local reference derivatives of shape function to global derivatives
FTensor::Index< 'j', 3 > j
FTensor::Tensor2< double *, 3, 3 > tInvJac
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpSetInvJacHdivAndHcurl(MatrixDouble3by3 &inv_jac)
FTensor::Index< 'i', 3 > i
FTensor::Index< 'k', 3 > k