v0.15.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 */
113
114template <int A> struct AssemblyTypeSelector {};
115
116template <>
118 Mat M, const EntitiesFieldData::EntData &row_data,
119 const EntitiesFieldData::EntData &col_data, const double *ptr,
120 InsertMode iora) {
121 return MatSetValues<EssentialBcStorage>(M, row_data, col_data, ptr, iora);
122}
123
124template <>
126 Vec V, const EntitiesFieldData::EntData &data, const double *ptr,
127 InsertMode iora) {
128 return VecSetValues<EssentialBcStorage>(V, data, ptr, iora);
129}
130
131/**
132 * @brief Form integrator integration types
133 * @ingroup mofem_forms
134 *
135 */
137
138/**
139 * @brief Scalar function type
140 * @ingroup mofem_forms
141 *
142 */
144 boost::function<double(const double, const double, const double)>;
145
146inline double scalar_fun_one(const double, const double, const double) {
147 return 1;
148}
149
150/**
151 * @brief Lambda function used to scale with time
152 *
153 */
154using TimeFun = boost::function<double(double)>;
155
156/**
157 * @brief Lambda function used to scale with time
158 *
159 */
160using FEFun = boost::function<double(const FEMethod *fe_ptr)>;
161
162/**
163 * @brief Constant function type
164 *
165 */
166using ConstantFun = boost::function<double()>;
167
168/**
169 * @brief Vector function type
170 * @ingroup mofem_forms
171 *
172 * @tparam DIM dimension of the return
173 */
174template <int DIM>
175using VectorFun = boost::function<FTensor::Tensor1<double, DIM>(
176 const double, const double, const double)>;
177
178template <AssemblyType A, typename EleOp> struct OpBaseImpl : public EleOp {
179 using OpType = typename EleOp::OpType;
181
182 OpBaseImpl(const std::string row_field_name, const std::string col_field_name,
183 const OpType type, boost::shared_ptr<Range> ents_ptr = nullptr)
184 : EleOp(row_field_name, col_field_name, type, false),
185 assembleTranspose(false), onlyTranspose(false), entsPtr(ents_ptr) {}
186
187 /**
188 * \brief Do calculations for the left hand side
189 * @param row_side row side number (local number) of entity on element
190 * @param col_side column side number (local number) of entity on element
191 * @param row_type type of row entity
192 * @param col_type type of column entity
193 * @param row_data data for row
194 * @param col_data data for column
195 * @return error code
196 */
197 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
198 EntityType col_type, EntData &row_data,
199 EntData &col_data);
200
201 /**
202 * @brief Do calculations for the right hand side
203 *
204 * @param row_side
205 * @param row_type
206 * @param row_data
207 * @return MoFEMErrorCode
208 */
209 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
210
211 TimeFun timeScalingFun; ///< assumes that time variable is set
212 FEFun feScalingFun; ///< assumes that time variable is set
213 boost::shared_ptr<Range> entsPtr; ///< Entities on which element is run
214
215 using MatSetValuesHook = boost::function<MoFEMErrorCode(
217 const EntitiesFieldData::EntData &row_data,
218 const EntitiesFieldData::EntData &col_data, MatrixDouble &m)>;
219
221
222protected:
223 using EleOp::EleOp;
224
225 template <int DIM>
229
230 template <int DIM>
232 getLocMat(const int rr) {
234 }
235
236 int nbRows; ///< number of dofs on rows
237 int nbCols; ///< number if dof on column
238 int nbIntegrationPts; ///< number of integration points
239 int nbRowBaseFunctions; ///< number or row base functions
240
241 int rowSide; ///< row side number
242 int colSide; ///< column side number
243 EntityType rowType; ///< row type
244 EntityType colType; ///< column type
245
248
249 MatrixDouble locMat; ///< local entity block matrix
250 MatrixDouble locMatTranspose; ///< local entity block matrix
251 VectorDouble locF; ///< local entity vector
252
253 /**
254 * \brief Integrate grad-grad operator
255 * @param row_data row data (consist base functions on row entity)
256 * @param col_data column data (consist base functions on column entity)
257 * @return error code
258 */
259 virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data) {
261 }
262
263 virtual MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data,
264 const bool trans);
265
266 /**
267 * \brief Class dedicated to integrate operator
268 * @param data entity data on element row
269 * @return error code
270 */
273 }
274
275 virtual MoFEMErrorCode aSsemble(EntData &data);
276
277 /** \brief Get number of base functions
278 *
279 * @param data
280 * @return number of base functions
281 */
283};
284
285template <AssemblyType A, typename EleOp>
288 [](ForcesAndSourcesCore::UserDataOperator *op_ptr,
289 const EntitiesFieldData::EntData &row_data,
290 const EntitiesFieldData::EntData &col_data, MatrixDouble &m) {
292 op_ptr->getKSPB(), row_data, col_data, m, ADD_VALUES);
293 };
294
295template <typename OpBase>
296struct OpBrokenBaseImpl; //< This is base operator for broken spaces.
297 //Implementation is in
298 //FormsBrokenSpaceConstraintImpl.hpp
299
300/**
301 * @brief Integrator forms
302 * @ingroup mofem_forms
303 *
304 * @tparam EleOp
305 */
306template <typename EleOp> struct FormsIntegrators {
307
309 using OpType = typename EleOp::OpType;
310
311 /**
312 * @brief Assembly methods
313 * @ingroup mofem_forms
314 *
315 * @tparam A
316 */
317 template <AssemblyType A> struct Assembly {
318
321
322 /**
323 * @brief Linear form
324 * @ingroup mofem_forms
325 *
326 * @tparam I
327 */
328 template <IntegrationType I> struct LinearForm;
329
330 /**
331 * @brief Bi linear form
332 * @ingroup mofem_forms
333 *
334 * @tparam I
335 */
336 template <IntegrationType I> struct BiLinearForm;
337
338 /**
339 * @brief Tri linear form
340 * @ingroup mofem_forms
341 *
342 * @tparam I
343 */
344 template <IntegrationType I> struct TriLinearForm;
345
346 }; // Assembly
347}; // namespace MoFEM
348
349template <AssemblyType A, typename EleOp>
350size_t
352 if (!data.getNSharedPtr(data.getBase())) {
353 #ifndef NDEBUG
354 MOFEM_LOG("SELF", Sev::warning)
355 << "Ptr to base " << ApproximationBaseNames[data.getBase()]
356 << " functions is null, returning 0";
357 #endif // NDEBUG
358 return 0;
359 }
360 auto nb_base_functions = data.getN().size2();
361 if (data.getBase() != USER_BASE) {
362 switch (data.getSpace()) {
363 case NOSPACE:
364 break;
365 case NOFIELD:
366 break;
367 case H1:
368 break;
369 case HCURL:
370 case HDIV:
371 nb_base_functions /= 3;
372#ifndef NDEBUG
373 if (data.getN().size2() % 3) {
374 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
375 "Number of base functions is not divisible by 3");
376 }
377#endif
378 break;
379 case L2:
380 break;
381 default:
382 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED,
383 "Space %s not implemented", FieldSpaceNames[data.getSpace()]);
384 }
385 }
386 return nb_base_functions;
387}
388
389template <AssemblyType A, typename EleOp>
391OpBaseImpl<A, EleOp>::doWork(int row_side, int col_side, EntityType row_type,
392 EntityType col_type,
394 EntitiesFieldData::EntData &col_data) {
396
397 if (entsPtr) {
398 if (entsPtr->find(this->getFEEntityHandle()) == entsPtr->end())
400 }
401
402 // get number of dofs on row
403 nbRows = row_data.getIndices().size();
404 // if no dofs on row, exit that work, nothing to do here
405 if (!nbRows)
407 rowSide = row_side;
408 rowType = row_type;
409 // get number of dofs on column
410 nbCols = col_data.getIndices().size();
411 // if no dofs on column, exit nothing to do here
412 if (!nbCols)
414 colSide = col_side;
415 colType = col_type;
416 // get number of integration points
417 nbIntegrationPts = EleOp::getGaussPts().size2();
418 // get row base functions
419 nbRowBaseFunctions = getNbOfBaseFunctions(row_data);
420
421 // set size of local entity bock
422 locMat.resize(nbRows, nbCols, false);
423 // clear matrix
424 locMat.clear();
425 // integrate local matrix for entity block
426 CHKERR this->iNtegrate(row_data, col_data);
427
428 // assemble local matrix
429 auto check_if_assemble_transpose = [&] {
430 if (this->sYmm) {
431 if (row_side != col_side || row_type != col_type)
432 return true;
433 else
434 return false;
435 } else if (assembleTranspose) {
436 return true;
437 }
438
439 return false;
440 };
441 CHKERR aSsemble(row_data, col_data, check_if_assemble_transpose());
443}
444
445template <AssemblyType A, typename EleOp>
446MoFEMErrorCode OpBaseImpl<A, EleOp>::doWork(int row_side, EntityType row_type,
447 EntData &row_data) {
449
450 if (entsPtr) {
451 if (entsPtr->find(this->getFEEntityHandle()) == entsPtr->end())
453 }
454
455 // get number of dofs on row
456 nbRows = row_data.getIndices().size();
457 rowSide = row_side;
458 rowType = row_type;
459
460 if (!nbRows)
462 // get number of integration points
463 nbIntegrationPts = EleOp::getGaussPts().size2();
464 // get row base functions
465 nbRowBaseFunctions = getNbOfBaseFunctions(row_data);
466 // resize and clear the right hand side vector
467 locF.resize(nbRows);
468 locF.clear();
469 // integrate local vector
470 CHKERR this->iNtegrate(row_data);
471 // assemble local vector
472 CHKERR this->aSsemble(row_data);
474}
475
476template <AssemblyType A, typename EleOp>
478 EntData &col_data,
479 const bool transpose) {
481
482 if (!this->timeScalingFun.empty())
483 this->locMat *= this->timeScalingFun(this->getFEMethod()->ts_t);
484 if (!this->feScalingFun.empty())
485 this->locMat *= this->feScalingFun(this->getFEMethod());
486
487 // Assemble transpose
488 if (transpose) {
489 this->locMatTranspose.resize(this->locMat.size2(), this->locMat.size1(),
490 false);
491 noalias(this->locMatTranspose) = trans(this->locMat);
492 CHKERR matSetValuesHook(this, col_data, row_data, this->locMatTranspose);
493 }
494
495 if (!this->onlyTranspose) {
496 // assemble local matrix
497 CHKERR matSetValuesHook(this, row_data, col_data, this->locMat);
498 }
499
501}
502
503template <AssemblyType A, typename EleOp>
505 if (!this->timeScalingFun.empty())
506 this->locF *= this->timeScalingFun(this->getFEMethod()->ts_t);
507 if (!this->feScalingFun.empty())
508 this->locF *= this->feScalingFun(this->getFEMethod());
509
510 return VecSetValues<AssemblyTypeSelector<A>>(this->getKSPf(), data,
511 this->locF, ADD_VALUES);
512}
513
514} // namespace MoFEM
515
516/**
517 * \defgroup mofem_forms Forms Integrators
518 *
519 * \brief Classes and functions used to evaluate fields at integration pts,
520 *jacobians, etc..
521 *
522 * \ingroup mofem_forces_and_sources
523 **/
524
525#endif //__FORMS_INTEGRATORS_HPP__
@ 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()
@ 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 ...
@ 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()
#define CHKERR
Inline error check.
IntegrationType
Form integrator integration types.
boost::function< double(const double, const double, const double)> ScalarFun
Scalar function type.
MoFEMErrorCode VecSetValues< EssentialBcStorage >(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Set values to vector in operator.
AssemblyType
[Storage and set boundary conditions]
boost::function< FTensor::Tensor1< double, DIM >( const double, const double, const double)> VectorFun
Vector function type.
@ BLOCK_PRECONDITIONER_SCHUR
#define MOFEM_LOG(channel, severity)
Log.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasVector< int > VectorInt
Definition Types.hpp:67
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
FTensor::Tensor2< FTensor::PackPtr< double *, S >, DIM1, DIM2 > getFTensor2FromArray(MatrixDouble &data, const size_t rr, const size_t cc=0)
boost::function< double(double)> TimeFun
Lambda function used to scale with time.
boost::function< double(const FEMethod *fe_ptr)> FEFun
Lambda function used to scale with time.
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
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.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
boost::function< double()> ConstantFun
Constant function type.
auto getFTensor1FromArray(VectorDouble &data)
Get FTensor1 from array.
constexpr auto field_name
FTensor::Index< 'm', 3 > m
Data on single entity (This is passed as argument to DataOperator::doWork)
FieldApproximationBase & getBase()
Get approximation base.
virtual boost::shared_ptr< MatrixDouble > & getNSharedPtr(const FieldApproximationBase base, const BaseDerivatives derivative)
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
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.
boost::function< MoFEMErrorCode( ForcesAndSourcesCore::UserDataOperator *op_ptr, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, MatrixDouble &m)> MatSetValuesHook
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()
int nbCols
number if dof on column
int nbRowBaseFunctions
number or row base functions
EntityType rowType
row type
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.
OpSetBc(std::string field_name, bool yes_set, boost::shared_ptr< std::vector< unsigned char > > boundary_marker)
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.
OpUnSetBc(std::string field_name)