v0.9.2
ElasticOps.hpp
Go to the documentation of this file.
1 /* This file is part of MoFEM.
2  * MoFEM is free software: you can redistribute it and/or modify it under
3  * the terms of the GNU Lesser General Public License as published by the
4  * Free Software Foundation, either version 3 of the License, or (at your
5  * option) any later version.
6  *
7  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
8  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
9  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
10  * License for more details.
11  *
12  * You should have received a copy of the GNU Lesser General Public
13  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
14 
15 namespace OpElasticTools {
16 
17 //! [Common data]
18 struct CommonData {
20  boost::shared_ptr<MatrixDouble> mGradPtr;
21  boost::shared_ptr<MatrixDouble> mStrainPtr;
22  boost::shared_ptr<MatrixDouble> mStressPtr;
23 };
24 //! [Common data]
25 
30 
31 //! [Operators definitions]
32 typedef boost::function<FTensor::Tensor1<double, 2>(const double, const double)>
34 
35 struct OpStrain : public DomainEleOp {
36  OpStrain(const std::string field_name,
37  boost::shared_ptr<CommonData> common_data_ptr);
38  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
39 
40 private:
41  boost::shared_ptr<CommonData> commonDataPtr;
42 };
43 
44 struct OpStress : public DomainEleOp {
45  OpStress(const std::string field_name,
46  boost::shared_ptr<CommonData> common_data_ptr);
47  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
48 
49 private:
50  boost::shared_ptr<CommonData> commonDataPtr;
51 };
52 
54  OpInternalForceRhs(const std::string field_name,
55  boost::shared_ptr<CommonData> common_data_ptr);
56  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
57 
58 private:
59  boost::shared_ptr<CommonData> commonDataPtr;
60 };
61 
62 struct OpForceRhs : public DomainEleOp {
63  OpForceRhs(const std::string field_name,
64  boost::shared_ptr<CommonData> common_data_ptr,
65  VectorFun body_force);
66  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
67 
68 private:
70  boost::shared_ptr<CommonData> commonDataPtr;
71 };
72 
74  OpStiffnessMatrixLhs(const std::string row_field_name,
75  const std::string col_field_name,
76  boost::shared_ptr<CommonData> common_data_ptr);
77  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
78  EntityType col_type, EntData &row_data,
79  EntData &col_data);
80 
81 private:
83  boost::shared_ptr<CommonData> commonDataPtr;
84 };
85 
86 struct OpPostProcElastic : public DomainEleOp {
87  OpPostProcElastic(const std::string field_name,
88  moab::Interface &post_proc_mesh,
89  std::vector<EntityHandle> &map_gauss_pts,
90  boost::shared_ptr<CommonData> common_data_ptr);
91  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
92 
93 private:
95  std::vector<EntityHandle> &mapGaussPts;
96  boost::shared_ptr<CommonData> commonDataPtr;
97 };
98 //! [Operators definitions]
99 
100 OpStrain::OpStrain(const std::string field_name,
101  boost::shared_ptr<CommonData> common_data_ptr)
102  : DomainEleOp(field_name, DomainEleOp::OPROW),
103  commonDataPtr(common_data_ptr) {
104  // Opetor is only executed for vertices
105  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
106 }
107 
108 //! [Calculate strain]
109 MoFEMErrorCode OpStrain::doWork(int side, EntityType type, EntData &data) {
111  const size_t nb_gauss_pts = commonDataPtr->mGradPtr->size2();
112  commonDataPtr->mStrainPtr->resize(3, nb_gauss_pts);
113  auto t_grad = getFTensor2FromMat<2, 2>(*(commonDataPtr->mGradPtr));
114  auto t_strain = getFTensor2SymmetricFromMat<2>(*(commonDataPtr->mStrainPtr));
115 
116  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
117  t_strain(i, j) = (t_grad(i, j) || t_grad(j, i)) / 2;
118  ++t_grad;
119  ++t_strain;
120  }
121 
123 }
124 //! [Calculate strain]
125 OpStress::OpStress(const std::string field_name,
126  boost::shared_ptr<CommonData> common_data_ptr)
127  : DomainEleOp(field_name, DomainEleOp::OPROW),
128  commonDataPtr(common_data_ptr) {
129  // Opetor is only executed for vertices
130  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
131 }
132 
133 //! [Calculate stress]
134 MoFEMErrorCode OpStress::doWork(int side, EntityType type, EntData &data) {
136  const size_t nb_gauss_pts = commonDataPtr->mStrainPtr->size2();
137  commonDataPtr->mStressPtr->resize(3, nb_gauss_pts);
138  auto &t_D = commonDataPtr->tD;
139  auto t_strain = getFTensor2SymmetricFromMat<2>(*(commonDataPtr->mStrainPtr));
140  auto t_stress = getFTensor2SymmetricFromMat<2>(*(commonDataPtr->mStressPtr));
141 
142  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
143  t_stress(i, j) = t_D(i, j, k, l) * t_strain(k, l);
144  ++t_strain;
145  ++t_stress;
146  }
147 
149 }
150 //! [Calculate stress]
151 
153  const std::string field_name, boost::shared_ptr<CommonData> common_data_ptr)
154  : DomainEleOp(field_name, DomainEleOp::OPROW),
155  commonDataPtr(common_data_ptr) {}
156 
157 //! [Internal force]
159  EntData &data) {
161 
162  const size_t nb_dofs = data.getIndices().size();
163  if (nb_dofs) {
164 
165  const size_t nb_base_functions = data.getN().size2();
166  if (2 * nb_base_functions < nb_dofs)
167  SETERRQ(
168  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
169  "Number of DOFs is larger than number of base functions on entity");
170 
171  const size_t nb_gauss_pts = data.getN().size1();
172  std::array<double, MAX_DOFS_ON_ENTITY> nf;
173  std::fill(&nf[0], &nf[nb_dofs], 0);
174 
175  auto t_w = getFTensor0IntegrationWeight();
176  auto t_stress =
177  getFTensor2SymmetricFromMat<2>(*(commonDataPtr->mStressPtr));
178  auto t_diff_base = data.getFTensor1DiffN<2>();
179 
180  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
181 
182  double alpha = getMeasure() * t_w;
183  FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2> t_nf{&nf[0], &nf[1]};
184 
185  size_t bb = 0;
186  for (; bb != nb_dofs / 2; ++bb) {
187  t_nf(i) += alpha * t_diff_base(j) * t_stress(i, j);
188  ++t_diff_base;
189  ++t_nf;
190  }
191  for (; bb < nb_base_functions; ++bb)
192  ++t_diff_base;
193 
194  ++t_stress;
195  ++t_w;
196  }
197 
198  CHKERR VecSetValues(getKSPf(), data, nf.data(), ADD_VALUES);
199  }
200 
202 }
203 //! [Internal force]
204 
205 OpForceRhs::OpForceRhs(const std::string field_name,
206  boost::shared_ptr<CommonData> common_data_ptr,
207  VectorFun body_force)
208  : DomainEleOp(field_name, DomainEleOp::OPROW),
209  commonDataPtr(common_data_ptr), funForce(body_force) {}
210 
211 //! [Body force]
212 MoFEMErrorCode OpForceRhs::doWork(int side, EntityType type, EntData &data) {
214 
215  const size_t nb_dofs = data.getIndices().size();
216  if (nb_dofs) {
217 
218  const size_t nb_base_functions = data.getN().size2();
219  if (2 * nb_base_functions < nb_dofs)
220  SETERRQ(
221  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
222  "Number of DOFs is larger than number of base functions on entity");
223 
224  const size_t nb_gauss_pts = data.getN().size1();
225  std::array<double, MAX_DOFS_ON_ENTITY> nf;
226  std::fill(&nf[0], &nf[nb_dofs], 0);
227 
228  auto t_w = getFTensor0IntegrationWeight();
229  auto t_base = data.getFTensor0N();
230  auto t_coords = getFTensor1CoordsAtGaussPts();
231 
232  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
233 
234  double alpha = getMeasure() * t_w;
235  auto t_force = funForce(t_coords(0), t_coords(1));
236 
237  FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2> t_nf{&nf[0], &nf[1]};
238  size_t bb = 0;
239  for (; bb != nb_dofs / 2; ++bb) {
240  t_nf(i) += alpha * t_base * t_force(i);
241  ++t_base;
242  ++t_nf;
243  }
244  for (; bb < nb_base_functions; ++bb)
245  ++t_base;
246 
247  ++t_w;
248  ++t_coords;
249  }
250 
251  if ((getDataCtx() & PetscData::CtxSetTime).any())
252  for (int dd = 0; dd != nb_dofs; ++dd)
253  nf[dd] *= getTStime();
254 
255  CHKERR VecSetValues(getKSPf(), data, nf.data(), ADD_VALUES);
256  }
257 
259 }
260 //! [Body force]
261 
263  const std::string row_field_name, const std::string col_field_name,
264  boost::shared_ptr<CommonData> common_data_ptr)
265  : DomainEleOp(row_field_name, col_field_name, DomainEleOp::OPROWCOL),
266  commonDataPtr(common_data_ptr) {
267  sYmm = false;
268 }
269 
270 //! [Stiffness]
271 MoFEMErrorCode OpStiffnessMatrixLhs::doWork(int row_side, int col_side,
272  EntityType row_type,
273  EntityType col_type,
274  EntData &row_data,
275  EntData &col_data) {
277 
278  const size_t nb_row_dofs = row_data.getIndices().size();
279  const size_t nb_col_dofs = col_data.getIndices().size();
280 
281  if (nb_row_dofs && nb_col_dofs) {
282 
283  locK.resize(nb_row_dofs, nb_col_dofs, false);
284 
285  const size_t nb_integration_pts = row_data.getN().size1();
286  const size_t nb_row_base_funcs = row_data.getN().size2();
287  auto t_row_diff_base = row_data.getFTensor1DiffN<2>();
288  auto t_w = getFTensor0IntegrationWeight();
289  auto &t_D = commonDataPtr->tD;
290 
291  locK.clear();
292  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
293  double alpha = getMeasure() * t_w;
294 
295  size_t rr = 0;
296  for (; rr != nb_row_dofs / 2; ++rr) {
297 
299 
300  &locK(2 * rr + 0, 0), &locK(2 * rr + 0, 1),
301 
302  &locK(2 * rr + 1, 0), &locK(2 * rr + 1, 1)};
303  auto t_col_diff_base = col_data.getFTensor1DiffN<2>(gg, 0);
304 
305  for (size_t cc = 0; cc != nb_col_dofs / 2; ++cc) {
306  t_a(i, k) += alpha * (t_D(i, j, k, l) *
307  (t_row_diff_base(j) * t_col_diff_base(l)));
308  ++t_col_diff_base;
309  ++t_a;
310  }
311 
312  ++t_row_diff_base;
313  }
314  for (; rr != nb_row_base_funcs; ++rr)
315  ++t_row_diff_base;
316 
317  ++t_w;
318  }
319 
320  CHKERR MatSetValues(getKSPB(), row_data, col_data, &locK(0, 0), ADD_VALUES);
321  }
322 
324 }
325 //! [Stiffness]
326 
328  const std::string field_name, moab::Interface &post_proc_mesh,
329  std::vector<EntityHandle> &map_gauss_pts,
330  boost::shared_ptr<CommonData> common_data_ptr)
331  : DomainEleOp(field_name, DomainEleOp::OPROW), postProcMesh(post_proc_mesh),
332  mapGaussPts(map_gauss_pts), commonDataPtr(common_data_ptr) {
333  // Opetor is only executed for vertices
334  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
335 }
336 
337 //! [Postprocessing]
339  EntData &data) {
341 
342  auto get_tag = [&](const std::string name) {
343  std::array<double, 9> def;
344  std::fill(def.begin(), def.end(), 0);
345  Tag th;
346  CHKERR postProcMesh.tag_get_handle(name.c_str(), 9, MB_TYPE_DOUBLE, th,
347  MB_TAG_CREAT | MB_TAG_SPARSE,
348  def.data());
349  return th;
350  };
351 
352  MatrixDouble3by3 mat(3, 3);
353 
354  auto set_matrix = [&](auto &t) -> MatrixDouble3by3 & {
355  mat.clear();
356  for (size_t r = 0; r != 2; ++r)
357  for (size_t c = 0; c != 2; ++c)
358  mat(r, c) = t(r, c);
359  return mat;
360  };
361 
362  auto set_matrix_symm = [&](auto &t) -> MatrixDouble3by3 & {
363  mat.clear();
364  for (size_t r = 0; r != 2; ++r)
365  for (size_t c = 0; c != 2; ++c)
366  mat(r, c) = t(r, c);
367  return mat;
368  };
369 
370  auto set_plain_stress_strain = [&](auto &mat, auto &t) -> MatrixDouble3by3 & {
371  mat(2, 2) = -poisson_ratio * (t(0, 0) + t(1, 1));
372  return mat;
373  };
374 
375  auto set_tag = [&](auto th, auto gg, MatrixDouble3by3 &mat) {
376  return postProcMesh.tag_set_data(th, &mapGaussPts[gg], 1,
377  &*mat.data().begin());
378  };
379 
380  auto th_grad = get_tag("GRAD");
381  auto th_strain = get_tag("STRAIN");
382  auto th_stress = get_tag("STRESS");
383 
384  size_t nb_gauss_pts = data.getN().size1();
385  auto t_grad = getFTensor2FromMat<2, 2>(*(commonDataPtr->mGradPtr));
386  auto t_strain = getFTensor2SymmetricFromMat<2>(*(commonDataPtr->mStrainPtr));
387  auto t_stress = getFTensor2SymmetricFromMat<2>(*(commonDataPtr->mStressPtr));
388 
389  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
390  CHKERR set_tag(th_grad, gg, set_matrix(t_grad));
391  CHKERR set_tag(
392  th_strain, gg,
393  set_plain_stress_strain(set_matrix_symm(t_strain), t_stress));
394  CHKERR set_tag(th_stress, gg, set_matrix_symm(t_stress));
395  ++t_grad;
396  ++t_strain;
397  ++t_stress;
398  }
399 
401 }
402 //! [Postprocessing]
403 
404 
406 
407  Range forceEdges;
409 
410  OpEdgeForceRhs(const std::string field_name, const Range &force_edges, const VectorDouble &force_vec)
411  : BoundaryEleOp(field_name, OPROW),
412  forceEdges(force_edges), forceVec(force_vec) {}
413 
414 
415  MoFEMErrorCode doWork(int side, EntityType type,
417 
419  const int nb_dofs = data.getIndices().size();
420  if (nb_dofs == 0)
422 
424  if (forceEdges.find(ent) == forceEdges.end()) {
426  }
427 
428  std::array<double, MAX_DOFS_ON_ENTITY> nF;
429  std::fill(&nF[0], &nF[nb_dofs], 0);
430 
432  const int nb_gauss_pts = data.getN().size1();
433  auto t_w = getFTensor0IntegrationWeight();
434  auto t_base = data.getFTensor0N();
435 
436  for (int gg = 0; gg != nb_gauss_pts; ++gg) {
437 
438  FTensor::Tensor1<FTensor::PackPtr<double *, 2>, 2> t_nf(&nF[0], &nF[1]);
439  for (int bb = 0; bb != nb_dofs / 2; ++bb) {
440  t_nf(i) += (t_w * t_base * getMeasure()) * t_force(i);
441  ++t_nf;
442  ++t_base;
443  }
444  ++t_w;
445  }
446 
447  if ((getDataCtx() & PetscData::CtxSetTime).any())
448  for (int dd = 0; dd != nb_dofs; ++dd)
449  nF[dd] *= getTStime();
450 
451  CHKERR VecSetValues(getKSPf(), data, nF.data(), ADD_VALUES);
452 
454  }
455 };
456 
457 
458 }; // namespace OpElasticTools
OpStrain(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
[Operators definitions]
Definition: ElasticOps.hpp:100
FTensor::Index< 'j', 2 > j
Definition: ElasticOps.hpp:27
OpForceRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr, VectorFun body_force)
[Internal force]
Definition: ElasticOps.hpp:205
boost::shared_ptr< MatrixDouble > mStrainPtr
Definition: ElasticOps.hpp:21
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Internal force]
Definition: ElasticOps.hpp:158
boost::shared_ptr< CommonData > commonDataPtr
Definition: ElasticOps.hpp:50
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
boost::shared_ptr< MatrixDouble > mStressPtr
Definition: ElasticOps.hpp:22
const PetscData::Switches & getDataCtx() const
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate stress]
Definition: ElasticOps.hpp:134
std::vector< EntityHandle > & mapGaussPts
Definition: ElasticOps.hpp:95
constexpr double poisson_ratio
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
OpPostProcElastic(const std::string field_name, moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, boost::shared_ptr< CommonData > common_data_ptr)
[Stiffness]
Definition: ElasticOps.hpp:327
const VectorInt & getIndices() const
Get global indices of dofs on entity.
bool sYmm
If true assume that matrix is symmetric structure.
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
boost::shared_ptr< CommonData > commonDataPtr
Definition: ElasticOps.hpp:96
DataForcesAndSourcesCore::EntData EntData
FTensor::Index< 'k', 2 > k
Definition: ElasticOps.hpp:28
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
OpStiffnessMatrixLhs(const std::string row_field_name, const std::string col_field_name, boost::shared_ptr< CommonData > common_data_ptr)
[Body force]
Definition: ElasticOps.hpp:262
auto getFTensor0IntegrationWeight()
Get integration weights.
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
boost::shared_ptr< MatrixDouble > mGradPtr
Definition: ElasticOps.hpp:20
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ElasticOps.hpp:26
OpEdgeForceRhs(const std::string field_name, const Range &force_edges, const VectorDouble &force_vec)
Definition: ElasticOps.hpp:410
MatrixBoundedArray< double, 9 > MatrixDouble3by3
Definition: Types.hpp:101
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Body force]
Definition: ElasticOps.hpp:212
#define CHKERR
Inline error check.
Definition: definitions.h:602
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Calculate strain]
Definition: ElasticOps.hpp:109
boost::shared_ptr< CommonData > commonDataPtr
Definition: ElasticOps.hpp:41
FTensor::Ddg< double, 2, 2 > tD
Definition: ElasticOps.hpp:19
boost::shared_ptr< CommonData > commonDataPtr
Definition: ElasticOps.hpp:59
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Definition: ElasticOps.hpp:415
boost::function< FTensor::Tensor1< double, 2 >const double, const double)> VectorFun
[Operators definitions]
Definition: ElasticOps.hpp:33
FTensor::Index< 'l', 2 > l
Definition: ElasticOps.hpp:29
Definition: ElasticOps.hpp:86
boost::shared_ptr< CommonData > commonDataPtr
Definition: ElasticOps.hpp:83
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
[Postprocessing]
Definition: ElasticOps.hpp:338
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
moab::Interface & postProcMesh
Definition: ElasticOps.hpp:94
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:73
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
OpInternalForceRhs(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
[Calculate stress]
Definition: ElasticOps.hpp:152
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
[Stiffness]
Definition: ElasticOps.hpp:271
OpStress(const std::string field_name, boost::shared_ptr< CommonData > common_data_ptr)
[Calculate strain]
Definition: ElasticOps.hpp:125
boost::shared_ptr< CommonData > commonDataPtr
Definition: ElasticOps.hpp:70