v0.13.2
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 enum AssembleTo { AMat, BMat }; ///< Where to assemble
208
209 /**
210 * @brief Where to assemble
211 *
212 * \note Default assembly is BMat, i.e. preconditioner matrix
213 *
214 * @param a_to
215 * @return MoFEMErrorCode
216 */
217 inline void setAssembleTo(AssembleTo a_to) { assembleTo = a_to; };
218
219protected:
221
222 /**
223 * @brief Select matrix
224 *
225 * @return auto
226 */
227 inline auto matSelector() {
228 switch (assembleTo) {
229 case AMat:
230 return this->getKSPA();
231 case BMat:
232 return this->getKSPB();
233 }
234 return this->getKSPB();
235 };
236
237
238 template <int DIM>
240 return getFTensor1FromArray<DIM, DIM>(locF);
241 }
242
243 template <int DIM>
245 getLocMat(const int rr) {
246 return getFTensor2FromArray<DIM, DIM, DIM>(locMat, rr);
247 }
248
249 int nbRows; ///< number of dofs on rows
250 int nbCols; ///< number if dof on column
251 int nbIntegrationPts; ///< number of integration points
252 int nbRowBaseFunctions; ///< number or row base functions
253
254 int rowSide; ///< row side number
255 int colSide; ///< column side number
256 EntityType rowType; ///< row type
257 EntityType colType; ///< column type
258
261
262 MatrixDouble locMat; ///< local entity block matrix
263 MatrixDouble locMatTranspose; ///< local entity block matrix
264 VectorDouble locF; ///< local entity vector
265
266 /**
267 * \brief Integrate grad-grad operator
268 * @param row_data row data (consist base functions on row entity)
269 * @param col_data column data (consist base functions on column entity)
270 * @return error code
271 */
272 virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data) {
274 }
275
276 virtual MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data,
277 const bool trans);
278
279 /**
280 * \brief Class dedicated to integrate operator
281 * @param data entity data on element row
282 * @return error code
283 */
286 }
287
289};
290
291/**
292 * @brief Integrator forms
293 * @ingroup mofem_forms
294 *
295 * @tparam EleOp
296 */
297template <typename EleOp> struct FormsIntegrators {
298
300 using OpType = typename EleOp::OpType;
301
302 /**
303 * @brief Assembly methods
304 * @ingroup mofem_forms
305 *
306 * @tparam A
307 */
308 template <AssemblyType A> struct Assembly {
309
311
312 /**
313 * @brief Linear form
314 * @ingroup mofem_forms
315 *
316 * @tparam I
317 */
318 template <IntegrationType I> struct LinearForm;
319
320 /**
321 * @brief Bi linear form
322 * @ingroup mofem_forms
323 *
324 * @tparam I
325 */
326 template <IntegrationType I> struct BiLinearForm;
327
328 }; // Assembly
329}; // namespace MoFEM
330
331template <AssemblyType A, typename EleOp>
333OpBaseImpl<A, EleOp>::doWork(int row_side, int col_side, EntityType row_type,
334 EntityType col_type,
336 EntitiesFieldData::EntData &col_data) {
338
339 if (entsPtr) {
340 if (entsPtr->find(this->getFEEntityHandle()) == entsPtr->end())
342 }
343
344 // get number of dofs on row
345 nbRows = row_data.getIndices().size();
346 // if no dofs on row, exit that work, nothing to do here
347 if (!nbRows)
349 rowSide = row_side;
350 rowType = row_type;
351 // get number of dofs on column
352 nbCols = col_data.getIndices().size();
353 // if no dofs on column, exit nothing to do here
354 if (!nbCols)
356 colSide = col_side;
357 colType = col_type;
358 // get number of integration points
359 nbIntegrationPts = EleOp::getGaussPts().size2();
360 // get row base functions
361 nbRowBaseFunctions = row_data.getN().size2();
362 // set size of local entity bock
363 locMat.resize(nbRows, nbCols, false);
364 // clear matrix
365 locMat.clear();
366 // integrate local matrix for entity block
367 CHKERR this->iNtegrate(row_data, col_data);
368
369 // assemble local matrix
370 auto check_if_assemble_transpose = [&] {
371 if (this->sYmm) {
372 if (row_side != col_side || row_type != col_type)
373 return true;
374 else
375 return false;
376 } else if (assembleTranspose) {
377 return true;
378 }
379
380 return false;
381 };
382 CHKERR aSsemble(row_data, col_data, check_if_assemble_transpose());
384}
385
386template <AssemblyType A, typename EleOp>
388 EntData &row_data) {
390
391 if (entsPtr) {
392 if (entsPtr->find(this->getFEEntityHandle()) == entsPtr->end())
394 }
395
396 // get number of dofs on row
397 nbRows = row_data.getIndices().size();
398 rowSide = row_side;
399 rowType = row_type;
400
401 if (!nbRows)
403 // get number of integration points
404 nbIntegrationPts = EleOp::getGaussPts().size2();
405 // get row base functions
406 nbRowBaseFunctions = row_data.getN().size2();
407 // resize and clear the right hand side vector
408 locF.resize(nbRows);
409 locF.clear();
410 // integrate local vector
411 CHKERR this->iNtegrate(row_data);
412 // assemble local vector
413 CHKERR this->aSsemble(row_data);
415}
416
417template <AssemblyType A, typename EleOp>
419 EntData &col_data,
420 const bool transpose) {
422
423 if (!this->timeScalingFun.empty())
424 this->locMat *= this->timeScalingFun(this->getFEMethod()->ts_t);
425 if (!this->feScalingFun.empty())
426 this->locMat *= this->feScalingFun(this->getFEMethod());
427
428 // Assemble transpose
429 if (transpose) {
430 this->locMatTranspose.resize(this->locMat.size2(), this->locMat.size1(),
431 false);
432 noalias(this->locMatTranspose) = trans(this->locMat);
433 CHKERR MatSetValues<AssemblyTypeSelector<A>>(
434 matSelector(), col_data, row_data, this->locMatTranspose, ADD_VALUES);
435 }
436
437 if (!this->onlyTranspose) {
438 // assemble local matrix
439 CHKERR MatSetValues<AssemblyTypeSelector<A>>(
440 matSelector(), row_data, col_data, this->locMat, ADD_VALUES);
441 }
442
444}
445
446template <AssemblyType A, typename EleOp>
448 if (!this->timeScalingFun.empty())
449 this->locF *= this->timeScalingFun(this->getFEMethod()->ts_t);
450 if (!this->feScalingFun.empty())
451 this->locF *= this->feScalingFun(this->getFEMethod());
452
453 return VecSetValues<AssemblyTypeSelector<A>>(this->getKSPf(), data,
454 this->locF, ADD_VALUES);
455}
456
457} // namespace MoFEM
458
459/**
460 * \defgroup mofem_forms Forms Integrators
461 *
462 * \brief Classes and functions used to evaluate fields at integration pts,
463 *jacobians, etc..
464 *
465 * \ingroup mofem_forces_and_sources
466 **/
467
468#endif //__FORMS_INTEGRATORS_HPP__
static Index< 'M', 3 > M
#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_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
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
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)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
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
void setAssembleTo(AssembleTo a_to)
Where to assemble.
int colSide
column side number
boost::shared_ptr< Range > entsPtr
Entities on which element is run.
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.
auto matSelector()
Select matrix.
MatrixDouble locMat
local entity block matrix
FTensor::Tensor2< FTensor::PackPtr< double *, DIM >, DIM, DIM > getLocMat(const int rr)
FEFun feScalingFun
assumes that time variable is set
enum AssembleTo assembleTo
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
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.