v0.13.1
Loading...
Searching...
No Matches
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
13
14#ifndef __DATAOPERATORS_HPP
15#define __DATAOPERATORS_HPP
16
17using namespace boost::numeric;
18
19namespace MoFEM {
20
21/** \brief base operator to do operations at Gauss Pt. level
22 * \ingroup mofem_forces_and_sources
23 */
25
26 DataOperator(const bool symm = true);
27
28 virtual ~DataOperator() = default;
29
30 using DoWorkLhsHookFunType = boost::function<MoFEMErrorCode(
31 DataOperator *op_ptr, int row_side, int col_side, EntityType row_type,
32 EntityType col_type, EntitiesFieldData::EntData &row_data,
34
36
37 /** \brief Operator for bi-linear form, usually to calculate values on
38 * left hand side
39 */
40 virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
41 EntityType col_type,
45 if (doWorkLhsHook) {
46 CHKERR doWorkLhsHook(this, row_side, col_side, row_type, col_type,
47 row_data, col_data);
48 } else {
49 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
50 "doWork function is not implemented for this operator");
51 }
53 }
54
55 virtual MoFEMErrorCode opLhs(EntitiesFieldData &row_data,
56 EntitiesFieldData &col_data);
57
58 using DoWorkRhsHookFunType = boost::function<MoFEMErrorCode(
59 DataOperator *op_ptr, int side, EntityType type,
61
63
64 /** \brief Operator for linear form, usually to calculate values on right hand
65 * side
66 */
67 virtual MoFEMErrorCode doWork(int side, EntityType type,
70 if (doWorkRhsHook) {
71 CHKERR doWorkRhsHook(this, side, type, data);
72 } else {
73 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
74 "doWork function is not implemented for this operator");
75 }
77 }
78
80 const bool error_if_no_base = false);
81
82 bool sYmm; ///< If true assume that matrix is symmetric structure
83
84 std::array<bool, MBMAXTYPE>
85 doEntities; ///< If true operator is executed for entity.
86
87 // Deprecated variables. Use doEntities instead. I keep them for back
88 // compatibility with some older modules. It will be removed in some future.
89
90 bool &doVertices; ///< \deprectaed If false skip vertices
91 bool &doEdges; ///< \deprectaed If false skip edges
92 bool &doQuads; ///< \deprectaed
93 bool &doTris; ///< \deprectaed
94 bool &doTets; ///< \deprectaed
95 bool &doPrisms; ///< \deprectaed
96
97 /**
98 * \brief Get if operator uses symmetry of DOFs or not
99 *
100 * If symmetry is used, only not repeating combinations of entities are
101 * looped. For an example pair of (Vertex, Edge_0) and (Edge_0, Vertex) will
102 * calculate the same matrices only transposed. Implementing that this can be
103 * exploited by integrating only one pair.
104 *
105 * @return true if symmetry
106 */
107 inline bool getSymm() const { return sYmm; }
108
109 /// set if operator is executed taking in account symmetry
110 inline void setSymm() { sYmm = true; }
111
112 /// unset if operator is executed for non symmetric problem
113 inline void unSetSymm() { sYmm = false; }
114
115private:
116 template <bool Symm>
117 inline MoFEMErrorCode opLhs(EntitiesFieldData &row_data,
118 EntitiesFieldData &col_data);
119
120 template <bool ErrorIfNoBase>
122 const std::array<bool, MBMAXTYPE> &do_entities);
123};
124
125/**
126 * \brief Transform local reference derivatives of shape function to global
127 derivatives
128
129 * \ingroup mofem_forces_and_sources
130 */
132
136
138 : tInvJac(&inv_jac(0, 0), &inv_jac(0, 1), &inv_jac(0, 2), &inv_jac(1, 0),
139 &inv_jac(1, 1), &inv_jac(1, 2), &inv_jac(2, 0), &inv_jac(2, 1),
140 &inv_jac(2, 2)) {}
141
143 MoFEMErrorCode doWork(int side, EntityType type,
145};
146
147/**
148 * \brief brief Transform local reference derivatives of shape function to
149 global derivatives
150
151 * \ingroup mofem_forces_and_sources
152 */
154
159
161 : tInvJac(&inv_jac(0, 0), &inv_jac(0, 1), &inv_jac(0, 2), &inv_jac(1, 0),
162 &inv_jac(1, 1), &inv_jac(1, 2), &inv_jac(2, 0), &inv_jac(2, 1),
163 &inv_jac(2, 2)) {}
164
166 MoFEMErrorCode doWork(int side, EntityType type,
168};
169
170/** \brief apply contravariant (Piola) transfer to Hdiv space
171
172Contravariant Piola transformation
173\f[
174\psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
175\left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
176=
177\frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
178\f]
179
180* \ingroup mofem_forces_and_sources
181
182*/
184
185 double &vOlume;
186
191
193 : vOlume(volume),
194 // jAc(jac),
195 tJac(&jac(0, 0), &jac(0, 1), &jac(0, 2), &jac(1, 0), &jac(1, 1),
196 &jac(1, 2), &jac(2, 0), &jac(2, 1), &jac(2, 2)) {}
197
200
201 MoFEMErrorCode doWork(int side, EntityType type,
203};
204
205/** \brief apply covariant transfer to Hcurl space
206
207Contravariant Piola transformation
208\f[
209\psi_i|_t = \frac{1}{\textrm{det}(J)}J_{ij}\hat{\psi}_j\\
210\left.\frac{\partial \psi_i}{\partial \xi_j}\right|_t
211=
212\frac{1}{\textrm{det}(J)}J_{ik}\frac{\partial \hat{\psi}_k}{\partial \xi_j}
213\f]
214
215
216* \ingroup mofem_forces_and_sources
217*/
219
224
226 : tInvJac(&inv_jac(0, 0), &inv_jac(0, 1), &inv_jac(0, 2), &inv_jac(1, 0),
227 &inv_jac(1, 1), &inv_jac(1, 2), &inv_jac(2, 0), &inv_jac(2, 1),
228 &inv_jac(2, 2)) {}
229
232
233 MoFEMErrorCode doWork(int side, EntityType type,
235};
236
237/** \brief Get field values and gradients at Gauss points
238 * \ingroup mofem_forces_and_sources
239 */
240template <int RANK, int DIM> struct OpGetDataAndGradient : public DataOperator {
241
244
246 MatrixDouble &data_grad_at_gauss_pt)
247 : dataAtGaussPts(data_at_gauss_pt),
248 dataGradAtGaussPts(data_grad_at_gauss_pt) {}
249
250 /**
251 * Return tensor associated with matrix storing values
252 */
253 template <int R>
255 THROW_MESSAGE("Not implemented");
256 }
257
258 /**
259 * Return tensor associated with matrix storing gradient values
260 */
261 template <int R, int D>
263 THROW_MESSAGE("Not implemented");
264 }
265
266 /**
267 * \brief Calculate gradient and values at integration points
268 * @param side side of entity on element
269 * @param type type of entity
270 * @param data data stored on entity (dofs values, dofs indices, etc.)
271 * @return error code
272 */
276 const int nb_base_functions = data.getN().size2();
277 bool constant_diff = false;
278 if (type == MBVERTEX && data.getDiffN().size1() * data.getDiffN().size2() ==
279 DIM * nb_base_functions) {
280 constant_diff = true;
281 }
282 const int nb_dofs = data.getFieldData().size();
283 for (unsigned int gg = 0; gg < data.getN().size1(); gg++) {
284 double *data_ptr, *n_ptr, *diff_n_ptr;
285 n_ptr = &data.getN()(gg, 0);
286 if (constant_diff) {
287 diff_n_ptr = &data.getDiffN()(0, 0);
288 } else {
289 diff_n_ptr = &data.getDiffN()(gg, 0);
290 }
291 data_ptr = &*data.getFieldData().data().begin();
292 for (int rr = 0; rr < RANK; rr++, data_ptr++) {
293 dataAtGaussPts(gg, rr) +=
294 cblas_ddot(nb_dofs / RANK, n_ptr, 1, data_ptr, RANK);
295 double *diff_n_ptr2 = diff_n_ptr;
296 for (unsigned int dd = 0; dd < DIM; dd++, diff_n_ptr2++) {
297 dataGradAtGaussPts(gg, DIM * rr + dd) +=
298 cblas_ddot(nb_dofs / RANK, diff_n_ptr2, DIM, data_ptr, RANK);
299 }
300 }
301 }
303 }
304
307
309
310 if (data.getFieldData().size() == 0) {
312 }
313
314 unsigned int nb_dofs = data.getFieldData().size();
315 if (nb_dofs == 0)
317
318 if (nb_dofs % RANK != 0) {
319 SETERRQ4(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
320 "data inconsistency, type %d, side %d, nb_dofs %d, rank %d",
321 type, side, nb_dofs, RANK);
322 }
323 if (nb_dofs / RANK > data.getN().size2()) {
324 std::cerr << side << " " << type << " "
325 << ApproximationBaseNames[data.getBase()] << std::endl;
326 std::cerr << data.getN() << std::endl;
327 SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
328 "data inconsistency nb_dofs >= data.N.size2(), i.e. %u >= %u",
329 nb_dofs, data.getN().size2());
330 }
331
332 if (type == MBVERTEX) {
333 dataAtGaussPts.resize(data.getN().size1(), RANK, false);
334 dataGradAtGaussPts.resize(data.getN().size1(), RANK * DIM, false);
335 dataAtGaussPts.clear();
336 dataGradAtGaussPts.clear();
337 }
338
339 CHKERR calculateValAndGrad(side, type, data);
340
342 }
343};
344
345/**
346 * \brief Specialization for field with 3 coefficients in 3 dimension
347 */
348template <>
349template <>
351OpGetDataAndGradient<3, 3>::getValAtGaussPtsTensor<3>(MatrixDouble &data);
352
353/**
354 * \brief Specialization for field with 3 coefficients in 3 dimension
355 */
356template <>
357template <>
359OpGetDataAndGradient<3, 3>::getGradAtGaussPtsTensor<3, 3>(MatrixDouble &data);
360
361/**
362 * \brief Specialization for field with 3 coefficients in 3 dimension
363 */
364template <>
366 int side, EntityType type, EntitiesFieldData::EntData &data);
367
368/**
369 * \brief Specialization for field with for scalar field in 3 dimension
370 */
371template <>
373 int side, EntityType type, EntitiesFieldData::EntData &data);
374
375/** \brief calculate normals at Gauss points of triangle element
376 * \ingroup mofem_forces_and_sources
377 */
379
388
390 MatrixDouble &coords_at_gaussptf3, MatrixDouble &normals_at_gaussptf3,
391 MatrixDouble &tangent1_at_gaussptf3, MatrixDouble &tangent2_at_gaussptf3,
392 MatrixDouble &coords_at_gaussptf4, MatrixDouble &normals_at_gaussptf4,
393 MatrixDouble &tangent1_at_gaussptf4, MatrixDouble &tangent2_at_gaussptf4)
394 : cOords_at_GaussPtF3(coords_at_gaussptf3),
395 nOrmals_at_GaussPtF3(normals_at_gaussptf3),
396 tAngent1_at_GaussPtF3(tangent1_at_gaussptf3),
397 tAngent2_at_GaussPtF3(tangent2_at_gaussptf3),
398 cOords_at_GaussPtF4(coords_at_gaussptf4),
399 nOrmals_at_GaussPtF4(normals_at_gaussptf4),
400 tAngent1_at_GaussPtF4(tangent1_at_gaussptf4),
401 tAngent2_at_GaussPtF4(tangent2_at_gaussptf4) {}
402
404 MoFEMErrorCode doWork(int side, EntityType type,
406
408};
409
410/** \brief transform Hdiv base fluxes from reference element to physical
411 * triangle \ingroup mofem_forces_and_sources
412 *
413 * \deprecated It is used in contact elements. Contact elements should be
414 * minified to work as face element,
415 */
417
418 const VectorDouble
419 *normalRawPtr; //< Normal of the element for linear geometry
420 const MatrixDouble *normalsAtGaussPtsRawPtr; //< Normals at integration points
421 // for nonlinear geometry
422
423 /**
424 * @brief Shift in vector for linear geometry
425 *
426 * Normal can have size larger than three, for example normal for contact
427 * prism and flat prims element have six comonents, for top and bottom
428 * triangle of the prims.
429 *
430 */
432
434 const MatrixDouble &normals_at_pts,
435 const int normal_shift = 0)
436 : normalRawPtr(&normal), normalsAtGaussPtsRawPtr(&normals_at_pts),
437 normalShift(normal_shift) {}
438
440 const VectorDouble *normal_raw_ptr = nullptr,
441 const MatrixDouble *normals_at_pts_ptr = nullptr,
442 const int normal_shift = 0)
443 : normalRawPtr(normal_raw_ptr),
444 normalsAtGaussPtsRawPtr(normals_at_pts_ptr), normalShift(normal_shift) {
445 }
446
447 MoFEMErrorCode doWork(int side, EntityType type,
449};
450
451/** \brief transform Hcurl base fluxes from reference element to physical
452 * triangle \ingroup mofem_forces_and_sources
453 *
454 * \deprecated It is used in contact elements. Contact elements should be
455 * minified to work as face element,
456 */
458
465
467 const MatrixDouble &normals_at_pts,
468 const VectorDouble &tangent0,
469 const MatrixDouble &tangent0_at_pts,
470 const VectorDouble &tangent1,
471 const MatrixDouble &tangent1_at_pts)
472 : nOrmal(normal), normalsAtGaussPts(normals_at_pts), tAngent0(tangent0),
473 tangent0AtGaussPt(tangent0_at_pts), tAngent1(tangent1),
474 tangent1AtGaussPt(tangent1_at_pts) {}
475
476 MoFEMErrorCode doWork(int side, EntityType type,
478};
479
480/** \brief Calculate tangent vector on edge form HO geometry approximation
481 * \ingroup mofem_forces_and_sources
482 */
484
486
488
489 MoFEMErrorCode doWork(int side, EntityType type,
491};
492
493/** \brief transform Hcurl base fluxes from reference element to physical edge
494 * \ingroup mofem_forces_and_sources
495 */
497
500
502 const MatrixDouble &tangent_at_pts)
503 : tAngent(tangent), tangentAtGaussPt(tangent_at_pts) {}
504
505 MoFEMErrorCode doWork(int side, EntityType type,
507};
508
509} // namespace MoFEM
510
511#endif //__DATAOPERATORS_HPP
512
513/**
514 * \defgroup mofem_forces_and_sources Forces and sources
515 **/
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
static const char *const ApproximationBaseNames[]
Definition: definitions.h:72
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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
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