v0.13.0
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 /* This file is part of MoFEM.
8  * MoFEM is free software: you can redistribute it and/or modify it under
9  * the terms of the GNU Lesser General Public License as published by the
10  * Free Software Foundation, either version 3 of the License, or (at your
11  * option) any later version.
12  *
13  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
14  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  * License for more details.
17  *
18  * You should have received a copy of the GNU Lesser General Public
19  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
20 
21 #ifndef __FORMS_INTEGRATORS_HPP__
22 #define __FORMS_INTEGRATORS_HPP__
23 
24 namespace MoFEM {
25 
26 //! [Storage and set boundary conditions]
27 
29  EssentialBcStorage(VectorInt &indices) : entityIndices(indices) {}
32  map<std::string, std::vector<boost::shared_ptr<EssentialBcStorage>>>;
34 };
35 
36 /**
37  * @brief Set indices on entities on finite element
38  * @ingroup mofem_forms
39  *
40  * If index is marked, its value is set to -1. DOF with such index is not
41  * assembled into the system.
42  *
43  * Indices are stored on the entity.
44  *
45  */
47  OpSetBc(std::string field_name, bool yes_set,
48  boost::shared_ptr<std::vector<unsigned char>> boundary_marker);
49  MoFEMErrorCode doWork(int side, EntityType type,
51 
52 public:
53  bool yesSet;
54  boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
55 };
56 
58  OpUnSetBc(std::string field_name);
59  MoFEMErrorCode doWork(int side, EntityType type,
61 };
62 
63 /**
64  * @brief Set values to vector in operator
65  * @ingroup mofem_forms
66  *
67  * @param V
68  * @param data
69  * @param ptr
70  * @param iora
71  * @return MoFEMErrorCode
72  */
73 template <>
76  const EntitiesFieldData::EntData &data,
77  const double *ptr, InsertMode iora);
78 
79 /**
80  * @brief Set values to matrix in operator
81  *
82  * @param M
83  * @param row_data
84  * @param col_data
85  * @param ptr
86  * @param iora
87  * @return MoFEMErrorCode
88  */
89 template <>
91  Mat M, const EntitiesFieldData::EntData &row_data,
92  const EntitiesFieldData::EntData &col_data, const double *ptr,
93  InsertMode iora);
94 
95 //! [Storage and set boundary conditions]
96 
97 /**
98  * @brief Form integrator assembly types
99  * @ingroup mofem_forms
100  *
101  */
103 
104 /**
105  * @brief Form integrator integration types
106  * @ingroup mofem_forms
107  *
108  */
110 
111 /**
112  * @brief Scalar function type
113  * @ingroup mofem_forms
114  *
115  */
116 using ScalarFun =
117  boost::function<double(const double, const double, const double)>;
118 
119 /**
120  * @brief Constant function type
121  *
122  */
123 using ConstantFun = boost::function<double()>;
124 
125 /**
126  * @brief Vector function type
127  * @ingroup mofem_forms
128  *
129  * @tparam DIM dimension of the return
130  */
131 template <int DIM>
132 using VectorFun = boost::function<FTensor::Tensor1<double, DIM>(
133  const double, const double, const double)>;
134 
135 template <AssemblyType A, typename EleOp> struct OpBaseImpl : public EleOp {
136  using OpType = typename EleOp::OpType;
138 
139  OpBaseImpl(const std::string row_field_name, const std::string col_field_name,
140  const OpType type)
141  : EleOp(row_field_name, col_field_name, type, false),
142  assembleTranspose(false), onlyTranspose(false) {}
143 
144  /**
145  * \brief Do calculations for the left hand side
146  * @param row_side row side number (local number) of entity on element
147  * @param col_side column side number (local number) of entity on element
148  * @param row_type type of row entity
149  * @param col_type type of column entity
150  * @param row_data data for row
151  * @param col_data data for column
152  * @return error code
153  */
154  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
155  EntityType col_type, EntData &row_data,
156  EntData &col_data);
157 
158  /**
159  * @brief Do calculations for the right hand side
160  *
161  * @param row_side
162  * @param row_type
163  * @param row_data
164  * @return MoFEMErrorCode
165  */
166  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
167 
168 protected:
169  template <int DIM>
171  return getFTensor1FromArray<DIM, DIM>(locF);
172  }
173 
174  template <int DIM>
176  getLocMat(const int rr) {
177  return getFTensor2FromArray<DIM, DIM, DIM>(locMat, rr);
178  }
179 
180  int nbRows; ///< number of dofs on rows
181  int nbCols; ///< number if dof on column
182  int nbIntegrationPts; ///< number of integration points
183  int nbRowBaseFunctions; ///< number or row base functions
184 
185  int rowSide; ///< row side number
186  int colSide; ///< column side number
187  EntityType rowType; ///< row type
188  EntityType colType; ///< column type
189 
192 
193  MatrixDouble locMat; ///< local entity block matrix
194  MatrixDouble locMatTranspose; ///< local entity block matrix
195  VectorDouble locF; ///< local entity vector
196 
197  /**
198  * \brief Integrate grad-grad operator
199  * @param row_data row data (consist base functions on row entity)
200  * @param col_data column data (consist base functions on column entity)
201  * @return error code
202  */
203  virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data) {
204  return MOFEM_NOT_IMPLEMENTED;
205  }
206 
207  virtual MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data,
208  const bool trans) = 0;
209 
210  /**
211  * \brief Class dedicated to integrate operator
212  * @param data entity data on element row
213  * @return error code
214  */
216  return MOFEM_NOT_IMPLEMENTED;
217  }
218 
219  virtual MoFEMErrorCode aSsemble(EntData &data) = 0;
220 };
221 
222 /**
223  * @brief Integrator forms
224  * @ingroup mofem_forms
225  *
226  * @tparam EleOp
227  */
228 template <typename EleOp> struct FormsIntegrators {
229 
231  using OpType = typename EleOp::OpType;
232 
233  /**
234  * @brief Assembly methods
235  * @ingroup mofem_forms
236  *
237  * @tparam A
238  */
239  template <AssemblyType A> struct Assembly {
240 
242 
243  /**
244  * @brief Linear form
245  * @ingroup mofem_forms
246  *
247  * @tparam I
248  */
249  template <IntegrationType I> struct LinearForm;
250 
251  /**
252  * @brief Bi linear form
253  * @ingroup mofem_forms
254  *
255  * @tparam I
256  */
257  template <IntegrationType I> struct BiLinearForm;
258 
259  }; // Assembly
260 }; // namespace MoFEM
261 
262 template <AssemblyType A, typename EleOp>
264 OpBaseImpl<A, EleOp>::doWork(int row_side, int col_side, EntityType row_type,
265  EntityType col_type,
266  EntitiesFieldData::EntData &row_data,
267  EntitiesFieldData::EntData &col_data) {
269  // get number of dofs on row
270  nbRows = row_data.getIndices().size();
271  // if no dofs on row, exit that work, nothing to do here
272  if (!nbRows)
274  rowSide = row_side;
275  rowType = row_type;
276  // get number of dofs on column
277  nbCols = col_data.getIndices().size();
278  // if no dofs on Columbia, exit nothing to do here
279  if (!nbCols)
281  colSide = col_side;
282  colType = col_type;
283  // get number of integration points
284  nbIntegrationPts = EleOp::getGaussPts().size2();
285  // get row base functions
286  nbRowBaseFunctions = row_data.getN().size2();
287  // set size of local entity bock
288  locMat.resize(nbRows, nbCols, false);
289  // clear matrix
290  locMat.clear();
291  // integrate local matrix for entity block
292  CHKERR this->iNtegrate(row_data, col_data);
293 
294  // assemble local matrix
295  auto check_if_assemble_transpose = [&] {
296  if (this->sYmm) {
297  if (row_side != col_side || row_type != col_type)
298  return true;
299  else
300  return false;
301  } else if (assembleTranspose) {
302  return true;
303  }
304 
305  return false;
306  };
307  CHKERR aSsemble(row_data, col_data, check_if_assemble_transpose());
309 }
310 
311 template <AssemblyType A, typename EleOp>
312 MoFEMErrorCode OpBaseImpl<A, EleOp>::doWork(int row_side, EntityType row_type,
313  EntData &row_data) {
315  // get number of dofs on row
316  nbRows = row_data.getIndices().size();
317  rowSide = row_side;
318  rowType = row_type;
319 
320  if (!nbRows)
322  // get number of integration points
323  nbIntegrationPts = EleOp::getGaussPts().size2();
324  // get row base functions
325  nbRowBaseFunctions = row_data.getN().size2();
326  // resize and clear the right hand side vector
327  locF.resize(nbRows);
328  locF.clear();
329  // integrate local vector
330  CHKERR this->iNtegrate(row_data);
331  // assemble local vector
332  CHKERR this->aSsemble(row_data);
334 }
335 
336 template <typename EleOp>
337 struct OpBaseImpl<PETSC, EleOp> : public OpBaseImpl<LAST_ASSEMBLE, EleOp> {
339 
340 protected:
342  EntitiesFieldData::EntData &col_data,
343  const bool transpose) {
345 
346  // Assemble transpose
347  if (transpose) {
348  this->locMatTranspose.resize(this->locMat.size2(), this->locMat.size1(),
349  false);
350  noalias(this->locMatTranspose) = trans(this->locMat);
352  this->getKSPB(), col_data, row_data,
353  &*this->locMatTranspose.data().begin(), ADD_VALUES);
354  }
355 
356  if (!this->onlyTranspose) {
357  // assemble local matrix
359  this->getKSPB(), row_data, col_data, &*this->locMat.data().begin(),
360  ADD_VALUES);
361  }
362 
364  }
365 
368  this->getKSPf(), data, &*this->locF.data().begin(), ADD_VALUES);
369  }
370 };
371 
372 } // namespace MoFEM
373 
374 /**
375  * \defgroup mofem_forms Forms Integrators
376  *
377  * \brief Classes and functions used to evaluate fields at integration pts,
378  *jacobians, etc..
379  *
380  * \ingroup mofem_forces_and_sources
381  **/
382 
383 #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:460
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:45
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
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:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
UBlasVector< int > VectorInt
Definition: Types.hpp:78
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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.
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)
typename EleOp::OpType OpType
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data, const bool transpose)
MoFEMErrorCode aSsemble(EntitiesFieldData::EntData &data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Do calculations for the right hand side.
VectorDouble locF
local entity vector
typename EleOp::OpType OpType
MatrixDouble locMatTranspose
local entity block matrix
int rowSide
row side number
int colSide
column side number
int nbRows
number of dofs on rows
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.
EntityType colType
column type
int nbIntegrationPts
number of integration points
virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate grad-grad operator.
FTensor::Tensor1< FTensor::PackPtr< double *, DIM >, DIM > getNf()
virtual MoFEMErrorCode iNtegrate(EntData &data)
Class dedicated to integrate operator.
MatrixDouble locMat
local entity block matrix
OpBaseImpl(const std::string row_field_name, const std::string col_field_name, const OpType type)
virtual MoFEMErrorCode aSsemble(EntData &data)=0
int nbCols
number if dof on column
FTensor::Tensor2< FTensor::PackPtr< double *, DIM >, DIM, DIM > getLocMat(const int rr)
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.
OpSetBc(std::string field_name, bool yes_set, boost::shared_ptr< std::vector< unsigned char >> boundary_marker)
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.
OpUnSetBc(std::string field_name)