v0.9.1
HookeElement.hpp
Go to the documentation of this file.
1 /**
2  * \file HookeElement.hpp
3  * \example HookeElement.hpp
4  *
5  * \brief Operators and data structures for linear elastic
6  * analysis
7  *
8  * Implemention of operators for Hooke material. Implementation is extended to
9  * the case when the mesh is moving as results of topological changes, also the
10  * calculation of material forces and associated tangent matrices are added to
11  * implementation.
12  *
13  * In other words, spatial deformation is small but topological changes large.
14  */
15 
16 /*
17  * This file is part of MoFEM.
18  * MoFEM is free software: you can redistribute it and/or modify it under
19  * the terms of the GNU Lesser General Public License as published by the
20  * Free Software Foundation, either version 3 of the License, or (at your
21  * option) any later version.
22  *
23  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
24  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
25  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
26  * License for more details.
27  *
28  * You should have received a copy of the GNU Lesser General Public
29  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
30 
31 #ifndef __HOOKE_ELEMENT_HPP
32 #define __HOOKE_ELEMENT_HPP
33 
34 #ifndef __BASICFINITEELEMENTS_HPP__
35 #include <BasicFiniteElements.hpp>
36 #endif // __BASICFINITEELEMENTS_HPP__
37 
38 #ifndef __NONLINEAR_ELASTIC_HPP
39 
40 struct NonlinearElasticElement {
41 
42  /** \brief data for calculation heat conductivity and heat capacity elements
43  * \ingroup nonlinear_elastic_elem
44  */
45  struct BlockData {
46  int iD;
47  double E;
48  double PoissonRatio;
49  Range tEts; ///< constrains elements in block set
50  Range forcesOnlyOnEntitiesRow;
51  Range forcesOnlyOnEntitiesCol;
52  };
53 };
54 
55 #endif // __NONLINEAR_ELASTIC_HPP
56 
57 /** \brief structure grouping operators and data used for calculation of
58  * nonlinear elastic element \ingroup nonlinear_elastic_elem
59  *
60  * In order to assemble matrices and right hand vectors, the loops over
61  * elements, entities over that elements and finally loop over integration
62  * points are executed.
63  *
64  * Following implementation separate those three categories of loops and to each
65  * loop attach operator.
66  *
67  */
68 
69 #ifndef __CONVECTIVE_MASS_ELEMENT_HPP
70 struct ConvectiveMassElement {
71  /** \brief data for calculation inertia forces
72  * \ingroup user_modules
73  */
74  struct BlockData {
75  double rho0; ///< reference density
76  VectorDouble a0; ///< constant acceleration
77  Range tEts; ///< elements in block set
78  };
79 }
80 
81 #endif //__CONVECTIVE_MASS_ELEMENT_HPP
82 struct HookeElement {
83 
84  using BlockData = NonlinearElasticElement::BlockData;
85  using MassBlockData = ConvectiveMassElement::BlockData;
86 
87  using EntData = DataForcesAndSourcesCore::EntData;
88  using UserDataOperator = ForcesAndSourcesCore::UserDataOperator;
89  using VolUserDataOperator =
90  VolumeElementForcesAndSourcesCore::UserDataOperator;
91 
92  struct DataAtIntegrationPts {
93 
94  boost::shared_ptr<MatrixDouble> smallStrainMat;
95  boost::shared_ptr<MatrixDouble> hMat;
96  boost::shared_ptr<MatrixDouble> FMat;
97 
98  boost::shared_ptr<MatrixDouble> HMat;
99  boost::shared_ptr<VectorDouble> detHVec;
100  boost::shared_ptr<MatrixDouble> invHMat;
101 
102  boost::shared_ptr<MatrixDouble> cauchyStressMat;
103  boost::shared_ptr<MatrixDouble> stiffnessMat;
104  boost::shared_ptr<VectorDouble> energyVec;
105  boost::shared_ptr<MatrixDouble> eshelbyStressMat;
106 
107  boost::shared_ptr<MatrixDouble> eshelbyStress_dx;
108 
109  DataAtIntegrationPts() {
110 
111  smallStrainMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
112  hMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
113  FMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
114 
115  HMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
116  detHVec = boost::shared_ptr<VectorDouble>(new VectorDouble());
117  invHMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
118 
119  cauchyStressMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
120  stiffnessMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
121  energyVec = boost::shared_ptr<VectorDouble>(new VectorDouble());
122  eshelbyStressMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
123  stiffnessMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
124 
125  eshelbyStress_dx = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
126  }
127 
128  Range forcesOnlyOnEntitiesRow;
129  Range forcesOnlyOnEntitiesCol;
130  };
131 
132  template <bool D = false>
133  struct OpCalculateStrain : public VolUserDataOperator {
134 
135  OpCalculateStrain(const std::string row_field, const std::string col_field,
136  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
137 
138  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
139 
140  private:
141  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
142  };
143 
144  struct OpCalculateStrainAle : public VolUserDataOperator {
145 
146  OpCalculateStrainAle(const std::string row_field,
147  const std::string col_field,
148  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
149 
150  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
151 
152  private:
153  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
154  };
155 
156 #define MAT_TO_DDG(SM) \
157  &(*SM)(0, 0), &(*SM)(1, 0), &(*SM)(2, 0), &(*SM)(3, 0), &(*SM)(4, 0), \
158  &(*SM)(5, 0), &(*SM)(6, 0), &(*SM)(7, 0), &(*SM)(8, 0), &(*SM)(9, 0), \
159  &(*SM)(10, 0), &(*SM)(11, 0), &(*SM)(12, 0), &(*SM)(13, 0), \
160  &(*SM)(14, 0), &(*SM)(15, 0), &(*SM)(16, 0), &(*SM)(17, 0), \
161  &(*SM)(18, 0), &(*SM)(19, 0), &(*SM)(20, 0), &(*SM)(21, 0), \
162  &(*SM)(22, 0), &(*SM)(23, 0), &(*SM)(24, 0), &(*SM)(25, 0), \
163  &(*SM)(26, 0), &(*SM)(27, 0), &(*SM)(28, 0), &(*SM)(29, 0), \
164  &(*SM)(30, 0), &(*SM)(31, 0), &(*SM)(32, 0), &(*SM)(33, 0), \
165  &(*SM)(34, 0), &(*SM)(35, 0)
166 
167  template <int S = 0> struct OpCalculateStress : public VolUserDataOperator {
168 
169  OpCalculateStress(const std::string row_field, const std::string col_field,
170  boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
171 
172  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
173 
174  protected:
175  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
176  };
177 
178  struct OpCalculateEnergy : public VolUserDataOperator {
179 
180  OpCalculateEnergy(const std::string row_field, const std::string col_field,
181  boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
182  Vec ghost_vec = PETSC_NULL);
183 
184  ~OpCalculateEnergy();
185 
186  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
187 
188  protected:
189  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
190  Vec ghostVec;
191  };
192 
193  struct OpCalculateEshelbyStress : public VolUserDataOperator {
194 
195  OpCalculateEshelbyStress(
196  const std::string row_field, const std::string col_field,
197  boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
198 
199  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
200 
201  protected:
202  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
203  };
204 
205  template <int S = 0>
206  struct OpCalculateHomogeneousStiffness : public VolUserDataOperator {
207 
208  OpCalculateHomogeneousStiffness(
209  const std::string row_field, const std::string col_field,
210  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
211  boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
212 
213  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
214 
215  protected:
216  boost::shared_ptr<map<int, BlockData>>
217  blockSetsPtr; ///< Structure keeping data about problem, like
218  ///< material parameters
219  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
220  };
221 
222  /** * @brief Assemble mass matrix for elastic element TODO: CHANGE FORMULA *
223  * \f[
224  * {\bf{M}} = \int\limits_\Omega
225  * \f]
226  *
227  */
228  struct OpCalculateMassMatrix
229  : public VolumeElementForcesAndSourcesCore::UserDataOperator {
230 
231  MatrixDouble locK;
232  MatrixDouble translocK;
233  BlockData &dAta;
234  MassBlockData &massData;
235 
236  boost::shared_ptr<DataAtIntegrationPts> commonData;
237 
238  OpCalculateMassMatrix(const std::string row_field,
239  const std::string col_field, BlockData &data,
240  MassBlockData &mass_data,
241  boost::shared_ptr<DataAtIntegrationPts> &common_data,
242  bool symm = true)
243  : VolumeElementForcesAndSourcesCore::UserDataOperator(
244  row_field, col_field, OPROWCOL, symm),
245  commonData(common_data), dAta(data), massData(mass_data) {}
246 
247  PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
248  EntityType col_type,
249  DataForcesAndSourcesCore::EntData &row_data,
250  DataForcesAndSourcesCore::EntData &col_data) {
251  MoFEMFunctionBegin;
252 
253  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
254  return FTensor::Tensor2<double *, 3, 3>(
255  &m(3 * r + 0, 3 * c + 0), &m(3 * r + 0, 3 * c + 1),
256  &m(3 * r + 0, 3 * c + 2), &m(3 * r + 1, 3 * c + 0),
257  &m(3 * r + 1, 3 * c + 1), &m(3 * r + 1, 3 * c + 2),
258  &m(3 * r + 2, 3 * c + 0), &m(3 * r + 2, 3 * c + 1),
259  &m(3 * r + 2, 3 * c + 2));
260  };
261 
262  const int row_nb_dofs = row_data.getIndices().size();
263  if (!row_nb_dofs)
264  MoFEMFunctionReturnHot(0);
265  const int col_nb_dofs = col_data.getIndices().size();
266  if (!col_nb_dofs)
267  MoFEMFunctionReturnHot(0);
268  if (dAta.tEts.find(getFEEntityHandle()) == dAta.tEts.end()) {
269  MoFEMFunctionReturnHot(0);
270  }
271  if (massData.tEts.find(getFEEntityHandle()) == massData.tEts.end()) {
272  MoFEMFunctionReturnHot(0);
273  }
274 
275  const bool diagonal_block =
276  (row_type == col_type) && (row_side == col_side);
277  // get number of integration points
278  // Set size can clear local tangent matrix
279  locK.resize(row_nb_dofs, col_nb_dofs, false);
280  locK.clear();
281 
282  const int row_nb_gauss_pts = row_data.getN().size1();
283  const int row_nb_base_functions = row_data.getN().size2();
284 
285  FTensor::Index<'i', 3> i;
286  FTensor::Index<'j', 3> j;
287  FTensor::Index<'k', 3> k;
288  FTensor::Index<'l', 3> l;
289 
290  double density = massData.rho0;
291 
292  // get integration weights
293  auto t_w = getFTensor0IntegrationWeight();
294 
295  // integrate local matrix for entity block
296  for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
297 
298  auto t_row_base_func = row_data.getFTensor0N(gg, 0);
299 
300  // Get volume and integration weight
301  double w = getVolume() * t_w;
302 
303  for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
304  auto t_col_base_func = col_data.getFTensor0N(gg, 0);
305  for (int col_bb = 0; col_bb != col_nb_dofs / 3; col_bb++) {
306  auto t_assemble = get_tensor2(locK, row_bb, col_bb);
307  t_assemble(i, j) += density * t_row_base_func * t_col_base_func * w;
308  // Next base function for column
309  ++t_col_base_func;
310  }
311  // Next base function for row
312  ++t_row_base_func;
313  }
314  // Next integration point for getting weight
315  ++t_w;
316  }
317 
318  CHKERR MatSetValues(getKSPB(), row_data, col_data, &locK(0, 0),
319  ADD_VALUES);
320 
321  // is symmetric
322  if (row_type != col_type || row_side != col_side) {
323  translocK.resize(col_nb_dofs, row_nb_dofs, false);
324  noalias(translocK) = trans(locK);
325 
326  CHKERR MatSetValues(getKSPB(), col_data, row_data, &translocK(0, 0),
327  ADD_VALUES);
328  }
329 
330  MoFEMFunctionReturn(0);
331  }
332  };
333 
334  struct OpCalculateStiffnessScaledByDensityField : public VolUserDataOperator {
335  protected:
336  boost::shared_ptr<map<int, BlockData>>
337  blockSetsPtr; ///< Structure keeping data about problem, like
338  ///< material parameters
339  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
340 
341  boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
342  const double rhoN; ///< exponent n in E(p) = E * (p / p_0)^n
343  const double rHo0; ///< p_0 reference density in E(p) = E * (p / p_0)^n
344  // // where p is density, E - youngs modulus
345  public:
346  OpCalculateStiffnessScaledByDensityField(
347  const std::string row_field, const std::string col_field,
348  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
349  boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
350  boost::shared_ptr<VectorDouble> rho_at_gauss_pts, const double rho_n,
351  const double rho_0);
352 
353  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
354  };
355 
356  struct OpAssemble : public VolUserDataOperator {
357 
358  OpAssemble(const std::string row_field, const std::string col_field,
359  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
360  const char type, bool symm = false);
361 
362  /**
363  * \brief Do calculations for give operator
364  * @param row_side row side number (local number) of entity on element
365  * @param col_side column side number (local number) of entity on element
366  * @param row_type type of row entity MBVERTEX, MBEDGE, MBTRI or MBTET
367  * @param col_type type of column entity MBVERTEX, MBEDGE, MBTRI or MBTET
368  * @param row_data data for row
369  * @param col_data data for column
370  * @return error code
371  */
372  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
373  EntityType col_type, EntData &row_data,
374  EntData &col_data);
375 
376  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
377 
378  protected:
379  // Finite element stiffness sub-matrix K_ij
380  MatrixDouble K;
381  MatrixDouble transK;
382  VectorDouble nF;
383 
384  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
385 
386  VectorInt rowIndices;
387  VectorInt colIndices;
388 
389  int nbRows; ///< number of dofs on rows
390  int nbCols; ///< number if dof on column
391  int nbIntegrationPts; ///< number of integration points
392  bool isDiag; ///< true if this block is on diagonal
393 
394  virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
395 
396  virtual MoFEMErrorCode iNtegrate(EntData &row_data);
397 
398  /**
399  * \brief Assemble local entity block matrix
400  * @param row_data row data (consist base functions on row entity)
401  * @param col_data column data (consist base functions on column
402  * entity)
403  * @return error code
404  */
405  MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data);
406 
407  /**
408  * \brief Assemble local entity right-hand vector
409  * @param row_data row data (consist base functions on row entity)
410  * @param col_data column data (consist base functions on column
411  * entity)
412  * @return error code
413  */
414  MoFEMErrorCode aSsemble(EntData &row_data);
415  };
416 
417  struct OpRhs_dx : public OpAssemble {
418 
419  OpRhs_dx(const std::string row_field, const std::string col_field,
420  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
421 
422  protected:
423  MoFEMErrorCode iNtegrate(EntData &row_data);
424  };
425 
426  template <int S = 0> struct OpLhs_dx_dx : public OpAssemble {
427 
428  OpLhs_dx_dx(const std::string row_field, const std::string col_field,
429  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
430 
431  protected:
432  /**
433  * \brief Integrate B^T D B operator
434  * @param row_data row data (consist base functions on row entity)
435  * @param col_data column data (consist base functions on column entity)
436  * @return error code
437  */
438  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
439  };
440 
441  struct OpAleRhs_dx : public OpAssemble {
442 
443  OpAleRhs_dx(const std::string row_field, const std::string col_field,
444  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
445 
446  protected:
447  MoFEMErrorCode iNtegrate(EntData &row_data);
448  };
449 
450  template <int S = 0> struct OpAleLhs_dx_dx : public OpAssemble {
451 
452  OpAleLhs_dx_dx(const std::string row_field, const std::string col_field,
453  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
454 
455  protected:
456  /**
457  * \brief Integrate B^T D B operator
458  * @param row_data row data (consist base functions on row entity)
459  * @param col_data column data (consist base functions on column entity)
460  * @return error code
461  */
462  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
463  };
464 
465  template <int S = 0> struct OpAleLhs_dx_dX : public OpAssemble {
466 
467  OpAleLhs_dx_dX(const std::string row_field, const std::string col_field,
468  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
469 
470  protected:
471  /**
472  * \brief Integrate tangent stiffness for spatial momentum
473  * @param row_data row data (consist base functions on row entity)
474  * @param col_data column data (consist base functions on column entity)
475  * @return error code
476  */
477  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
478  };
479 
480  struct OpAleLhsWithDensity_dx_dX : public OpAssemble {
481 
482  boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
483  boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
484  const double rhoN;
485  const double rHo0;
486 
487  OpAleLhsWithDensity_dx_dX(
488  const std::string row_field, const std::string col_field,
489  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
490  boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
491  boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
492  const double rho_n, const double rho_0);
493 
494  protected:
495  /**
496  * \brief Integrate tangent stiffness for spatial momentum
497  * @param row_data row data (consist base functions on row entity)
498  * @param col_data column data (consist base functions on column entity)
499  * @return error code
500  */
501  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
502  };
503 
504  struct OpAleLhsWithDensity_dX_dX : public OpAssemble {
505 
506  boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
507  boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
508  const double rhoN;
509  const double rHo0;
510 
511  OpAleLhsWithDensity_dX_dX(
512  const std::string row_field, const std::string col_field,
513  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
514  boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
515  boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
516  const double rho_n, const double rho_0);
517 
518  protected:
519  /**
520  * \brief Integrate tangent stiffness for material momentum
521  * @param row_data row data (consist base functions on row entity)
522  * @param col_data column data (consist base functions on column entity)
523  * @return error code
524  */
525  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
526  };
527 
528  struct OpAleRhs_dX : public OpAssemble {
529 
530  OpAleRhs_dX(const std::string row_field, const std::string col_field,
531  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
532 
533  protected:
534  MoFEMErrorCode iNtegrate(EntData &row_data);
535  };
536 
537  template <int S = 0> struct OpAleLhs_dX_dX : public OpAssemble {
538 
539  OpAleLhs_dX_dX(const std::string row_field, const std::string col_field,
540  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
541 
542  protected:
543  /**
544  * \brief Integrate tangent stiffness for material momentum
545  * @param row_data row data (consist base functions on row entity)
546  * @param col_data column data (consist base functions on column entity)
547  * @return error code
548  */
549  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
550  };
551 
552  template <int S = 0> struct OpAleLhsPre_dX_dx : public VolUserDataOperator {
553 
554  OpAleLhsPre_dX_dx(const std::string row_field, const std::string col_field,
555  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
556 
557  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
558 
559  private:
560  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
561  };
562 
563  struct OpAleLhs_dX_dx : public OpAssemble {
564 
565  OpAleLhs_dX_dx(const std::string row_field, const std::string col_field,
566  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
567  : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
568 
569  protected:
570  /**
571  * \brief Integrate tangent stiffness for material momentum
572  * @param row_data row data (consist base functions on row entity)
573  * @param col_data column data (consist base functions on column entity)
574  * @return error code
575  */
576  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
577  };
578 
579  template <class ELEMENT>
580  struct OpPostProcHookeElement : public ELEMENT::UserDataOperator {
581  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
582  map<int, BlockData>
583  &blockSetsPtr; // FIXME: (works only with the first block)
584  moab::Interface &postProcMesh;
585  std::vector<EntityHandle> &mapGaussPts;
586  bool isALE;
587 
588  OpPostProcHookeElement(const string row_field,
589  boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
590  map<int, BlockData> &block_sets_ptr,
591  moab::Interface &post_proc_mesh,
592  std::vector<EntityHandle> &map_gauss_pts,
593  bool is_ale = false);
594 
595  MoFEMErrorCode doWork(int side, EntityType type,
596  DataForcesAndSourcesCore::EntData &data);
597  };
598 
599  static MoFEMErrorCode
600  setBlocks(MoFEM::Interface &m_field,
601  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr);
602 
603  static MoFEMErrorCode
604  addElasticElement(MoFEM::Interface &m_field,
605  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
606  const std::string element_name, const std::string x_field,
607  const std::string X_field, const bool ale);
608 
609  static MoFEMErrorCode
610  setOperators(boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr,
611  boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr,
612  boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
613  const std::string x_field, const std::string X_field,
614  const bool ale, const bool field_disp,
615  const EntityType type = MBTET);
616 
617  static MoFEMErrorCode
618  calculateEnergy(DM dm, boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
619  const std::string x_field, const std::string X_field,
620  const bool ale, const bool field_disp, Vec *v_energy_ptr);
621 
622 private:
623  MatrixDouble invJac;
624 };
625 
626 template <bool D>
627 HookeElement::OpCalculateStrain<D>::OpCalculateStrain(
628  const std::string row_field, const std::string col_field,
629  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
630  : VolUserDataOperator(row_field, col_field, OPROW, true),
631  dataAtPts(data_at_pts) {
632  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
633 }
634 
635 template <bool D>
636 MoFEMErrorCode HookeElement::OpCalculateStrain<D>::doWork(int row_side,
637  EntityType row_type,
638  EntData &row_data) {
639  MoFEMFunctionBegin;
640  FTensor::Index<'i', 3> i;
641  FTensor::Index<'j', 3> j;
642  // get number of integration points
643  const int nb_integration_pts = getGaussPts().size2();
644  dataAtPts->smallStrainMat->resize(6, nb_integration_pts, false);
645  auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
646  auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hMat));
647 
648  for (int gg = 0; gg != nb_integration_pts; ++gg) {
649  t_strain(i, j) = (t_h(i, j) || t_h(j, i)) / 2.;
650 
651  // If displacement field, not field o spatial positons is given
652  if (!D) {
653  t_strain(0, 0) -= 1;
654  t_strain(1, 1) -= 1;
655  t_strain(2, 2) -= 1;
656  }
657 
658  ++t_strain;
659  ++t_h;
660  }
661  MoFEMFunctionReturn(0);
662 }
663 
664 template <int S>
665 HookeElement::OpAleLhs_dx_dx<S>::OpAleLhs_dx_dx(
666  const std::string row_field, const std::string col_field,
667  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
668  : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
669 
670 template <int S>
671 MoFEMErrorCode HookeElement::OpAleLhs_dx_dx<S>::iNtegrate(EntData &row_data,
672  EntData &col_data) {
673  MoFEMFunctionBegin;
674 
675  // get sub-block (3x3) of local stiffens matrix, here represented by
676  // second order tensor
677  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
678  return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
679  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
680  &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
681  &m(r + 2, c + 2));
682  };
683 
684  FTensor::Index<'i', 3> i;
685  FTensor::Index<'j', 3> j;
686  FTensor::Index<'k', 3> k;
687  FTensor::Index<'l', 3> l;
688 
689  // get element volume
690  double vol = getVolume();
691 
692  // get intergrayion weights
693  auto t_w = getFTensor0IntegrationWeight();
694 
695  // get derivatives of base functions on rows
696  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
697  const int row_nb_base_fun = row_data.getN().size2();
698 
699  // Elastic stiffness tensor (4th rank tensor with minor and major
700  // symmetry)
701  FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
702  MAT_TO_DDG(dataAtPts->stiffnessMat));
703 
704  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
705  auto &det_H = *dataAtPts->detHVec;
706 
707  // iterate over integration points
708  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
709 
710  // calculate scalar weight times element volume
711  double a = t_w * vol * det_H[gg];
712 
713  // iterate over row base functions
714  int rr = 0;
715  for (; rr != nbRows / 3; ++rr) {
716 
717  // get sub matrix for the row
718  auto t_m = get_tensor2(K, 3 * rr, 0);
719 
720  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
721  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
722 
723  FTensor::Christof<double, 3, 3> t_rowD;
724  // I mix up the indices here so that it behaves like a
725  // Dg. That way I don't have to have a separate wrapper
726  // class Christof_Expr, which simplifies things.
727  t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base_pulled(i));
728 
729  // get derivatives of base functions for columns
730  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
731 
732  // iterate column base functions
733  for (int cc = 0; cc != nbCols / 3; ++cc) {
734 
735  FTensor::Tensor1<double, 3> t_col_diff_base_pulled;
736  t_col_diff_base_pulled(j) = t_col_diff_base(i) * t_invH(i, j);
737 
738  // integrate block local stiffens matrix
739  t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base_pulled(k);
740 
741  // move to next column base function
742  ++t_col_diff_base;
743 
744  // move to next block of local stiffens matrix
745  ++t_m;
746  }
747 
748  // move to next row base function
749  ++t_row_diff_base;
750  }
751 
752  for (; rr != row_nb_base_fun; ++rr)
753  ++t_row_diff_base;
754 
755  // move to next integration weight
756  ++t_w;
757  ++t_D;
758  ++t_invH;
759  }
760 
761  MoFEMFunctionReturn(0);
762 }
763 
764 template <int S>
765 HookeElement::OpCalculateHomogeneousStiffness<S>::
766  OpCalculateHomogeneousStiffness(
767  const std::string row_field, const std::string col_field,
768  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
769  boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
770  : VolUserDataOperator(row_field, col_field, OPROW, true),
771  blockSetsPtr(block_sets_ptr), dataAtPts(data_at_pts) {
772  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
773 }
774 
775 template <int S>
776 MoFEMErrorCode HookeElement::OpCalculateHomogeneousStiffness<S>::doWork(
777  int row_side, EntityType row_type, EntData &row_data) {
778  MoFEMFunctionBegin;
779 
780  for (auto &m : (*blockSetsPtr)) {
781  if (m.second.tEts.find(getFEEntityHandle()) != m.second.tEts.end()) {
782 
783  dataAtPts->stiffnessMat->resize(36, 1, false);
784  FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
785  MAT_TO_DDG(dataAtPts->stiffnessMat));
786  const double young = m.second.E;
787  const double poisson = m.second.PoissonRatio;
788 
789  // coefficient used in intermediate calculation
790  const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
791 
792  FTensor::Index<'i', 3> i;
793  FTensor::Index<'j', 3> j;
794  FTensor::Index<'k', 3> k;
795  FTensor::Index<'l', 3> l;
796 
797  t_D(i, j, k, l) = 0.;
798 
799  t_D(0, 0, 0, 0) = 1 - poisson;
800  t_D(1, 1, 1, 1) = 1 - poisson;
801  t_D(2, 2, 2, 2) = 1 - poisson;
802 
803  t_D(0, 1, 0, 1) = 0.5 * (1 - 2 * poisson);
804  t_D(0, 2, 0, 2) = 0.5 * (1 - 2 * poisson);
805  t_D(1, 2, 1, 2) = 0.5 * (1 - 2 * poisson);
806 
807  t_D(0, 0, 1, 1) = poisson;
808  t_D(1, 1, 0, 0) = poisson;
809  t_D(0, 0, 2, 2) = poisson;
810  t_D(2, 2, 0, 0) = poisson;
811  t_D(1, 1, 2, 2) = poisson;
812  t_D(2, 2, 1, 1) = poisson;
813 
814  t_D(i, j, k, l) *= coefficient;
815 
816  break;
817  }
818  }
819 
820  MoFEMFunctionReturn(0);
821 }
822 
823 template <int S>
824 HookeElement::OpCalculateStress<S>::OpCalculateStress(
825  const std::string row_field, const std::string col_field,
826  boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
827  : VolUserDataOperator(row_field, col_field, OPROW, true),
828  dataAtPts(data_at_pts) {
829  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
830 }
831 
832 template <int S>
833 MoFEMErrorCode HookeElement::OpCalculateStress<S>::doWork(int row_side,
834  EntityType row_type,
835  EntData &row_data) {
836  MoFEMFunctionBegin;
837  // get number of integration points
838  const int nb_integration_pts = getGaussPts().size2();
839  auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
840  dataAtPts->cauchyStressMat->resize(6, nb_integration_pts, false);
841  auto t_cauchy_stress =
842  getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
843 
844  FTensor::Index<'i', 3> i;
845  FTensor::Index<'j', 3> j;
846  FTensor::Index<'k', 3> k;
847  FTensor::Index<'l', 3> l;
848 
849  // elastic stiffness tensor (4th rank tensor with minor and major
850  // symmetry)
851  FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
852  MAT_TO_DDG(dataAtPts->stiffnessMat));
853  for (int gg = 0; gg != nb_integration_pts; ++gg) {
854  t_cauchy_stress(i, j) = t_D(i, j, k, l) * t_strain(k, l);
855  ++t_strain;
856  ++t_cauchy_stress;
857  ++t_D;
858  }
859  MoFEMFunctionReturn(0);
860 }
861 
862 template <int S>
863 HookeElement::OpLhs_dx_dx<S>::OpLhs_dx_dx(
864  const std::string row_field, const std::string col_field,
865  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
866  : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
867 
868 template <int S>
869 MoFEMErrorCode HookeElement::OpLhs_dx_dx<S>::iNtegrate(EntData &row_data,
870  EntData &col_data) {
871  MoFEMFunctionBegin;
872 
873  // get sub-block (3x3) of local stiffens matrix, here represented by
874  // second order tensor
875  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
876  return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
877  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
878  &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
879  &m(r + 2, c + 2));
880  };
881 
882  FTensor::Index<'i', 3> i;
883  FTensor::Index<'j', 3> j;
884  FTensor::Index<'k', 3> k;
885  FTensor::Index<'l', 3> l;
886 
887  // get element volume
888  double vol = getVolume();
889 
890  // get intergrayion weights
891  auto t_w = getFTensor0IntegrationWeight();
892 
893  // get derivatives of base functions on rows
894  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
895  const int row_nb_base_fun = row_data.getN().size2();
896 
897  // Elastic stiffness tensor (4th rank tensor with minor and major
898  // symmetry)
899  FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
900  MAT_TO_DDG(dataAtPts->stiffnessMat));
901 
902  // iterate over integration points
903  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
904 
905  // calculate scalar weight times element volume
906  double a = t_w * vol;
907  if (getHoGaussPtsDetJac().size()) {
908  a *= getHoGaussPtsDetJac()[gg];
909  }
910 
911  // iterate over row base functions
912  int rr = 0;
913  for (; rr != nbRows / 3; ++rr) {
914 
915  // get sub matrix for the row
916  auto t_m = get_tensor2(K, 3 * rr, 0);
917 
918  // get derivatives of base functions for columns
919  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
920 
921  FTensor::Christof<double, 3, 3> t_rowD;
922  // I mix up the indices here so that it behaves like a
923  // Dg. That way I don't have to have a separate wrapper
924  // class Christof_Expr, which simplifies things.
925  t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base(i));
926 
927  // iterate column base functions
928  for (int cc = 0; cc != nbCols / 3; ++cc) {
929 
930  // integrate block local stiffens matrix
931  t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base(k);
932 
933  // move to next column base function
934  ++t_col_diff_base;
935 
936  // move to next block of local stiffens matrix
937  ++t_m;
938  }
939 
940  // move to next row base function
941  ++t_row_diff_base;
942  }
943 
944  for (; rr != row_nb_base_fun; ++rr)
945  ++t_row_diff_base;
946 
947  // move to next integration weight
948  ++t_w;
949  ++t_D;
950  }
951 
952  MoFEMFunctionReturn(0);
953 }
954 
955 template <int S>
956 HookeElement::OpAleLhs_dx_dX<S>::OpAleLhs_dx_dX(
957  const std::string row_field, const std::string col_field,
958  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
959  : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
960 
961 template <int S>
962 MoFEMErrorCode HookeElement::OpAleLhs_dx_dX<S>::iNtegrate(EntData &row_data,
963  EntData &col_data) {
964  MoFEMFunctionBegin;
965 
966  // get sub-block (3x3) of local stiffens matrix, here represented by
967  // second order tensor
968  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
969  return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
970  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
971  &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
972  &m(r + 2, c + 2));
973  };
974 
975  FTensor::Index<'i', 3> i;
976  FTensor::Index<'j', 3> j;
977  FTensor::Index<'k', 3> k;
978  FTensor::Index<'l', 3> l;
979  FTensor::Index<'m', 3> m;
980  FTensor::Index<'n', 3> n;
981 
982  // get element volume
983  double vol = getVolume();
984 
985  // get intergrayion weights
986  auto t_w = getFTensor0IntegrationWeight();
987 
988  // get derivatives of base functions on rows
989  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
990  const int row_nb_base_fun = row_data.getN().size2();
991 
992  // Elastic stiffness tensor (4th rank tensor with minor and major
993  // symmetry)
994  FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
995  MAT_TO_DDG(dataAtPts->stiffnessMat));
996 
997  auto t_cauchy_stress =
998  getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
999  auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1000  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1001  auto &det_H = *dataAtPts->detHVec;
1002 
1003  // iterate over integration points
1004  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1005 
1006  // calculate scalar weight times element volume
1007  double a = t_w * vol * det_H[gg];
1008 
1009  FTensor::Tensor4<double, 3, 3, 3, 3> t_F_dX;
1010  t_F_dX(i, j, k, l) = -(t_h(i, m) * t_invH(m, k)) * t_invH(l, j);
1011 
1012  // iterate over row base functions
1013  int rr = 0;
1014  for (; rr != nbRows / 3; ++rr) {
1015 
1016  // get sub matrix for the row
1017  auto t_m = get_tensor2(K, 3 * rr, 0);
1018 
1019  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1020  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1021 
1022  FTensor::Tensor1<double, 3> t_row_stress;
1023  t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_cauchy_stress(i, j);
1024 
1025  FTensor::Tensor3<double, 3, 3, 3> t_row_diff_base_pulled_dX;
1026  t_row_diff_base_pulled_dX(j, k, l) =
1027  -(t_invH(i, k) * t_row_diff_base(i)) * t_invH(l, j);
1028 
1029  FTensor::Tensor3<double, 3, 3, 3> t_row_dX_stress;
1030  t_row_dX_stress(i, k, l) =
1031  a * (t_row_diff_base_pulled_dX(j, k, l) * t_cauchy_stress(j, i));
1032 
1033  FTensor::Christof<double, 3, 3> t_row_D;
1034  t_row_D(l, j, k) = (a * t_row_diff_base_pulled(i)) * t_D(i, j, k, l);
1035 
1036  FTensor::Tensor3<double, 3, 3, 3> t_row_stress_dX;
1037  // FIXME: This operator is not implemented, doing operation by hand
1038  // t_row_stress_dX(i, m, n) = t_row_D(i, k, l) * t_F_dX(k, l, m, n);
1039  t_row_stress_dX(i, j, k) = 0;
1040  for (int ii = 0; ii != 3; ++ii)
1041  for (int mm = 0; mm != 3; ++mm)
1042  for (int nn = 0; nn != 3; ++nn) {
1043  auto &v = t_row_stress_dX(ii, mm, nn);
1044  for (int kk = 0; kk != 3; ++kk)
1045  for (int ll = 0; ll != 3; ++ll)
1046  v += t_row_D(ii, kk, ll) * t_F_dX(kk, ll, mm, nn);
1047  }
1048 
1049  // get derivatives of base functions for columns
1050  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1051 
1052  // iterate column base functions
1053  for (int cc = 0; cc != nbCols / 3; ++cc) {
1054 
1055  t_m(i, k) += t_row_stress(i) * (t_invH(j, k) * t_col_diff_base(j));
1056  t_m(i, k) += t_row_dX_stress(i, k, l) * t_col_diff_base(l);
1057  t_m(i, k) += t_row_stress_dX(i, k, l) * t_col_diff_base(l);
1058 
1059  // move to next column base function
1060  ++t_col_diff_base;
1061 
1062  // move to next block of local stiffens matrix
1063  ++t_m;
1064  }
1065 
1066  // move to next row base function
1067  ++t_row_diff_base;
1068  }
1069 
1070  for (; rr != row_nb_base_fun; ++rr)
1071  ++t_row_diff_base;
1072 
1073  // move to next integration weight
1074  ++t_w;
1075  ++t_D;
1076  ++t_cauchy_stress;
1077  ++t_invH;
1078  ++t_h;
1079  }
1080 
1081  MoFEMFunctionReturn(0);
1082 }
1083 
1084 template <int S>
1085 HookeElement::OpAleLhs_dX_dX<S>::OpAleLhs_dX_dX(
1086  const std::string row_field, const std::string col_field,
1087  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1088  : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
1089 
1090 template <int S>
1091 MoFEMErrorCode HookeElement::OpAleLhs_dX_dX<S>::iNtegrate(EntData &row_data,
1092  EntData &col_data) {
1093  MoFEMFunctionBegin;
1094 
1095  // get sub-block (3x3) of local stiffens matrix, here represented by
1096  // second order tensor
1097  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1098  return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
1099  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1100  &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1101  &m(r + 2, c + 2));
1102  };
1103 
1104  FTensor::Index<'i', 3> i;
1105  FTensor::Index<'j', 3> j;
1106  FTensor::Index<'k', 3> k;
1107  FTensor::Index<'l', 3> l;
1108  FTensor::Index<'m', 3> m;
1109  FTensor::Index<'n', 3> n;
1110 
1111  // get element volume
1112  double vol = getVolume();
1113 
1114  // get intergrayion weights
1115  auto t_w = getFTensor0IntegrationWeight();
1116 
1117  // get derivatives of base functions on rows
1118  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1119  const int row_nb_base_fun = row_data.getN().size2();
1120 
1121  // Elastic stiffness tensor (4th rank tensor with minor and major
1122  // symmetry)
1123  FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
1124  MAT_TO_DDG(dataAtPts->stiffnessMat));
1125  auto t_cauchy_stress =
1126  getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1127  auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
1128  auto t_eshelby_stress =
1129  getFTensor2FromMat<3, 3>(*dataAtPts->eshelbyStressMat);
1130  auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1131  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1132  auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
1133  auto &det_H = *dataAtPts->detHVec;
1134 
1135  // iterate over integration points
1136  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1137 
1138  // calculate scalar weight times element volume
1139  double a = t_w * vol * det_H[gg];
1140 
1141  FTensor::Tensor4<double, 3, 3, 3, 3> t_F_dX;
1142  t_F_dX(i, j, k, l) = -(t_h(i, m) * t_invH(m, k)) * t_invH(l, j);
1143 
1144  FTensor::Tensor4<double, 3, 3, 3, 3> t_D_strain_dX;
1145  t_D_strain_dX(i, j, m, n) = 0.;
1146  for (int ii = 0; ii != 3; ++ii)
1147  for (int jj = 0; jj != 3; ++jj)
1148  for (int ll = 0; ll != 3; ++ll)
1149  for (int kk = 0; kk != 3; ++kk) {
1150  auto &v = t_D_strain_dX(ii, jj, kk, ll);
1151  for (int mm = 0; mm != 3; ++mm)
1152  for (int nn = 0; nn != 3; ++nn)
1153  v += t_D(ii, jj, mm, nn) * t_F_dX(mm, nn, kk, ll);
1154  }
1155 
1156  FTensor::Tensor4<double, 3, 3, 3, 3> t_eshelby_stress_dX;
1157  t_eshelby_stress_dX(i, j, m, n) = t_F(k, i) * t_D_strain_dX(k, j, m, n);
1158 
1159  for (int ii = 0; ii != 3; ++ii)
1160  for (int jj = 0; jj != 3; ++jj)
1161  for (int mm = 0; mm != 3; ++mm)
1162  for (int nn = 0; nn != 3; ++nn) {
1163  auto &v = t_eshelby_stress_dX(ii, jj, mm, nn);
1164  for (int kk = 0; kk != 3; ++kk)
1165  v += t_F_dX(kk, ii, mm, nn) * t_cauchy_stress(kk, jj);
1166  }
1167 
1168  t_eshelby_stress_dX(i, j, k, l) *= -1;
1169 
1170  FTensor::Tensor2<double, 3, 3> t_energy_dX;
1171  t_energy_dX(k, l) = t_F_dX(i, j, k, l) * t_cauchy_stress(i, j);
1172  t_energy_dX(k, l) +=
1173  (t_strain(m, n) * t_D(m, n, i, j)) * t_F_dX(i, j, k, l);
1174  t_energy_dX(k, l) /= 2.;
1175 
1176  for (int kk = 0; kk != 3; ++kk)
1177  for (int ll = 0; ll != 3; ++ll) {
1178  auto v = t_energy_dX(kk, ll);
1179  for (int ii = 0; ii != 3; ++ii)
1180  t_eshelby_stress_dX(ii, ii, kk, ll) += v;
1181  }
1182 
1183  // iterate over row base functions
1184  int rr = 0;
1185  for (; rr != nbRows / 3; ++rr) {
1186 
1187  // get sub matrix for the row
1188  auto t_m = get_tensor2(K, 3 * rr, 0);
1189 
1190  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1191  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1192 
1193  FTensor::Tensor1<double, 3> t_row_stress;
1194  t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
1195 
1196  FTensor::Tensor3<double, 3, 3, 3> t_row_diff_base_pulled_dX;
1197  t_row_diff_base_pulled_dX(j, k, l) =
1198  -(t_row_diff_base(i) * t_invH(i, k)) * t_invH(l, j);
1199 
1200  FTensor::Tensor3<double, 3, 3, 3> t_row_dX_stress;
1201  t_row_dX_stress(i, k, l) =
1202  a * (t_row_diff_base_pulled_dX(j, k, l) * t_eshelby_stress(i, j));
1203 
1204  FTensor::Tensor3<double, 3, 3, 3> t_row_stress_dX;
1205  t_row_stress_dX(i, m, n) =
1206  a * t_row_diff_base_pulled(j) * t_eshelby_stress_dX(i, j, m, n);
1207 
1208  // get derivatives of base functions for columns
1209  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1210 
1211  // iterate column base functions
1212  for (int cc = 0; cc != nbCols / 3; ++cc) {
1213 
1214  t_m(i, k) += t_row_stress(i) * (t_invH(j, k) * t_col_diff_base(j));
1215  t_m(i, k) += t_row_dX_stress(i, k, l) * t_col_diff_base(l);
1216  t_m(i, k) += t_row_stress_dX(i, k, l) * t_col_diff_base(l);
1217 
1218  // move to next column base function
1219  ++t_col_diff_base;
1220 
1221  // move to next block of local stiffens matrix
1222  ++t_m;
1223  }
1224 
1225  // move to next row base function
1226  ++t_row_diff_base;
1227  }
1228 
1229  for (; rr != row_nb_base_fun; ++rr)
1230  ++t_row_diff_base;
1231 
1232  // move to next integration weight
1233  ++t_w;
1234  ++t_D;
1235  ++t_cauchy_stress;
1236  ++t_strain;
1237  ++t_eshelby_stress;
1238  ++t_h;
1239  ++t_invH;
1240  ++t_F;
1241  }
1242 
1243  MoFEMFunctionReturn(0);
1244 }
1245 
1246 template <int S>
1247 HookeElement::OpAleLhsPre_dX_dx<S>::OpAleLhsPre_dX_dx(
1248  const std::string row_field, const std::string col_field,
1249  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1250  : VolUserDataOperator(row_field, col_field, OPROW, true),
1251  dataAtPts(data_at_pts) {
1252  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1253 }
1254 
1255 template <int S>
1256 MoFEMErrorCode HookeElement::OpAleLhsPre_dX_dx<S>::doWork(int row_side,
1257  EntityType row_type,
1258  EntData &row_data) {
1259  MoFEMFunctionBegin;
1260 
1261  const int nb_integration_pts = row_data.getN().size1();
1262 
1263  auto get_eshelby_stress_dx = [this, nb_integration_pts]() {
1264  FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>
1265  t_eshelby_stress_dx;
1266  dataAtPts->eshelbyStress_dx->resize(81, nb_integration_pts, false);
1267  int mm = 0;
1268  for (int ii = 0; ii != 3; ++ii)
1269  for (int jj = 0; jj != 3; ++jj)
1270  for (int kk = 0; kk != 3; ++kk)
1271  for (int ll = 0; ll != 3; ++ll)
1272  t_eshelby_stress_dx.ptr(ii, jj, kk, ll) =
1273  &(*dataAtPts->eshelbyStress_dx)(mm++, 0);
1274  return t_eshelby_stress_dx;
1275  };
1276 
1277  auto t_eshelby_stress_dx = get_eshelby_stress_dx();
1278 
1279  FTensor::Index<'i', 3> i;
1280  FTensor::Index<'j', 3> j;
1281  FTensor::Index<'k', 3> k;
1282  FTensor::Index<'l', 3> l;
1283  FTensor::Index<'m', 3> m;
1284  FTensor::Index<'n', 3> n;
1285 
1286  // Elastic stiffness tensor (4th rank tensor with minor and major
1287  // symmetry)
1288  FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
1289  MAT_TO_DDG(dataAtPts->stiffnessMat));
1290  auto t_cauchy_stress =
1291  getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1292  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1293  auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
1294 
1295  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1296 
1297  t_eshelby_stress_dx(i, j, m, n) =
1298  (t_F(k, i) * t_D(k, j, m, l)) * t_invH(n, l);
1299  for (int ii = 0; ii != 3; ++ii)
1300  for (int jj = 0; jj != 3; ++jj)
1301  for (int mm = 0; mm != 3; ++mm)
1302  for (int nn = 0; nn != 3; ++nn) {
1303  auto &v = t_eshelby_stress_dx(ii, jj, mm, nn);
1304  v += t_invH(nn, ii) * t_cauchy_stress(mm, jj);
1305  }
1306  t_eshelby_stress_dx(i, j, k, l) *= -1;
1307 
1308  FTensor::Tensor2<double, 3, 3> t_energy_dx;
1309  t_energy_dx(m, n) = t_invH(n, j) * t_cauchy_stress(m, j);
1310 
1311  for (int mm = 0; mm != 3; ++mm)
1312  for (int nn = 0; nn != 3; ++nn) {
1313  auto v = t_energy_dx(mm, nn);
1314  for (int ii = 0; ii != 3; ++ii)
1315  t_eshelby_stress_dx(ii, ii, mm, nn) += v;
1316  }
1317 
1318  ++t_D;
1319  ++t_invH;
1320  ++t_cauchy_stress;
1321  ++t_eshelby_stress_dx;
1322  ++t_F;
1323  }
1324 
1325  MoFEMFunctionReturn(0);
1326 }
1327 
1328 template <class ELEMENT>
1329 HookeElement::OpPostProcHookeElement<ELEMENT>::OpPostProcHookeElement(
1330  const string row_field, boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1331  map<int, BlockData> &block_sets_ptr, moab::Interface &post_proc_mesh,
1332  std::vector<EntityHandle> &map_gauss_pts, bool is_ale)
1333  : ELEMENT::UserDataOperator(row_field, UserDataOperator::OPROW),
1334  dataAtPts(data_at_pts), blockSetsPtr(block_sets_ptr),
1335  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts), isALE(is_ale) {}
1336 
1337 template <class ELEMENT>
1338 MoFEMErrorCode HookeElement::OpPostProcHookeElement<ELEMENT>::doWork(
1339  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1340  MoFEMFunctionBegin;
1341 
1342  if (type != MBVERTEX) {
1343  MoFEMFunctionReturnHot(0);
1344  }
1345 
1346  auto tensor_to_tensor = [](const auto &t1, auto &t2) {
1347  t2(0, 0) = t1(0, 0);
1348  t2(1, 1) = t1(1, 1);
1349  t2(2, 2) = t1(2, 2);
1350  t2(0, 1) = t2(1, 0) = t1(1, 0);
1351  t2(0, 2) = t2(2, 0) = t1(2, 0);
1352  t2(1, 2) = t2(2, 1) = t1(2, 1);
1353  };
1354 
1355  double def_VAL[9];
1356  bzero(def_VAL, 9 * sizeof(double));
1357 
1358  Tag th_stress;
1359  CHKERR postProcMesh.tag_get_handle("STRESS", 9, MB_TYPE_DOUBLE, th_stress,
1360  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
1361  Tag th_psi;
1362  CHKERR postProcMesh.tag_get_handle("ENERGY", 1, MB_TYPE_DOUBLE, th_psi,
1363  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
1364 
1365  const int nb_integration_pts = mapGaussPts.size();
1366 
1367  FTensor::Index<'i', 3> i;
1368  FTensor::Index<'j', 3> j;
1369  FTensor::Index<'k', 3> k;
1370  FTensor::Index<'l', 3> l;
1371 
1372  auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hMat));
1373  auto t_H = getFTensor2FromMat<3, 3>(*dataAtPts->HMat);
1374 
1375  dataAtPts->stiffnessMat->resize(36, 1, false);
1376  FTensor::Ddg<FTensor::PackPtr<double *, 1>, 3, 3> t_D(
1377  MAT_TO_DDG(dataAtPts->stiffnessMat));
1378  for (auto &m : (blockSetsPtr)) {
1379  const double young = m.second.E;
1380  const double poisson = m.second.PoissonRatio;
1381 
1382  const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
1383 
1384  t_D(i, j, k, l) = 0.;
1385  t_D(0, 0, 0, 0) = t_D(1, 1, 1, 1) = t_D(2, 2, 2, 2) = 1 - poisson;
1386  t_D(0, 1, 0, 1) = t_D(0, 2, 0, 2) = t_D(1, 2, 1, 2) =
1387  0.5 * (1 - 2 * poisson);
1388  t_D(0, 0, 1, 1) = t_D(1, 1, 0, 0) = t_D(0, 0, 2, 2) = t_D(2, 2, 0, 0) =
1389  t_D(1, 1, 2, 2) = t_D(2, 2, 1, 1) = poisson;
1390  t_D(i, j, k, l) *= coefficient;
1391 
1392  break; // FIXME: calculates only first block
1393  }
1394 
1395  double detH = 0.;
1396  FTensor::Tensor2<double, 3, 3> t_invH;
1397  FTensor::Tensor2<double, 3, 3> t_F;
1398  FTensor::Tensor2<double, 3, 3> t_stress;
1399  FTensor::Tensor2<double, 3, 3> t_small_strain;
1400  FTensor::Tensor2_symmetric<double, 3> t_stress_symm;
1401  FTensor::Tensor2_symmetric<double, 3> t_small_strain_symm;
1402 
1403  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1404 
1405  if (!isALE) {
1406  t_small_strain_symm(i, j) = (t_h(i, j) || t_h(j, i)) / 2.;
1407  } else {
1408  CHKERR determinantTensor3by3(t_H, detH);
1409  CHKERR invertTensor3by3(t_H, detH, t_invH);
1410  t_F(i, j) = t_h(i, k) * t_invH(k, j);
1411  t_small_strain_symm(i, j) = (t_F(i, j) || t_F(j, i)) / 2.;
1412 
1413  t_small_strain_symm(0, 0) -= 1;
1414  t_small_strain_symm(1, 1) -= 1;
1415  t_small_strain_symm(2, 2) -= 1;
1416  ++t_H;
1417  }
1418 
1419  // symmetric tensors need improvement
1420  t_stress_symm(i, j) = t_D(i, j, k, l) * t_small_strain_symm(k, l);
1421  tensor_to_tensor(t_stress_symm, t_stress);
1422 
1423  const double psi = 0.5 * t_stress_symm(i, j) * t_small_strain_symm(i, j);
1424 
1425  CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &psi);
1426  CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
1427  &t_stress(0, 0));
1428 
1429  ++t_h;
1430  }
1431 
1432  MoFEMFunctionReturn(0);
1433 }
1434 
1435 #endif // __HOOKE_ELEMENT_HPP
OpCalculateStiffnessScaledByDensityField::rhoAtGaussPtsPtr
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
Definition: HookeElement.hpp:341
aSsemble
MoFEMErrorCode aSsemble(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Assemble local entity block matrix.
OpAleLhsWithDensity_dx_dX
Definition: HookeElement.hpp:480
OpRhs_dx
Definition: HookeElement.hpp:417
OpPostProcHookeElement::isALE
bool isALE
Definition: HookeElement.hpp:586
OpAleLhsWithDensity_dX_dX::rhoGradAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
Definition: HookeElement.hpp:507
OpAleLhsWithDensity_dx_dX::rhoAtGaussPtsPtr
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
Definition: HookeElement.hpp:482
iNtegrate
MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Integrate B^T D B operator.
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:26
DataAtIntegrationPts::smallStrainMat
boost::shared_ptr< MatrixDouble > smallStrainMat
Definition: HookeElement.hpp:94
OpCalculateStrainAle::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:153
setOperators
static MoFEMErrorCode setOperators(boost::shared_ptr< ForcesAndSourcesCore > fe_lhs_ptr, boost::shared_ptr< ForcesAndSourcesCore > fe_rhs_ptr, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, const EntityType type=MBTET)
OpAleLhsWithDensity_dX_dX::rHo0
const double rHo0
Definition: HookeElement.hpp:509
MoFEM::determinantTensor3by3
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:482
OpLhs_dx_dx
Definition: HookeElement.hpp:426
OpCalculateStiffnessScaledByDensityField::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:339
DataAtIntegrationPts::eshelbyStress_dx
boost::shared_ptr< MatrixDouble > eshelbyStress_dx
Definition: HookeElement.hpp:107
EntData
DataForcesAndSourcesCore::EntData EntData
Definition: HookeElement.hpp:87
OpAssemble::nbCols
int nbCols
number if dof on column
Definition: HookeElement.hpp:390
OpAleLhs_dx_dX
Definition: HookeElement.hpp:465
OpCalculateStiffnessScaledByDensityField::blockSetsPtr
boost::shared_ptr< map< int, BlockData > > blockSetsPtr
Definition: HookeElement.hpp:337
MoFEM::Types::MatrixDouble
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
OpPlasticTools::n
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
DataAtIntegrationPts::hMat
boost::shared_ptr< MatrixDouble > hMat
Definition: HookeElement.hpp:95
invJac
MatrixDouble invJac
Definition: HookeElement.hpp:623
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
ConvectiveMassElement
structure grouping operators and data used for calculation of mass (convective) element \ nonlinear_e...
Definition: ConvectiveMassElement.hpp:42
OpAssemble::K
MatrixDouble K
Definition: HookeElement.hpp:380
OpCalculateStrain::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:141
MoFEM::DataForcesAndSourcesCore::EntData::getFTensor1DiffN
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: DataStructures.cpp:415
OpAleLhsWithDensity_dx_dX::rhoGradAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
Definition: HookeElement.hpp:483
DataAtIntegrationPts::stiffnessMat
boost::shared_ptr< MatrixDouble > stiffnessMat
Definition: HookeElement.hpp:103
OpAleLhs_dX_dX
Definition: HookeElement.hpp:537
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
OpAleRhs_dx
Definition: HookeElement.hpp:441
OpAssemble::isDiag
bool isDiag
true if this block is on diagonal
Definition: HookeElement.hpp:392
OpCalculateStress::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:175
DataAtIntegrationPts::invHMat
boost::shared_ptr< MatrixDouble > invHMat
Definition: HookeElement.hpp:100
OpAleLhsWithDensity_dx_dX::rHo0
const double rHo0
Definition: HookeElement.hpp:485
MoFEM::invertTensor3by3
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:463
OpCalculateStiffnessScaledByDensityField::rhoN
const double rhoN
exponent n in E(p) = E * (p / p_0)^n
Definition: HookeElement.hpp:342
OpPostProcHookeElement::mapGaussPts
std::vector< EntityHandle > & mapGaussPts
Definition: HookeElement.hpp:585
DataAtIntegrationPts::HMat
boost::shared_ptr< MatrixDouble > HMat
Definition: HookeElement.hpp:98
OpAleLhs_dX_dx::OpAleLhs_dX_dx
OpAleLhs_dX_dx(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
Definition: HookeElement.hpp:565
MoFEM::DataForcesAndSourcesCore::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: DataStructures.hpp:1148
NonlinearElasticElement::BlockData::tEts
Range tEts
constrains elements in block set
Definition: HookeElement.hpp:49
OpAssemble::transK
MatrixDouble transK
Definition: HookeElement.hpp:381
DataAtIntegrationPts::FMat
boost::shared_ptr< MatrixDouble > FMat
Definition: HookeElement.hpp:96
OpPostProcHookeElement::blockSetsPtr
map< int, BlockData > & blockSetsPtr
Definition: HookeElement.hpp:583
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:88
OpAssemble::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:384
OpPostProcHookeElement
Definition: HookeElement.hpp:580
OpCalculateMassMatrix::OpCalculateMassMatrix
OpCalculateMassMatrix(const std::string row_field, const std::string col_field, BlockData &data, MassBlockData &mass_data, boost::shared_ptr< DataAtIntegrationPts > &common_data, bool symm=true)
Definition: HookeElement.hpp:238
OpAssemble::nbIntegrationPts
int nbIntegrationPts
number of integration points
Definition: HookeElement.hpp:391
OpCalculateHomogeneousStiffness::blockSetsPtr
boost::shared_ptr< map< int, BlockData > > blockSetsPtr
Definition: HookeElement.hpp:217
OpPlasticTools::m
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
ConvectiveMassElement::BlockData::a0
VectorDouble a0
constant acceleration
Definition: ConvectiveMassElement.hpp:118
OpCalculateEshelbyStress::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:202
DataAtIntegrationPts::eshelbyStressMat
boost::shared_ptr< MatrixDouble > eshelbyStressMat
Definition: HookeElement.hpp:105
DataAtIntegrationPts::detHVec
boost::shared_ptr< VectorDouble > detHVec
Definition: HookeElement.hpp:99
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: DataStructures.hpp:1550
OpAleLhs_dX_dx
Definition: HookeElement.hpp:563
OpAssemble::nbRows
int nbRows
number of dofs on rows
Definition: HookeElement.hpp:389
OpCalculateMassMatrix::commonData
boost::shared_ptr< DataAtIntegrationPts > commonData
Definition: HookeElement.hpp:236
OpAleLhsWithDensity_dX_dX
Definition: HookeElement.hpp:504
BasicFiniteElements.hpp
MAT_TO_DDG
#define MAT_TO_DDG(SM)
Definition: HookeElement.hpp:156
OpAssemble::colIndices
VectorInt colIndices
Definition: HookeElement.hpp:387
OpAleLhsWithDensity_dX_dX::rhoAtGaussPtsPtr
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
Definition: HookeElement.hpp:506
setBlocks
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr)
addElasticElement
static MoFEMErrorCode addElasticElement(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr, const std::string element_name, const std::string x_field, const std::string X_field, const bool ale)
OpCalculateStiffnessScaledByDensityField::rHo0
const double rHo0
p_0 reference density in E(p) = E * (p / p_0)^n
Definition: HookeElement.hpp:343
OpAleLhsWithDensity_dX_dX::rhoN
const double rhoN
Definition: HookeElement.hpp:508
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:602
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
OpAleLhsWithDensity_dx_dX::rhoN
const double rhoN
Definition: HookeElement.hpp:484
OpContactTools::k
FTensor::Index< 'k', 2 > k
Definition: ContactOps.hpp:28
OpPostProcHookeElement::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:581
OpCalculateHomogeneousStiffness::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:219
MoFEM::Types::VectorInt
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:72
OpAssemble::nF
VectorDouble nF
Definition: HookeElement.hpp:382
calculateEnergy
static MoFEMErrorCode calculateEnergy(DM dm, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, Vec *v_energy_ptr)
DataAtIntegrationPts::energyVec
boost::shared_ptr< VectorDouble > energyVec
Definition: HookeElement.hpp:104
OpAleLhs_dx_dx
Definition: HookeElement.hpp:450
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1788
BlockData
struct ConvectiveMassElement BlockData
MoFEM::DataForcesAndSourcesCore::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: DataStructures.hpp:51
OpContactTools::j
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
OpCalculateEnergy::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:189
OpAleLhsPre_dX_dx
Definition: HookeElement.hpp:552
OpAssemble::rowIndices
VectorInt rowIndices
Definition: HookeElement.hpp:386
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
OpCalculateMassMatrix::doWork
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Definition: HookeElement.hpp:247
OpPostProcHookeElement::postProcMesh
moab::Interface & postProcMesh
Definition: HookeElement.hpp:584
MoFEM::Types::VectorDouble
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:73
MoFEM::DataForcesAndSourcesCore::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: DataStructures.hpp:1249
OpAleRhs_dX
Definition: HookeElement.hpp:528
NonlinearElasticElement::BlockData
data for calculation heat conductivity and heat capacity elements
Definition: HookeElement.hpp:45
MoFEM::ForcesAndSourcesCore::UserDataOperator
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
Definition: ForcesAndSourcesCore.hpp:84
ConvectiveMassElement::BlockData
data for calculation inertia forces
Definition: ConvectiveMassElement.hpp:116
OpCalculateMassMatrix::massData
MassBlockData & massData
Definition: HookeElement.hpp:234
OpContactTools::l
FTensor::Index< 'l', 2 > l
Definition: ContactOps.hpp:29
OpAleLhsPre_dX_dx::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:560
DataAtIntegrationPts::cauchyStressMat
boost::shared_ptr< MatrixDouble > cauchyStressMat
Definition: HookeElement.hpp:102
OpCalculateMassMatrix
Assemble mass matrix for elastic element TODO: CHANGE FORMULA * .
Definition: HookeElement.hpp:228
MoFEM::DataForcesAndSourcesCore::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: DataStructures.hpp:1420
OpContactTools::w
double w(const double g, const double t)
Definition: ContactOps.hpp:160
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic elementIn order to as...
Definition: HookeElement.hpp:40
OpContactTools::i
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26
Generated on Fri May 15 2020 09:45:59 for MoFEM by Doxygen 1.8.15 and hosted at University of Glasgow