v0.13.1
FormsIntegrators.hpp
Go to the documentation of this file.
1/** \file FormsIntegrators.hpp
2 * \brief Forms inteegrators
3 * \ingroup mofem_form
4
5*/
6
7
8
9#ifndef __FORMS_INTEGRATORS_HPP__
10#define __FORMS_INTEGRATORS_HPP__
11
12namespace MoFEM {
13
14//! [Storage and set boundary conditions]
15
20 map<std::string, std::vector<boost::shared_ptr<EssentialBcStorage>>>;
22};
23
24/**
25 * @brief Set indices on entities on finite element
26 * @ingroup mofem_forms
27 *
28 * If index is marked, its value is set to -1. DOF with such index is not
29 * assembled into the system.
30 *
31 * Indices are stored on the entity.
32 *
33 */
35 OpSetBc(std::string field_name, bool yes_set,
36 boost::shared_ptr<std::vector<unsigned char>> boundary_marker);
39
40public:
41 bool yesSet;
42 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
43};
44
46 OpUnSetBc(std::string field_name);
49};
50
51/**
52 * @brief Set values to vector in operator
53 * @ingroup mofem_forms
54 *
55 * @param V
56 * @param data
57 * @param ptr
58 * @param iora
59 * @return MoFEMErrorCode
60 */
61template <>
65 const double *ptr, InsertMode iora);
66
67/**
68 * @brief Set values to matrix in operator
69 *
70 * @param M
71 * @param row_data
72 * @param col_data
73 * @param ptr
74 * @param iora
75 * @return MoFEMErrorCode
76 */
77template <>
79 Mat M, const EntitiesFieldData::EntData &row_data,
80 const EntitiesFieldData::EntData &col_data, const double *ptr,
81 InsertMode iora);
82
83//! [Storage and set boundary conditions]
84
85/**
86 * @brief Form integrator assembly types
87 * @ingroup mofem_forms
88 *
89 */
91
92/**
93 * @brief Form integrator integration types
94 * @ingroup mofem_forms
95 *
96 */
98
99/**
100 * @brief Scalar function type
101 * @ingroup mofem_forms
102 *
103 */
105 boost::function<double(const double, const double, const double)>;
106
107inline double scalar_fun_one(const double, const double, const double) {
108 return 1;
109}
110
111/**
112 * @brief Lambda function used to scale with time
113 *
114 */
115using TimeFun = boost::function<double(double)>;
116
117/**
118 * @brief Lambda function used to scale with time
119 *
120 */
121using FEFun = boost::function<double(const FEMethod *fe_ptr)>;
122
123/**
124 * @brief Constant function type
125 *
126 */
127using ConstantFun = boost::function<double()>;
128
129/**
130 * @brief Vector function type
131 * @ingroup mofem_forms
132 *
133 * @tparam DIM dimension of the return
134 */
135template <int DIM>
136using VectorFun = boost::function<FTensor::Tensor1<double, DIM>(
137 const double, const double, const double)>;
138
139template <AssemblyType A, typename EleOp> struct OpBaseImpl : public EleOp {
140 using OpType = typename EleOp::OpType;
142
143 OpBaseImpl(const std::string row_field_name, const std::string col_field_name,
144 const OpType type, boost::shared_ptr<Range> ents_ptr = nullptr)
145 : EleOp(row_field_name, col_field_name, type, false),
146 assembleTranspose(false), onlyTranspose(false), entsPtr(ents_ptr) {}
147
148 /**
149 * \brief Do calculations for the left hand side
150 * @param row_side row side number (local number) of entity on element
151 * @param col_side column side number (local number) of entity on element
152 * @param row_type type of row entity
153 * @param col_type type of column entity
154 * @param row_data data for row
155 * @param col_data data for column
156 * @return error code
157 */
158 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
159 EntityType col_type, EntData &row_data,
160 EntData &col_data);
161
162 /**
163 * @brief Do calculations for the right hand side
164 *
165 * @param row_side
166 * @param row_type
167 * @param row_data
168 * @return MoFEMErrorCode
169 */
170 MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
171
172 TimeFun timeScalingFun; ///< assumes that time variable is set
173 FEFun feScalingFun; ///< assumes that time variable is set
174 boost::shared_ptr<Range> entsPtr; ///< Entities on which element is run
175
176protected:
177 template <int DIM>
179 return getFTensor1FromArray<DIM, DIM>(locF);
180 }
181
182 template <int DIM>
184 getLocMat(const int rr) {
185 return getFTensor2FromArray<DIM, DIM, DIM>(locMat, rr);
186 }
187
188 int nbRows; ///< number of dofs on rows
189 int nbCols; ///< number if dof on column
190 int nbIntegrationPts; ///< number of integration points
191 int nbRowBaseFunctions; ///< number or row base functions
192
193 int rowSide; ///< row side number
194 int colSide; ///< column side number
195 EntityType rowType; ///< row type
196 EntityType colType; ///< column type
197
200
201 MatrixDouble locMat; ///< local entity block matrix
202 MatrixDouble locMatTranspose; ///< local entity block matrix
203 VectorDouble locF; ///< local entity vector
204
205 /**
206 * \brief Integrate grad-grad operator
207 * @param row_data row data (consist base functions on row entity)
208 * @param col_data column data (consist base functions on column entity)
209 * @return error code
210 */
211 virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data) {
213 }
214
215 virtual MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data,
216 const bool trans) = 0;
217
218 /**
219 * \brief Class dedicated to integrate operator
220 * @param data entity data on element row
221 * @return error code
222 */
225 }
226
227 virtual MoFEMErrorCode aSsemble(EntData &data) = 0;
228};
229
230/**
231 * @brief Integrator forms
232 * @ingroup mofem_forms
233 *
234 * @tparam EleOp
235 */
236template <typename EleOp> struct FormsIntegrators {
237
239 using OpType = typename EleOp::OpType;
240
241 /**
242 * @brief Assembly methods
243 * @ingroup mofem_forms
244 *
245 * @tparam A
246 */
247 template <AssemblyType A> struct Assembly {
248
250
251 /**
252 * @brief Linear form
253 * @ingroup mofem_forms
254 *
255 * @tparam I
256 */
257 template <IntegrationType I> struct LinearForm;
258
259 /**
260 * @brief Bi linear form
261 * @ingroup mofem_forms
262 *
263 * @tparam I
264 */
265 template <IntegrationType I> struct BiLinearForm;
266
267 }; // Assembly
268}; // namespace MoFEM
269
270template <AssemblyType A, typename EleOp>
272OpBaseImpl<A, EleOp>::doWork(int row_side, int col_side, EntityType row_type,
273 EntityType col_type,
275 EntitiesFieldData::EntData &col_data) {
277
278 if (entsPtr) {
279 if (entsPtr->find(this->getFEEntityHandle()) == entsPtr->end())
281 }
282
283 // get number of dofs on row
284 nbRows = row_data.getIndices().size();
285 // if no dofs on row, exit that work, nothing to do here
286 if (!nbRows)
288 rowSide = row_side;
289 rowType = row_type;
290 // get number of dofs on column
291 nbCols = col_data.getIndices().size();
292 // if no dofs on Columbia, exit nothing to do here
293 if (!nbCols)
295 colSide = col_side;
296 colType = col_type;
297 // get number of integration points
298 nbIntegrationPts = EleOp::getGaussPts().size2();
299 // get row base functions
300 nbRowBaseFunctions = row_data.getN().size2();
301 // set size of local entity bock
302 locMat.resize(nbRows, nbCols, false);
303 // clear matrix
304 locMat.clear();
305 // integrate local matrix for entity block
306 CHKERR this->iNtegrate(row_data, col_data);
307
308 // assemble local matrix
309 auto check_if_assemble_transpose = [&] {
310 if (this->sYmm) {
311 if (row_side != col_side || row_type != col_type)
312 return true;
313 else
314 return false;
315 } else if (assembleTranspose) {
316 return true;
317 }
318
319 return false;
320 };
321 CHKERR aSsemble(row_data, col_data, check_if_assemble_transpose());
323}
324
325template <AssemblyType A, typename EleOp>
327 EntData &row_data) {
329
330 if (entsPtr) {
331 if (entsPtr->find(this->getFEEntityHandle()) == entsPtr->end())
333 }
334
335 // get number of dofs on row
336 nbRows = row_data.getIndices().size();
337 rowSide = row_side;
338 rowType = row_type;
339
340 if (!nbRows)
342 // get number of integration points
343 nbIntegrationPts = EleOp::getGaussPts().size2();
344 // get row base functions
345 nbRowBaseFunctions = row_data.getN().size2();
346 // resize and clear the right hand side vector
347 locF.resize(nbRows);
348 locF.clear();
349 // integrate local vector
350 CHKERR this->iNtegrate(row_data);
351 // assemble local vector
352 CHKERR this->aSsemble(row_data);
354}
355
356template <typename EleOp>
357struct OpBaseImpl<PETSC, EleOp> : public OpBaseImpl<LAST_ASSEMBLE, EleOp> {
359
360protected:
363 const bool transpose) {
365
366 if (!this->timeScalingFun.empty())
367 this->locMat *= this->timeScalingFun(this->getFEMethod()->ts_t);
368 if (!this->feScalingFun.empty())
369 this->locMat *= this->feScalingFun(this->getFEMethod());
370
371 // Assemble transpose
372 if (transpose) {
373 this->locMatTranspose.resize(this->locMat.size2(), this->locMat.size1(),
374 false);
375 noalias(this->locMatTranspose) = trans(this->locMat);
377 this->getKSPB(), col_data, row_data,
378 &*this->locMatTranspose.data().begin(), ADD_VALUES);
379 }
380
381 if (!this->onlyTranspose) {
382 // assemble local matrix
384 this->getKSPB(), row_data, col_data, &*this->locMat.data().begin(),
385 ADD_VALUES);
386 }
387
389 }
390
392
393 if (!this->timeScalingFun.empty())
394 this->locF *= this->timeScalingFun(this->getFEMethod()->ts_t);
395 if (!this->feScalingFun.empty())
396 this->locF *= this->feScalingFun(this->getFEMethod());
397
399 this->getKSPf(), data, &*this->locF.data().begin(), ADD_VALUES);
400 }
401};
402
403} // namespace MoFEM
404
405/**
406 * \defgroup mofem_forms Forms Integrators
407 *
408 * \brief Classes and functions used to evaluate fields at integration pts,
409 *jacobians, etc..
410 *
411 * \ingroup mofem_forces_and_sources
412 **/
413
414#endif //__FORMS_INTEGRATORS_HPP__
static Index< 'M', 3 > M
EntitiesFieldData::EntData EntData
#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.
const FTensor::Tensor2< T, Dim, Dim > Vec
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: MoFEM.hpp:24
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
EssentialBcStorage(VectorInt &indices)
structure for User Loop Methods on finite elements
typename EleOp::OpType OpType
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data, const bool transpose)
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)
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.
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 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 MoFEMErrorCode aSsemble(EntData &data)=0
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
virtual MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data, const bool trans)=0
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)