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 
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
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 
86 
89  using VolUserDataOperator =
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 
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 
130  };
131 
132  template <bool D = false>
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 
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 
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 
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 
194 
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>
207 
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  */
230 
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)
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,
252 
253  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
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)
265  const int col_nb_dofs = col_data.getIndices().size();
266  if (!col_nb_dofs)
268  if (dAta.tEts.find(getFEEntityHandle()) == dAta.tEts.end()) {
270  }
271  if (massData.tEts.find(getFEEntityHandle()) == massData.tEts.end()) {
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 
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 
331  }
332  };
333 
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:
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
383 
384  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
385 
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 
481 
482  boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
483  boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
484  const double rhoN;
485  const double rHo0;
486 
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 
505 
506  boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
507  boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
508  const double rhoN;
509  const double rHo0;
510 
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>
581  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
582  map<int, BlockData>
583  &blockSetsPtr; // FIXME: (works only with the first block)
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,
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
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:
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) {
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  }
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>
672  EntData &col_data) {
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) {
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 
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)
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 
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 
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) {
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);
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 
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 
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) {
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 
848 
849  // elastic stiffness tensor (4th rank tensor with minor and major
850  // symmetry)
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  }
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>
870  EntData &col_data) {
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) {
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 
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)
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 
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 
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>
963  EntData &col_data) {
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) {
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 
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)
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 
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 
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 
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>
1092  EntData &col_data) {
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) {
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 
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)
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 
1142  t_F_dX(i, j, k, l) = -(t_h(i, m) * t_invH(m, k)) * t_invH(l, j);
1143 
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 
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) {
1260 
1261  const int nb_integration_pts = row_data.getN().size1();
1262 
1263  auto get_eshelby_stress_dx = [this, nb_integration_pts]() {
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 
1285 
1286  // Elastic stiffness tensor (4th rank tensor with minor and major
1287  // symmetry)
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 
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) {
1341 
1342  if (type != MBVERTEX) {
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 
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);
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.;
1399  FTensor::Tensor2<double, 3, 3> t_small_strain;
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 
1433 }
1434 
1435 #endif // __HOOKE_ELEMENT_HPP
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
MoFEMErrorCode aSsemble(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Assemble local entity block matrix.
bool isALE
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Integrate B^T D B operator.
Deprecated interface functions.
boost::shared_ptr< MatrixDouble > smallStrainMat
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
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)
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:482
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< MatrixDouble > eshelbyStress_dx
DataForcesAndSourcesCore::EntData EntData
int nbCols
number if dof on column
boost::shared_ptr< map< int, BlockData > > blockSetsPtr
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
boost::shared_ptr< MatrixDouble > hMat
MatrixDouble invJac
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
structure grouping operators and data used for calculation of mass (convective) element \ nonlinear_e...
MatrixDouble K
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > stiffnessMat
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
bool isDiag
true if this block is on diagonal
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< MatrixDouble > invHMat
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
const double rhoN
exponent n in E(p) = E * (p / p_0)^n
std::vector< EntityHandle > & mapGaussPts
boost::shared_ptr< MatrixDouble > HMat
OpAleLhs_dX_dx(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Range tEts
constrains elements in block set
MatrixDouble transK
boost::shared_ptr< MatrixDouble > FMat
map< int, BlockData > & blockSetsPtr
ForcesAndSourcesCore::UserDataOperator UserDataOperator
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
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)
int nbIntegrationPts
number of integration points
boost::shared_ptr< map< int, BlockData > > blockSetsPtr
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
VectorDouble a0
constant acceleration
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< MatrixDouble > eshelbyStressMat
boost::shared_ptr< VectorDouble > detHVec
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
int nbRows
number of dofs on rows
boost::shared_ptr< DataAtIntegrationPts > commonData
#define MAT_TO_DDG(SM)
VectorInt colIndices
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr)
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)
const double rHo0
p_0 reference density in E(p) = E * (p / p_0)^n
#define CHKERR
Inline error check.
Definition: definitions.h:602
FTensor::Index< 'k', 2 > k
Definition: ContactOps.hpp:28
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:72
VectorDouble nF
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)
boost::shared_ptr< VectorDouble > energyVec
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1788
struct ConvectiveMassElement BlockData
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
VectorInt rowIndices
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
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.
moab::Interface & postProcMesh
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....
data for calculation heat conductivity and heat capacity elements
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
data for calculation inertia forces
MassBlockData & massData
FTensor::Index< 'l', 2 > l
Definition: ContactOps.hpp:29
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< MatrixDouble > cauchyStressMat
Assemble mass matrix for elastic element TODO: CHANGE FORMULA * .
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
double w(const double g, const double t)
Definition: ContactOps.hpp:160
structure grouping operators and data used for calculation of nonlinear elastic elementIn order to as...
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26