v0.14.0
Loading...
Searching...
No Matches
FormsIntegrators.hpp
Go to the documentation of this file.
1/** \file FormsIntegrators.hpp
2 * \brief Forms integrators
3 * \ingroup mofem_form
4
5*/
6
7#ifndef __FORMS_INTEGRATORS_HPP__
8#define __FORMS_INTEGRATORS_HPP__
9
10namespace MoFEM {
11
12//! [Storage and set boundary conditions]
13
17 /**
18 * @brief Store modifed indices by field
19 *
20 * Hash map, key is field name, value is storage.
21 *
22 */
24 map<std::string, std::vector<boost::shared_ptr<EssentialBcStorage>>>;
26};
27
28/**
29 * @brief Set indices on entities on finite element
30 * @ingroup mofem_forms
31 *
32 * If index is marked, its value is set to -1. DOF with such index is not
33 * assembled into the system.
34 *
35 * Indices are stored on the entity.
36 *
37 */
39 OpSetBc(std::string field_name, bool yes_set,
40 boost::shared_ptr<std::vector<unsigned char>> boundary_marker);
41 MoFEMErrorCode doWork(int side, EntityType type,
43
44public:
45 bool yesSet;
46 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
47};
48
50 OpUnSetBc(std::string field_name);
51 MoFEMErrorCode doWork(int side, EntityType type,
53};
54
55/**
56 * @brief Set values to vector in operator
57 * @ingroup mofem_forms
58 *
59 * MoFEM::FieldEntity provides MoFEM::FieldEntity::getWeakStoragePtr() storage
60 * function which allows to transfer data between FEs or operators processing
61 * the same entities.
62 *
63 * When MoFEM::OpSetBc is pushed in weak storage indices taking in account
64 * indices which are skip to take boundary conditions are stored. Those entities
65 * are used by VecSetValues.
66 *
67 * @param V
68 * @param data
69 * @param ptr
70 * @param iora
71 * @return MoFEMErrorCode
72 */
73template <>
76 const double *ptr, InsertMode iora);
77
78/**
79 * @brief Set values to matrix in operator
80 *
81 * See MoFEM::VecSetValues<EssentialBcStorage> for explanation.
82 *
83 * @param M
84 * @param row_data
85 * @param col_data
86 * @param ptr
87 * @param iora
88 * @return MoFEMErrorCode
89 */
90template <>
93 const EntitiesFieldData::EntData &row_data,
94 const EntitiesFieldData::EntData &col_data,
95 const double *ptr, InsertMode iora);
96
97//! [Storage and set boundary conditions]
98
99/**
100 * @brief Form integrator assembly types
101 * @ingroup mofem_forms
102 *
103 */
105
106template <int A> struct AssemblyTypeSelector {};
107
108template <>
109inline MoFEMErrorCode MatSetValues<AssemblyTypeSelector<PETSC>>(
110 Mat M, const EntitiesFieldData::EntData &row_data,
111 const EntitiesFieldData::EntData &col_data, const double *ptr,
112 InsertMode iora) {
113 return MatSetValues<EssentialBcStorage>(M, row_data, col_data, ptr, iora);
114}
115
116template <>
117inline MoFEMErrorCode VecSetValues<AssemblyTypeSelector<PETSC>>(
118 Vec V, const EntitiesFieldData::EntData &data, const double *ptr,
119 InsertMode iora) {
120 return VecSetValues<EssentialBcStorage>(V, data, ptr, iora);
121}
122
123/**
124 * @brief Form integrator integration types
125 * @ingroup mofem_forms
126 *
127 */
129
130/**
131 * @brief Scalar function type
132 * @ingroup mofem_forms
133 *
134 */
136 boost::function<double(const double, const double, const double)>;
137
138inline double scalar_fun_one(const double, const double, const double) {
139 return 1;
140}
141
142/**
143 * @brief Lambda function used to scale with time
144 *
145 */
146using TimeFun = boost::function<double(double)>;
147
148/**
149 * @brief Lambda function used to scale with time
150 *
151 */
152using FEFun = boost::function<double(const FEMethod *fe_ptr)>;
153
154/**
155 * @brief Constant function type
156 *
157 */
158using ConstantFun = boost::function<double()>;
159
160/**
161 * @brief Vector function type
162 * @ingroup mofem_forms
163 *
164 * @tparam DIM dimension of the return
165 */
166template <int DIM>
167using VectorFun = boost::function<FTensor::Tensor1<double, DIM>(
168 const double, const double, const double)>;
169
170template <AssemblyType A, typename EleOp> struct OpBaseImpl : public EleOp {
171 using OpType = typename EleOp::OpType;
173
174 OpBaseImpl(const std::string row_field_name, const std::string col_field_name,
175 const OpType type, boost::shared_ptr<Range> ents_ptr = nullptr)
176 : EleOp(row_field_name, col_field_name, type, false),
177 assembleTranspose(false), onlyTranspose(false), entsPtr(ents_ptr) {}
178
179 /**
180 * \brief Do calculations for the left hand side
181 * @param row_side row side number (local number) of entity on element
182 * @param col_side column side number (local number) of entity on element
183 * @param row_type type of row entity
184 * @param col_type type of column entity
185 * @param row_data data for row
186 * @param col_data data for column
187 * @return error code
188 */
189 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
190 EntityType col_type, EntData &row_data,
191 EntData &col_data);
192
193 /**
194 * @brief Do calculations for the right hand side
195 *
196 * @param row_side
197 * @param row_type
198 * @param row_data
199 * @return MoFEMErrorCode
200 */
201 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
202
203 TimeFun timeScalingFun; ///< assumes that time variable is set
204 FEFun feScalingFun; ///< assumes that time variable is set
205 boost::shared_ptr<Range> entsPtr; ///< Entities on which element is run
206
207 using MatSetValuesHook = boost::function<MoFEMErrorCode(
209 const EntitiesFieldData::EntData &row_data,
210 const EntitiesFieldData::EntData &col_data, MatrixDouble &m)>;
211
213
214protected:
215 template <int DIM>
217 return getFTensor1FromArray<DIM, DIM>(locF);
218 }
219
220 template <int DIM>
222 getLocMat(const int rr) {
223 return getFTensor2FromArray<DIM, DIM, DIM>(locMat, rr);
224 }
225
226 int nbRows; ///< number of dofs on rows
227 int nbCols; ///< number if dof on column
228 int nbIntegrationPts; ///< number of integration points
229 int nbRowBaseFunctions; ///< number or row base functions
230
231 int rowSide; ///< row side number
232 int colSide; ///< column side number
233 EntityType rowType; ///< row type
234 EntityType colType; ///< column type
235
238
239 MatrixDouble locMat; ///< local entity block matrix
240 MatrixDouble locMatTranspose; ///< local entity block matrix
241 VectorDouble locF; ///< local entity vector
242
243 /**
244 * \brief Integrate grad-grad operator
245 * @param row_data row data (consist base functions on row entity)
246 * @param col_data column data (consist base functions on column entity)
247 * @return error code
248 */
249 virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data) {
251 }
252
253 virtual MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data,
254 const bool trans);
255
256 /**
257 * \brief Class dedicated to integrate operator
258 * @param data entity data on element row
259 * @return error code
260 */
263 }
264
266
267 /** \brief Get number of base functions
268 *
269 * @param data
270 * @return number of base functions
271 */
273};
274
275template <AssemblyType A, typename EleOp>
279 const EntitiesFieldData::EntData &row_data,
280 const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
281 return MatSetValues<AssemblyTypeSelector<A>>(
282 op_ptr->getKSPB(), row_data, col_data, m, ADD_VALUES);
283 };
284
285/**
286 * @brief Integrator forms
287 * @ingroup mofem_forms
288 *
289 * @tparam EleOp
290 */
291template <typename EleOp> struct FormsIntegrators {
292
294 using OpType = typename EleOp::OpType;
295
296 /**
297 * @brief Assembly methods
298 * @ingroup mofem_forms
299 *
300 * @tparam A
301 */
302 template <AssemblyType A> struct Assembly {
303
305
306 /**
307 * @brief Linear form
308 * @ingroup mofem_forms
309 *
310 * @tparam I
311 */
312 template <IntegrationType I> struct LinearForm;
313
314 /**
315 * @brief Bi linear form
316 * @ingroup mofem_forms
317 *
318 * @tparam I
319 */
320 template <IntegrationType I> struct BiLinearForm;
321
322 /**
323 * @brief Tri linear form
324 * @ingroup mofem_forms
325 *
326 * @tparam I
327 */
328 template <IntegrationType I> struct TriLinearForm;
329
330 }; // Assembly
331}; // namespace MoFEM
332
333template <AssemblyType A, typename EleOp>
334size_t
336 auto nb_base_functions = data.getN().size2();
337 if (data.getBase() != USER_BASE) {
338 switch (data.getSpace()) {
339 case NOSPACE:
340 break;
341 case NOFIELD:
342 break;
343 case H1:
344 break;
345 case HCURL:
346 case HDIV:
347 nb_base_functions /= 3;
348#ifndef NDEBUG
349 if (data.getN().size2() % 3) {
350 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
351 "Number of base functions is not divisible by 3");
352 }
353#endif
354 break;
355 case L2:
356 break;
357 default:
358 SETERRQ1(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
359 "Space %s not implemented", FieldSpaceNames[data.getSpace()]);
360 }
361 }
362 return nb_base_functions;
363}
364
365template <AssemblyType A, typename EleOp>
367OpBaseImpl<A, EleOp>::doWork(int row_side, int col_side, EntityType row_type,
368 EntityType col_type,
370 EntitiesFieldData::EntData &col_data) {
372
373 if (entsPtr) {
374 if (entsPtr->find(this->getFEEntityHandle()) == entsPtr->end())
376 }
377
378 // get number of dofs on row
379 nbRows = row_data.getIndices().size();
380 // if no dofs on row, exit that work, nothing to do here
381 if (!nbRows)
383 rowSide = row_side;
384 rowType = row_type;
385 // get number of dofs on column
386 nbCols = col_data.getIndices().size();
387 // if no dofs on column, exit nothing to do here
388 if (!nbCols)
390 colSide = col_side;
391 colType = col_type;
392 // get number of integration points
393 nbIntegrationPts = EleOp::getGaussPts().size2();
394 // get row base functions
395 nbRowBaseFunctions = getNbOfBaseFunctions(row_data);
396
397 // set size of local entity bock
398 locMat.resize(nbRows, nbCols, false);
399 // clear matrix
400 locMat.clear();
401 // integrate local matrix for entity block
402 CHKERR this->iNtegrate(row_data, col_data);
403
404 // assemble local matrix
405 auto check_if_assemble_transpose = [&] {
406 if (this->sYmm) {
407 if (row_side != col_side || row_type != col_type)
408 return true;
409 else
410 return false;
411 } else if (assembleTranspose) {
412 return true;
413 }
414
415 return false;
416 };
417 CHKERR aSsemble(row_data, col_data, check_if_assemble_transpose());
419}
420
421template <AssemblyType A, typename EleOp>
423 EntData &row_data) {
425
426 if (entsPtr) {
427 if (entsPtr->find(this->getFEEntityHandle()) == entsPtr->end())
429 }
430
431 // get number of dofs on row
432 nbRows = row_data.getIndices().size();
433 rowSide = row_side;
434 rowType = row_type;
435
436 if (!nbRows)
438 // get number of integration points
439 nbIntegrationPts = EleOp::getGaussPts().size2();
440 // get row base functions
441 nbRowBaseFunctions = getNbOfBaseFunctions(row_data);
442 // resize and clear the right hand side vector
443 locF.resize(nbRows);
444 locF.clear();
445 // integrate local vector
446 CHKERR this->iNtegrate(row_data);
447 // assemble local vector
448 CHKERR this->aSsemble(row_data);
450}
451
452template <AssemblyType A, typename EleOp>
454 EntData &col_data,
455 const bool transpose) {
457
458 if (!this->timeScalingFun.empty())
459 this->locMat *= this->timeScalingFun(this->getFEMethod()->ts_t);
460 if (!this->feScalingFun.empty())
461 this->locMat *= this->feScalingFun(this->getFEMethod());
462
463 // Assemble transpose
464 if (transpose) {
465 this->locMatTranspose.resize(this->locMat.size2(), this->locMat.size1(),
466 false);
467 noalias(this->locMatTranspose) = trans(this->locMat);
468 CHKERR matSetValuesHook(this, col_data, row_data, this->locMatTranspose);
469 }
470
471 if (!this->onlyTranspose) {
472 // assemble local matrix
473 CHKERR matSetValuesHook(this, row_data, col_data, this->locMat);
474 }
475
477}
478
479template <AssemblyType A, typename EleOp>
481 if (!this->timeScalingFun.empty())
482 this->locF *= this->timeScalingFun(this->getFEMethod()->ts_t);
483 if (!this->feScalingFun.empty())
484 this->locF *= this->feScalingFun(this->getFEMethod());
485
486 return VecSetValues<AssemblyTypeSelector<A>>(this->getKSPf(), data,
487 this->locF, ADD_VALUES);
488}
489
490} // namespace MoFEM
491
492/**
493 * \defgroup mofem_forms Forms Integrators
494 *
495 * \brief Classes and functions used to evaluate fields at integration pts,
496 *jacobians, etc..
497 *
498 * \ingroup mofem_forces_and_sources
499 **/
500
501#endif //__FORMS_INTEGRATORS_HPP__
static Index< 'M', 3 > M
@ USER_BASE
user implemented approximation base
Definition: definitions.h:68
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ NOFIELD
scalar or vector of scalars describe (no true field)
Definition: definitions.h:84
@ H1
continuous field
Definition: definitions.h:85
@ NOSPACE
Definition: definitions.h:83
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
static const char *const FieldSpaceNames[]
Definition: definitions.h:92
#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
#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
FTensor::Index< 'm', SPACE_DIM > m
IntegrationType
Form integrator integration types.
MoFEMErrorCode VecSetValues< EssentialBcStorage >(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Set values to vector in operator.
boost::function< double(const double, const double, const double)> ScalarFun
Scalar function type.
AssemblyType
[Storage and set boundary conditions]
boost::function< FTensor::Tensor1< double, DIM >(const double, const double, const double)> VectorFun
Vector function type.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
UBlasVector< int > VectorInt
Definition: Types.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
boost::function< double(const FEMethod *fe_ptr)> FEFun
Lambda function used to scale with time.
boost::function< double(double)> TimeFun
Lambda function used to scale with time.
double scalar_fun_one(const double, const double, const double)
MoFEMErrorCode MatSetValues< EssentialBcStorage >(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Set values to matrix in operator.
boost::function< double()> ConstantFun
Constant function type.
constexpr auto field_name
Data on single entity (This is passed as argument to DataOperator::doWork)
FieldApproximationBase & getBase()
Get approximation base.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FieldSpace & getSpace()
Get field space.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
[Storage and set boundary conditions]
static HashVectorStorage feStorage
[Storage and set boundary conditions]
map< std::string, std::vector< boost::shared_ptr< EssentialBcStorage > > > HashVectorStorage
Store modifed indices by field.
EssentialBcStorage(VectorInt &indices)
structure for User Loop Methods on finite elements
typename EleOp::OpType OpType
OpBaseImpl(const std::string row_field_name, const std::string col_field_name, const OpType type, boost::shared_ptr< Range > ents_ptr=nullptr)
VectorDouble locF
local entity vector
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Do calculations for the right hand side.
typename EleOp::OpType OpType
TimeFun timeScalingFun
assumes that time variable is set
MatrixDouble locMatTranspose
local entity block matrix
int rowSide
row side number
int colSide
column side number
boost::shared_ptr< Range > entsPtr
Entities on which element is run.
static MatSetValuesHook matSetValuesHook
int nbRows
number of dofs on rows
EntityType colType
column type
int nbIntegrationPts
number of integration points
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Do calculations for the left hand side.
virtual MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data, const bool trans)
virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate grad-grad operator.
virtual MoFEMErrorCode iNtegrate(EntData &data)
Class dedicated to integrate operator.
MatrixDouble locMat
local entity block matrix
virtual size_t getNbOfBaseFunctions(EntitiesFieldData::EntData &data)
Get number of base functions.
FTensor::Tensor2< FTensor::PackPtr< double *, DIM >, DIM, DIM > getLocMat(const int rr)
FEFun feScalingFun
assumes that time variable is set
FTensor::Tensor1< FTensor::PackPtr< double *, DIM >, DIM > getNf()
virtual MoFEMErrorCode aSsemble(EntData &data)
int nbCols
number if dof on column
int nbRowBaseFunctions
number or row base functions
EntityType rowType
row type
boost::function< MoFEMErrorCode(ForcesAndSourcesCore::UserDataOperator *op_ptr, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, MatrixDouble &m)> MatSetValuesHook
Set indices on entities on finite element.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.