v0.14.0
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 
18 #ifndef __HOOKE_ELEMENT_HPP
19 #define __HOOKE_ELEMENT_HPP
20 
21 #ifndef __BASICFINITEELEMENTS_HPP__
22 #include <BasicFiniteElements.hpp>
23 #endif // __BASICFINITEELEMENTS_HPP__
24 
25 #ifndef __NONLINEAR_ELASTIC_HPP
26 
28 
29  /** \brief data for calculation heat conductivity and heat capacity elements
30  * \ingroup nonlinear_elastic_elem
31  */
32  struct BlockData {
33  int iD;
34  double E;
35  double PoissonRatio;
36  Range tEts; ///< constrains elements in block set
39  };
40 };
41 
42 #endif // __NONLINEAR_ELASTIC_HPP
43 
44 /** \brief structure grouping operators and data used for calculation of
45  * nonlinear elastic element \ingroup nonlinear_elastic_elem
46  *
47  * In order to assemble matrices and right hand vectors, the loops over
48  * elements, entities over that elements and finally loop over integration
49  * points are executed.
50  *
51  * Following implementation separate those three categories of loops and to each
52  * loop attach operator.
53  *
54  */
55 
56 #ifndef __CONVECTIVE_MASS_ELEMENT_HPP
57 struct ConvectiveMassElement {
58  /** \brief data for calculation inertia forces
59  * \ingroup user_modules
60  */
61  struct BlockData {
62  double rho0; ///< reference density
63  VectorDouble a0; ///< constant acceleration
64  Range tEts; ///< elements in block set
65  };
66 }
67 
68 #endif //__CONVECTIVE_MASS_ELEMENT_HPP
69 struct HookeElement {
70 
73 
76  using VolUserDataOperator =
78 
80 
81  boost::shared_ptr<MatrixDouble> smallStrainMat;
82  boost::shared_ptr<MatrixDouble> hMat;
83  boost::shared_ptr<MatrixDouble> FMat;
84 
85  boost::shared_ptr<MatrixDouble> HMat;
86  boost::shared_ptr<VectorDouble> detHVec;
87  boost::shared_ptr<MatrixDouble> invHMat;
88 
89  boost::shared_ptr<MatrixDouble> cauchyStressMat;
90  boost::shared_ptr<MatrixDouble> stiffnessMat;
91  boost::shared_ptr<VectorDouble> energyVec;
92  boost::shared_ptr<MatrixDouble> eshelbyStressMat;
93 
94  boost::shared_ptr<MatrixDouble> eshelbyStress_dx;
95 
97 
98  smallStrainMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
99  hMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
100  FMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
101 
102  HMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
103  detHVec = boost::shared_ptr<VectorDouble>(new VectorDouble());
104  invHMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
105 
106  cauchyStressMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
107  stiffnessMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
108  energyVec = boost::shared_ptr<VectorDouble>(new VectorDouble());
109  eshelbyStressMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
110 
111  eshelbyStress_dx = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
112  }
113 
116  };
117 
118  template <bool D = false>
120 
121  OpCalculateStrain(const std::string row_field, const std::string col_field,
122  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
123 
124  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
125 
126  private:
127  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
128  };
129 
131 
132  OpCalculateStrainAle(const std::string row_field,
133  const std::string col_field,
134  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
135 
136  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
137 
138  private:
139  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
140  };
141 
142 #define MAT_TO_DDG(SM) \
143  &(*SM)(0, 0), &(*SM)(1, 0), &(*SM)(2, 0), &(*SM)(3, 0), &(*SM)(4, 0), \
144  &(*SM)(5, 0), &(*SM)(6, 0), &(*SM)(7, 0), &(*SM)(8, 0), &(*SM)(9, 0), \
145  &(*SM)(10, 0), &(*SM)(11, 0), &(*SM)(12, 0), &(*SM)(13, 0), \
146  &(*SM)(14, 0), &(*SM)(15, 0), &(*SM)(16, 0), &(*SM)(17, 0), \
147  &(*SM)(18, 0), &(*SM)(19, 0), &(*SM)(20, 0), &(*SM)(21, 0), \
148  &(*SM)(22, 0), &(*SM)(23, 0), &(*SM)(24, 0), &(*SM)(25, 0), \
149  &(*SM)(26, 0), &(*SM)(27, 0), &(*SM)(28, 0), &(*SM)(29, 0), \
150  &(*SM)(30, 0), &(*SM)(31, 0), &(*SM)(32, 0), &(*SM)(33, 0), \
151  &(*SM)(34, 0), &(*SM)(35, 0)
152 
153  template <int S = 0> struct OpCalculateStress : public VolUserDataOperator {
154 
155  OpCalculateStress(const std::string row_field, const std::string col_field,
156  boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
157 
158  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
159 
160  protected:
161  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
162  };
163 
165 
166  OpCalculateEnergy(const std::string row_field, const std::string col_field,
167  boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
168  SmartPetscObj<Vec> ghost_vec = SmartPetscObj<Vec>());
169 
170  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
171 
172  protected:
173  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
175  };
176 
178 
180  const std::string row_field, const std::string col_field,
181  boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
182 
183  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
184 
185  protected:
186  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
187  };
188 
189  template <int S = 0>
191 
193  const std::string row_field, const std::string col_field,
194  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
195  boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
196 
197  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
198 
199  protected:
200  boost::shared_ptr<map<int, BlockData>>
201  blockSetsPtr; ///< Structure keeping data about problem, like
202  ///< material parameters
203  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
204  };
205 
206  /** * @brief Assemble mass matrix for elastic element TODO: CHANGE FORMULA *
207  * \f[
208  * {\bf{M}} = \int\limits_\Omega
209  * \f]
210  *
211  */
214 
219 
220  boost::shared_ptr<DataAtIntegrationPts> commonData;
221 
222  OpCalculateMassMatrix(const std::string row_field,
223  const std::string col_field, BlockData &data,
224  MassBlockData &mass_data,
225  boost::shared_ptr<DataAtIntegrationPts> &common_data,
226  bool symm = true)
228  row_field, col_field, OPROWCOL, symm),
229  commonData(common_data), dAta(data), massData(mass_data) {}
230 
231  PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
232  EntityType col_type,
233  EntitiesFieldData::EntData &row_data,
234  EntitiesFieldData::EntData &col_data) {
236 
237  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
239  &m(3 * r + 0, 3 * c + 0), &m(3 * r + 0, 3 * c + 1),
240  &m(3 * r + 0, 3 * c + 2), &m(3 * r + 1, 3 * c + 0),
241  &m(3 * r + 1, 3 * c + 1), &m(3 * r + 1, 3 * c + 2),
242  &m(3 * r + 2, 3 * c + 0), &m(3 * r + 2, 3 * c + 1),
243  &m(3 * r + 2, 3 * c + 2));
244  };
245 
246  const int row_nb_dofs = row_data.getIndices().size();
247  if (!row_nb_dofs)
249  const int col_nb_dofs = col_data.getIndices().size();
250  if (!col_nb_dofs)
252  if (dAta.tEts.find(getFEEntityHandle()) == dAta.tEts.end()) {
254  }
255  if (massData.tEts.find(getFEEntityHandle()) == massData.tEts.end()) {
257  }
258 
259  const bool diagonal_block =
260  (row_type == col_type) && (row_side == col_side);
261  // get number of integration points
262  // Set size can clear local tangent matrix
263  locK.resize(row_nb_dofs, col_nb_dofs, false);
264  locK.clear();
265 
266  const int row_nb_gauss_pts = row_data.getN().size1();
267  const int row_nb_base_functions = row_data.getN().size2();
268 
273 
274  double density = massData.rho0;
275 
276  // get integration weights
277  auto t_w = getFTensor0IntegrationWeight();
278 
279  // integrate local matrix for entity block
280  for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
281 
282  auto t_row_base_func = row_data.getFTensor0N(gg, 0);
283 
284  // Get volume and integration weight
285  double w = getVolume() * t_w;
286 
287  for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
288  auto t_col_base_func = col_data.getFTensor0N(gg, 0);
289  for (int col_bb = 0; col_bb != col_nb_dofs / 3; col_bb++) {
290  auto t_assemble = get_tensor2(locK, row_bb, col_bb);
291  t_assemble(i, j) += density * t_row_base_func * t_col_base_func * w;
292  // Next base function for column
293  ++t_col_base_func;
294  }
295  // Next base function for row
296  ++t_row_base_func;
297  }
298  // Next integration point for getting weight
299  ++t_w;
300  }
301 
302  CHKERR MatSetValues(getKSPB(), row_data, col_data, &locK(0, 0),
303  ADD_VALUES);
304 
305  // is symmetric
306  if (row_type != col_type || row_side != col_side) {
307  translocK.resize(col_nb_dofs, row_nb_dofs, false);
308  noalias(translocK) = trans(locK);
309 
310  CHKERR MatSetValues(getKSPB(), col_data, row_data, &translocK(0, 0),
311  ADD_VALUES);
312  }
313 
315  }
316  };
317 
319  protected:
320  boost::shared_ptr<map<int, BlockData>>
321  blockSetsPtr; ///< Structure keeping data about problem, like
322  ///< material parameters
323  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
324 
325  boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
326  const double rhoN; ///< exponent n in E(p) = E * (p / p_0)^n
327  const double rHo0; ///< p_0 reference density in E(p) = E * (p / p_0)^n
328  // // where p is density, E - youngs modulus
329  public:
331  const std::string row_field, const std::string col_field,
332  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
333  boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
334  boost::shared_ptr<VectorDouble> rho_at_gauss_pts, const double rho_n,
335  const double rho_0);
336 
337  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
338  };
339 
340  struct OpAssemble : public VolUserDataOperator {
341 
342  OpAssemble(const std::string row_field, const std::string col_field,
343  boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
344  const char type, bool symm = false);
345 
346  /**
347  * \brief Do calculations for give operator
348  * @param row_side row side number (local number) of entity on element
349  * @param col_side column side number (local number) of entity on element
350  * @param row_type type of row entity MBVERTEX, MBEDGE, MBTRI or MBTET
351  * @param col_type type of column entity MBVERTEX, MBEDGE, MBTRI or MBTET
352  * @param row_data data for row
353  * @param col_data data for column
354  * @return error code
355  */
356  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
357  EntityType col_type, EntData &row_data,
358  EntData &col_data);
359 
360  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
361 
362  protected:
363  // Finite element stiffness sub-matrix K_ij
367 
368  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
369 
372 
373  int nbRows; ///< number of dofs on rows
374  int nbCols; ///< number if dof on column
375  int nbIntegrationPts; ///< number of integration points
376  bool isDiag; ///< true if this block is on diagonal
377 
378  virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
379 
380  virtual MoFEMErrorCode iNtegrate(EntData &row_data);
381 
382  /**
383  * \brief Assemble local entity block matrix
384  * @param row_data row data (consist base functions on row entity)
385  * @param col_data column data (consist base functions on column
386  * entity)
387  * @return error code
388  */
389  MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data);
390 
391  /**
392  * \brief Assemble local entity right-hand vector
393  * @param row_data row data (consist base functions on row entity)
394  * @param col_data column data (consist base functions on column
395  * entity)
396  * @return error code
397  */
398  MoFEMErrorCode aSsemble(EntData &row_data);
399  };
400 
401  struct OpRhs_dx : public OpAssemble {
402 
403  OpRhs_dx(const std::string row_field, const std::string col_field,
404  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
405 
406  protected:
407  MoFEMErrorCode iNtegrate(EntData &row_data);
408  };
409 
410  template <int S = 0> struct OpLhs_dx_dx : public OpAssemble {
411 
412  OpLhs_dx_dx(const std::string row_field, const std::string col_field,
413  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
414 
415  protected:
416  /**
417  * \brief Integrate B^T D B operator
418  * @param row_data row data (consist base functions on row entity)
419  * @param col_data column data (consist base functions on column entity)
420  * @return error code
421  */
422  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
423  };
424 
425  struct OpAleRhs_dx : public OpAssemble {
426 
427  OpAleRhs_dx(const std::string row_field, const std::string col_field,
428  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
429 
430  protected:
431  MoFEMErrorCode iNtegrate(EntData &row_data);
432  };
433 
434  template <int S = 0> struct OpAleLhs_dx_dx : public OpAssemble {
435 
436  OpAleLhs_dx_dx(const std::string row_field, const std::string col_field,
437  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
438 
439  protected:
440  /**
441  * \brief Integrate B^T D B operator
442  * @param row_data row data (consist base functions on row entity)
443  * @param col_data column data (consist base functions on column entity)
444  * @return error code
445  */
446  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
447  };
448 
449  template <int S = 0> struct OpAleLhs_dx_dX : public OpAssemble {
450 
451  OpAleLhs_dx_dX(const std::string row_field, const std::string col_field,
452  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
453 
454  protected:
455  /**
456  * \brief Integrate tangent stiffness for spatial momentum
457  * @param row_data row data (consist base functions on row entity)
458  * @param col_data column data (consist base functions on column entity)
459  * @return error code
460  */
461  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
462  };
463 
465 
466  boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
467  boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
468  const double rhoN;
469  const double rHo0;
470 
472  const std::string row_field, const std::string col_field,
473  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
474  boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
475  boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
476  const double rho_n, const double rho_0);
477 
478  protected:
479  /**
480  * \brief Integrate tangent stiffness for spatial momentum
481  * @param row_data row data (consist base functions on row entity)
482  * @param col_data column data (consist base functions on column entity)
483  * @return error code
484  */
485  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
486  };
487 
489 
490  boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
491  boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
492  const double rhoN;
493  const double rHo0;
494 
496  const std::string row_field, const std::string col_field,
497  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
498  boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
499  boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
500  const double rho_n, const double rho_0);
501 
502  protected:
503  /**
504  * \brief Integrate tangent stiffness for material momentum
505  * @param row_data row data (consist base functions on row entity)
506  * @param col_data column data (consist base functions on column entity)
507  * @return error code
508  */
509  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
510  };
511 
512  struct OpAleRhs_dX : public OpAssemble {
513 
514  OpAleRhs_dX(const std::string row_field, const std::string col_field,
515  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
516 
517  protected:
518  MoFEMErrorCode iNtegrate(EntData &row_data);
519  };
520 
521  template <int S = 0> struct OpAleLhs_dX_dX : public OpAssemble {
522 
523  OpAleLhs_dX_dX(const std::string row_field, const std::string col_field,
524  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
525 
526  protected:
527  /**
528  * \brief Integrate tangent stiffness for material momentum
529  * @param row_data row data (consist base functions on row entity)
530  * @param col_data column data (consist base functions on column entity)
531  * @return error code
532  */
533  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
534  };
535 
536  template <int S = 0> struct OpAleLhsPre_dX_dx : public VolUserDataOperator {
537 
538  OpAleLhsPre_dX_dx(const std::string row_field, const std::string col_field,
539  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
540 
541  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
542 
543  private:
544  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
545  };
546 
547  struct OpAleLhs_dX_dx : public OpAssemble {
548 
549  OpAleLhs_dX_dx(const std::string row_field, const std::string col_field,
550  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
551  : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
552 
553  protected:
554  /**
555  * \brief Integrate tangent stiffness for material momentum
556  * @param row_data row data (consist base functions on row entity)
557  * @param col_data column data (consist base functions on column entity)
558  * @return error code
559  */
560  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
561  };
562 
563  template <int S> struct OpAnalyticalInternalStrain_dx : public OpAssemble {
564 
565  typedef boost::function<
566 
568 
570 
571  )
572 
573  >
575 
577  const std::string row_field,
578  boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
579  StrainFunction strain_fun);
580 
581  protected:
582  MoFEMErrorCode iNtegrate(EntData &row_data);
584  };
585 
586  template <int S> struct OpAnalyticalInternalAleStrain_dX : public OpAssemble {
587 
588  typedef boost::function<
589 
591 
593 
594  )
595 
596  >
598 
600  const std::string row_field,
601  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
602  StrainFunction strain_fun,
603  boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr);
604 
605  protected:
606  MoFEMErrorCode iNtegrate(EntData &row_data);
608  boost::shared_ptr<MatrixDouble> matPosAtPtsPtr;
609  };
610 
611  template <int S> struct OpAnalyticalInternalAleStrain_dx : public OpAssemble {
612 
613  typedef boost::function<
614 
616 
618 
619  )
620 
621  >
623 
625  const std::string row_field,
626  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
627  StrainFunction strain_fun,
628  boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr);
629 
630  protected:
631  MoFEMErrorCode iNtegrate(EntData &row_data);
633  boost::shared_ptr<MatrixDouble> matPosAtPtsPtr;
634  };
635 
636  template <class ELEMENT>
638  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
639  map<int, BlockData>
640  &blockSetsPtr; // FIXME: (works only with the first block)
642  std::vector<EntityHandle> &mapGaussPts;
643  bool isALE;
645 
646  OpPostProcHookeElement(const string row_field,
647  boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
648  map<int, BlockData> &block_sets_ptr,
649  moab::Interface &post_proc_mesh,
650  std::vector<EntityHandle> &map_gauss_pts,
651  bool is_ale = false, bool is_field_disp = true);
652 
653  MoFEMErrorCode doWork(int side, EntityType type,
655  };
656 
657  static MoFEMErrorCode
658  setBlocks(MoFEM::Interface &m_field,
659  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr);
660 
661  static MoFEMErrorCode
663  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
664  const std::string element_name, const std::string x_field,
665  const std::string X_field, const bool ale);
666 
667  static MoFEMErrorCode
668  setOperators(boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr,
669  boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr,
670  boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
671  const std::string x_field, const std::string X_field,
672  const bool ale, const bool field_disp,
673  const EntityType type = MBTET,
674  boost::shared_ptr<DataAtIntegrationPts> data_at_pts = nullptr);
675 
676  static MoFEMErrorCode
677  calculateEnergy(DM dm, boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
678  const std::string x_field, const std::string X_field,
679  const bool ale, const bool field_disp,
680  SmartPetscObj<Vec> &v_energy_ptr);
681 
682 private:
684 };
685 
686 template <bool D>
687 HookeElement::OpCalculateStrain<D>::OpCalculateStrain(
688  const std::string row_field, const std::string col_field,
689  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
690  : VolUserDataOperator(row_field, col_field, OPROW, true),
691  dataAtPts(data_at_pts) {
692  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
693 }
694 
695 template <bool D>
696 MoFEMErrorCode HookeElement::OpCalculateStrain<D>::doWork(int row_side,
697  EntityType row_type,
698  EntData &row_data) {
702  // get number of integration points
703  const int nb_integration_pts = getGaussPts().size2();
704  dataAtPts->smallStrainMat->resize(6, nb_integration_pts, false);
705  auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
706  auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hMat));
707 
708  for (int gg = 0; gg != nb_integration_pts; ++gg) {
709  t_strain(i, j) = (t_h(i, j) || t_h(j, i)) / 2.;
710 
711  // If displacement field, not field o spatial positons is given
712  if (!D) {
713  t_strain(0, 0) -= 1;
714  t_strain(1, 1) -= 1;
715  t_strain(2, 2) -= 1;
716  }
717 
718  ++t_strain;
719  ++t_h;
720  }
722 }
723 
724 template <int S>
725 HookeElement::OpAleLhs_dx_dx<S>::OpAleLhs_dx_dx(
726  const std::string row_field, const std::string col_field,
727  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
728  : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
729 
730 template <int S>
731 MoFEMErrorCode HookeElement::OpAleLhs_dx_dx<S>::iNtegrate(EntData &row_data,
732  EntData &col_data) {
734 
735  // get sub-block (3x3) of local stiffens matrix, here represented by
736  // second order tensor
737  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
739  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
740  &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
741  &m(r + 2, c + 2));
742  };
743 
748 
749  // get element volume
750  double vol = getVolume();
751 
752  // get intergrayion weights
753  auto t_w = getFTensor0IntegrationWeight();
754 
755  // get derivatives of base functions on rows
756  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
757  const int row_nb_base_fun = row_data.getN().size2();
758 
759  // Elastic stiffness tensor (4th rank tensor with minor and major
760  // symmetry)
762  MAT_TO_DDG(dataAtPts->stiffnessMat));
763 
764  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
765  auto &det_H = *dataAtPts->detHVec;
766 
767  // iterate over integration points
768  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
769 
770  // calculate scalar weight times element volume
771  double a = t_w * vol * det_H[gg];
772 
773  // iterate over row base functions
774  int rr = 0;
775  for (; rr != nbRows / 3; ++rr) {
776 
777  // get sub matrix for the row
778  auto t_m = get_tensor2(K, 3 * rr, 0);
779 
780  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
781  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
782 
784  // I mix up the indices here so that it behaves like a
785  // Dg. That way I don't have to have a separate wrapper
786  // class Christof_Expr, which simplifies things.
787  t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base_pulled(i));
788 
789  // get derivatives of base functions for columns
790  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
791 
792  // iterate column base functions
793  for (int cc = 0; cc != nbCols / 3; ++cc) {
794 
795  FTensor::Tensor1<double, 3> t_col_diff_base_pulled;
796  t_col_diff_base_pulled(j) = t_col_diff_base(i) * t_invH(i, j);
797 
798  // integrate block local stiffens matrix
799  t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base_pulled(k);
800 
801  // move to next column base function
802  ++t_col_diff_base;
803 
804  // move to next block of local stiffens matrix
805  ++t_m;
806  }
807 
808  // move to next row base function
809  ++t_row_diff_base;
810  }
811 
812  for (; rr != row_nb_base_fun; ++rr)
813  ++t_row_diff_base;
814 
815  // move to next integration weight
816  ++t_w;
817  ++t_D;
818  ++t_invH;
819  }
820 
822 }
823 
824 template <int S>
825 HookeElement::OpCalculateHomogeneousStiffness<S>::
826  OpCalculateHomogeneousStiffness(
827  const std::string row_field, const std::string col_field,
828  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
829  boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
830  : VolUserDataOperator(row_field, col_field, OPROW, true),
831  blockSetsPtr(block_sets_ptr), dataAtPts(data_at_pts) {
832  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
833 }
834 
835 template <int S>
836 MoFEMErrorCode HookeElement::OpCalculateHomogeneousStiffness<S>::doWork(
837  int row_side, EntityType row_type, EntData &row_data) {
839 
840  for (auto &m : (*blockSetsPtr)) {
841  if (m.second.tEts.find(getFEEntityHandle()) != m.second.tEts.end()) {
842 
843  dataAtPts->stiffnessMat->resize(36, 1, false);
845  MAT_TO_DDG(dataAtPts->stiffnessMat));
846  const double young = m.second.E;
847  const double poisson = m.second.PoissonRatio;
848 
849  // coefficient used in intermediate calculation
850  const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
851 
856 
857  t_D(i, j, k, l) = 0.;
858 
859  t_D(0, 0, 0, 0) = 1 - poisson;
860  t_D(1, 1, 1, 1) = 1 - poisson;
861  t_D(2, 2, 2, 2) = 1 - poisson;
862 
863  t_D(0, 1, 0, 1) = 0.5 * (1 - 2 * poisson);
864  t_D(0, 2, 0, 2) = 0.5 * (1 - 2 * poisson);
865  t_D(1, 2, 1, 2) = 0.5 * (1 - 2 * poisson);
866 
867  t_D(0, 0, 1, 1) = poisson;
868  t_D(1, 1, 0, 0) = poisson;
869  t_D(0, 0, 2, 2) = poisson;
870  t_D(2, 2, 0, 0) = poisson;
871  t_D(1, 1, 2, 2) = poisson;
872  t_D(2, 2, 1, 1) = poisson;
873 
874  t_D(i, j, k, l) *= coefficient;
875 
876  break;
877  }
878  }
879 
881 }
882 
883 template <int S>
884 HookeElement::OpCalculateStress<S>::OpCalculateStress(
885  const std::string row_field, const std::string col_field,
886  boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
887  : VolUserDataOperator(row_field, col_field, OPROW, true),
888  dataAtPts(data_at_pts) {
889  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
890 }
891 
892 template <int S>
893 MoFEMErrorCode HookeElement::OpCalculateStress<S>::doWork(int row_side,
894  EntityType row_type,
895  EntData &row_data) {
897  // get number of integration points
898  const int nb_integration_pts = getGaussPts().size2();
899  auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
900  dataAtPts->cauchyStressMat->resize(6, nb_integration_pts, false);
901  auto t_cauchy_stress =
902  getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
903 
908 
909  // elastic stiffness tensor (4th rank tensor with minor and major
910  // symmetry)
912  MAT_TO_DDG(dataAtPts->stiffnessMat));
913  for (int gg = 0; gg != nb_integration_pts; ++gg) {
914  t_cauchy_stress(i, j) = t_D(i, j, k, l) * t_strain(k, l);
915  ++t_strain;
916  ++t_cauchy_stress;
917  ++t_D;
918  }
920 }
921 
922 template <int S>
923 HookeElement::OpLhs_dx_dx<S>::OpLhs_dx_dx(
924  const std::string row_field, const std::string col_field,
925  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
926  : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
927 
928 template <int S>
929 MoFEMErrorCode HookeElement::OpLhs_dx_dx<S>::iNtegrate(EntData &row_data,
930  EntData &col_data) {
932 
933  // get sub-block (3x3) of local stiffens matrix, here represented by
934  // second order tensor
935  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
937  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
938  &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
939  &m(r + 2, c + 2));
940  };
941 
946 
947  // get element volume
948  double vol = getVolume();
949 
950  // get intergrayion weights
951  auto t_w = getFTensor0IntegrationWeight();
952 
953  // get derivatives of base functions on rows
954  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
955  const int row_nb_base_fun = row_data.getN().size2();
956 
957  // Elastic stiffness tensor (4th rank tensor with minor and major
958  // symmetry)
960  MAT_TO_DDG(dataAtPts->stiffnessMat));
961 
962  // iterate over integration points
963  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
964 
965  // calculate scalar weight times element volume
966  double a = t_w * vol;
967 
968  // iterate over row base functions
969  int rr = 0;
970  for (; rr != nbRows / 3; ++rr) {
971 
972  // get sub matrix for the row
973  auto t_m = get_tensor2(K, 3 * rr, 0);
974 
975  // get derivatives of base functions for columns
976  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
977 
979  // I mix up the indices here so that it behaves like a
980  // Dg. That way I don't have to have a separate wrapper
981  // class Christof_Expr, which simplifies things.
982  t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base(i));
983 
984  // iterate column base functions
985  for (int cc = 0; cc != nbCols / 3; ++cc) {
986 
987  // integrate block local stiffens matrix
988  t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base(k);
989 
990  // move to next column base function
991  ++t_col_diff_base;
992 
993  // move to next block of local stiffens matrix
994  ++t_m;
995  }
996 
997  // move to next row base function
998  ++t_row_diff_base;
999  }
1000 
1001  for (; rr != row_nb_base_fun; ++rr)
1002  ++t_row_diff_base;
1003 
1004  // move to next integration weight
1005  ++t_w;
1006  ++t_D;
1007  }
1008 
1010 }
1011 
1012 template <int S>
1013 HookeElement::OpAleLhs_dx_dX<S>::OpAleLhs_dx_dX(
1014  const std::string row_field, const std::string col_field,
1015  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1016  : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
1017 
1018 template <int S>
1019 MoFEMErrorCode HookeElement::OpAleLhs_dx_dX<S>::iNtegrate(EntData &row_data,
1020  EntData &col_data) {
1022 
1023  // get sub-block (3x3) of local stiffens matrix, here represented by
1024  // second order tensor
1025  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1027  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1028  &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1029  &m(r + 2, c + 2));
1030  };
1031 
1038 
1039  // get element volume
1040  double vol = getVolume();
1041 
1042  // get intergrayion weights
1043  auto t_w = getFTensor0IntegrationWeight();
1044 
1045  // get derivatives of base functions on rows
1046  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1047  const int row_nb_base_fun = row_data.getN().size2();
1048 
1049  // Elastic stiffness tensor (4th rank tensor with minor and major
1050  // symmetry)
1052  MAT_TO_DDG(dataAtPts->stiffnessMat));
1053 
1054  auto t_cauchy_stress =
1055  getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1056  auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1057  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1058  auto &det_H = *dataAtPts->detHVec;
1059 
1060  // iterate over integration points
1061  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1062 
1063  // calculate scalar weight times element volume
1064  double a = t_w * vol * det_H[gg];
1065 
1067  t_F_dX(i, j, k, l) = -(t_h(i, m) * t_invH(m, k)) * t_invH(l, j);
1068 
1069  // iterate over row base functions
1070  int rr = 0;
1071  for (; rr != nbRows / 3; ++rr) {
1072 
1073  // get sub matrix for the row
1074  auto t_m = get_tensor2(K, 3 * rr, 0);
1075 
1076  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1077  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1078 
1079  FTensor::Tensor1<double, 3> t_row_stress;
1080  t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_cauchy_stress(i, j);
1081 
1082  FTensor::Tensor3<double, 3, 3, 3> t_row_diff_base_pulled_dX;
1083  t_row_diff_base_pulled_dX(j, k, l) =
1084  -(t_invH(i, k) * t_row_diff_base(i)) * t_invH(l, j);
1085 
1086  FTensor::Tensor3<double, 3, 3, 3> t_row_dX_stress;
1087  t_row_dX_stress(i, k, l) =
1088  a * (t_row_diff_base_pulled_dX(j, k, l) * t_cauchy_stress(j, i));
1089 
1091  t_row_D(l, j, k) = (a * t_row_diff_base_pulled(i)) * t_D(i, j, k, l);
1092 
1093  FTensor::Tensor3<double, 3, 3, 3> t_row_stress_dX;
1094  // FIXME: This operator is not implemented, doing operation by hand
1095  // t_row_stress_dX(i, m, n) = t_row_D(i, k, l) * t_F_dX(k, l, m, n);
1096  t_row_stress_dX(i, j, k) = 0;
1097  for (int ii = 0; ii != 3; ++ii)
1098  for (int mm = 0; mm != 3; ++mm)
1099  for (int nn = 0; nn != 3; ++nn) {
1100  auto &v = t_row_stress_dX(ii, mm, nn);
1101  for (int kk = 0; kk != 3; ++kk)
1102  for (int ll = 0; ll != 3; ++ll)
1103  v += t_row_D(ii, kk, ll) * t_F_dX(kk, ll, mm, nn);
1104  }
1105 
1106  // get derivatives of base functions for columns
1107  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1108 
1109  // iterate column base functions
1110  for (int cc = 0; cc != nbCols / 3; ++cc) {
1111 
1112  t_m(i, k) += t_row_stress(i) * (t_invH(j, k) * t_col_diff_base(j));
1113  t_m(i, k) += t_row_dX_stress(i, k, l) * t_col_diff_base(l);
1114  t_m(i, k) += t_row_stress_dX(i, k, l) * t_col_diff_base(l);
1115 
1116  // move to next column base function
1117  ++t_col_diff_base;
1118 
1119  // move to next block of local stiffens matrix
1120  ++t_m;
1121  }
1122 
1123  // move to next row base function
1124  ++t_row_diff_base;
1125  }
1126 
1127  for (; rr != row_nb_base_fun; ++rr)
1128  ++t_row_diff_base;
1129 
1130  // move to next integration weight
1131  ++t_w;
1132  ++t_D;
1133  ++t_cauchy_stress;
1134  ++t_invH;
1135  ++t_h;
1136  }
1137 
1139 }
1140 
1141 template <int S>
1142 HookeElement::OpAleLhs_dX_dX<S>::OpAleLhs_dX_dX(
1143  const std::string row_field, const std::string col_field,
1144  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1145  : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
1146 
1147 template <int S>
1148 MoFEMErrorCode HookeElement::OpAleLhs_dX_dX<S>::iNtegrate(EntData &row_data,
1149  EntData &col_data) {
1151 
1152  // get sub-block (3x3) of local stiffens matrix, here represented by
1153  // second order tensor
1154  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1156  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1157  &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1158  &m(r + 2, c + 2));
1159  };
1160 
1167 
1168  // get element volume
1169  double vol = getVolume();
1170 
1171  // get intergrayion weights
1172  auto t_w = getFTensor0IntegrationWeight();
1173 
1174  // get derivatives of base functions on rows
1175  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1176  const int row_nb_base_fun = row_data.getN().size2();
1177 
1178  // Elastic stiffness tensor (4th rank tensor with minor and major
1179  // symmetry)
1181  MAT_TO_DDG(dataAtPts->stiffnessMat));
1182  auto t_cauchy_stress =
1183  getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1184  auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
1185  auto t_eshelby_stress =
1186  getFTensor2FromMat<3, 3>(*dataAtPts->eshelbyStressMat);
1187  auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1188  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1189  auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
1190  auto &det_H = *dataAtPts->detHVec;
1191 
1192  // iterate over integration points
1193  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1194 
1195  // calculate scalar weight times element volume
1196  double a = t_w * vol * det_H[gg];
1197 
1199  t_F_dX(i, j, k, l) = -(t_h(i, m) * t_invH(m, k)) * t_invH(l, j);
1200 
1202  t_D_strain_dX(i, j, m, n) = 0.;
1203  for (int ii = 0; ii != 3; ++ii)
1204  for (int jj = 0; jj != 3; ++jj)
1205  for (int ll = 0; ll != 3; ++ll)
1206  for (int kk = 0; kk != 3; ++kk) {
1207  auto &v = t_D_strain_dX(ii, jj, kk, ll);
1208  for (int mm = 0; mm != 3; ++mm)
1209  for (int nn = 0; nn != 3; ++nn)
1210  v += t_D(ii, jj, mm, nn) * t_F_dX(mm, nn, kk, ll);
1211  }
1212 
1213  FTensor::Tensor4<double, 3, 3, 3, 3> t_eshelby_stress_dX;
1214  t_eshelby_stress_dX(i, j, m, n) = t_F(k, i) * t_D_strain_dX(k, j, m, n);
1215 
1216  for (int ii = 0; ii != 3; ++ii)
1217  for (int jj = 0; jj != 3; ++jj)
1218  for (int mm = 0; mm != 3; ++mm)
1219  for (int nn = 0; nn != 3; ++nn) {
1220  auto &v = t_eshelby_stress_dX(ii, jj, mm, nn);
1221  for (int kk = 0; kk != 3; ++kk)
1222  v += t_F_dX(kk, ii, mm, nn) * t_cauchy_stress(kk, jj);
1223  }
1224 
1225  t_eshelby_stress_dX(i, j, k, l) *= -1;
1226 
1227  FTensor::Tensor2<double, 3, 3> t_energy_dX;
1228  t_energy_dX(k, l) = t_F_dX(i, j, k, l) * t_cauchy_stress(i, j);
1229  t_energy_dX(k, l) +=
1230  (t_strain(m, n) * t_D(m, n, i, j)) * t_F_dX(i, j, k, l);
1231  t_energy_dX(k, l) /= 2.;
1232 
1233  for (int kk = 0; kk != 3; ++kk)
1234  for (int ll = 0; ll != 3; ++ll) {
1235  auto v = t_energy_dX(kk, ll);
1236  for (int ii = 0; ii != 3; ++ii)
1237  t_eshelby_stress_dX(ii, ii, kk, ll) += v;
1238  }
1239 
1240  // iterate over row base functions
1241  int rr = 0;
1242  for (; rr != nbRows / 3; ++rr) {
1243 
1244  // get sub matrix for the row
1245  auto t_m = get_tensor2(K, 3 * rr, 0);
1246 
1247  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1248  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1249 
1250  FTensor::Tensor1<double, 3> t_row_stress;
1251  t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
1252 
1253  FTensor::Tensor3<double, 3, 3, 3> t_row_diff_base_pulled_dX;
1254  t_row_diff_base_pulled_dX(j, k, l) =
1255  -(t_row_diff_base(i) * t_invH(i, k)) * t_invH(l, j);
1256 
1257  FTensor::Tensor3<double, 3, 3, 3> t_row_dX_stress;
1258  t_row_dX_stress(i, k, l) =
1259  a * (t_row_diff_base_pulled_dX(j, k, l) * t_eshelby_stress(i, j));
1260 
1261  FTensor::Tensor3<double, 3, 3, 3> t_row_stress_dX;
1262  t_row_stress_dX(i, m, n) =
1263  a * t_row_diff_base_pulled(j) * t_eshelby_stress_dX(i, j, m, n);
1264 
1265  // get derivatives of base functions for columns
1266  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1267 
1268  // iterate column base functions
1269  for (int cc = 0; cc != nbCols / 3; ++cc) {
1270 
1271  t_m(i, k) += t_row_stress(i) * (t_invH(j, k) * t_col_diff_base(j));
1272  t_m(i, k) += t_row_dX_stress(i, k, l) * t_col_diff_base(l);
1273  t_m(i, k) += t_row_stress_dX(i, k, l) * t_col_diff_base(l);
1274 
1275  // move to next column base function
1276  ++t_col_diff_base;
1277 
1278  // move to next block of local stiffens matrix
1279  ++t_m;
1280  }
1281 
1282  // move to next row base function
1283  ++t_row_diff_base;
1284  }
1285 
1286  for (; rr != row_nb_base_fun; ++rr)
1287  ++t_row_diff_base;
1288 
1289  // move to next integration weight
1290  ++t_w;
1291  ++t_D;
1292  ++t_cauchy_stress;
1293  ++t_strain;
1294  ++t_eshelby_stress;
1295  ++t_h;
1296  ++t_invH;
1297  ++t_F;
1298  }
1299 
1301 }
1302 
1303 template <int S>
1304 HookeElement::OpAleLhsPre_dX_dx<S>::OpAleLhsPre_dX_dx(
1305  const std::string row_field, const std::string col_field,
1306  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1307  : VolUserDataOperator(row_field, col_field, OPROW, true),
1308  dataAtPts(data_at_pts) {
1309  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1310 }
1311 
1312 template <int S>
1313 MoFEMErrorCode HookeElement::OpAleLhsPre_dX_dx<S>::doWork(int row_side,
1314  EntityType row_type,
1315  EntData &row_data) {
1317 
1318  const int nb_integration_pts = row_data.getN().size1();
1319 
1320  auto get_eshelby_stress_dx = [this, nb_integration_pts]() {
1322  t_eshelby_stress_dx;
1323  dataAtPts->eshelbyStress_dx->resize(81, nb_integration_pts, false);
1324  int mm = 0;
1325  for (int ii = 0; ii != 3; ++ii)
1326  for (int jj = 0; jj != 3; ++jj)
1327  for (int kk = 0; kk != 3; ++kk)
1328  for (int ll = 0; ll != 3; ++ll)
1329  t_eshelby_stress_dx.ptr(ii, jj, kk, ll) =
1330  &(*dataAtPts->eshelbyStress_dx)(mm++, 0);
1331  return t_eshelby_stress_dx;
1332  };
1333 
1334  auto t_eshelby_stress_dx = get_eshelby_stress_dx();
1335 
1342 
1343  // Elastic stiffness tensor (4th rank tensor with minor and major
1344  // symmetry)
1346  MAT_TO_DDG(dataAtPts->stiffnessMat));
1347  auto t_cauchy_stress =
1348  getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1349  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1350  auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
1351 
1352  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1353 
1354  t_eshelby_stress_dx(i, j, m, n) =
1355  (t_F(k, i) * t_D(k, j, m, l)) * t_invH(n, l);
1356  for (int ii = 0; ii != 3; ++ii)
1357  for (int jj = 0; jj != 3; ++jj)
1358  for (int mm = 0; mm != 3; ++mm)
1359  for (int nn = 0; nn != 3; ++nn) {
1360  auto &v = t_eshelby_stress_dx(ii, jj, mm, nn);
1361  v += t_invH(nn, ii) * t_cauchy_stress(mm, jj);
1362  }
1363  t_eshelby_stress_dx(i, j, k, l) *= -1;
1364 
1365  FTensor::Tensor2<double, 3, 3> t_energy_dx;
1366  t_energy_dx(m, n) = t_invH(n, j) * t_cauchy_stress(m, j);
1367 
1368  for (int mm = 0; mm != 3; ++mm)
1369  for (int nn = 0; nn != 3; ++nn) {
1370  auto v = t_energy_dx(mm, nn);
1371  for (int ii = 0; ii != 3; ++ii)
1372  t_eshelby_stress_dx(ii, ii, mm, nn) += v;
1373  }
1374 
1375  ++t_D;
1376  ++t_invH;
1377  ++t_cauchy_stress;
1378  ++t_eshelby_stress_dx;
1379  ++t_F;
1380  }
1381 
1383 }
1384 
1385 template <class ELEMENT>
1386 HookeElement::OpPostProcHookeElement<ELEMENT>::OpPostProcHookeElement(
1387  const string row_field, boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1388  map<int, BlockData> &block_sets_ptr, moab::Interface &post_proc_mesh,
1389  std::vector<EntityHandle> &map_gauss_pts, bool is_ale, bool is_field_disp)
1390  : ELEMENT::UserDataOperator(row_field, UserDataOperator::OPROW),
1391  dataAtPts(data_at_pts), blockSetsPtr(block_sets_ptr),
1392  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts), isALE(is_ale),
1393  isFieldDisp(is_field_disp) {}
1394 
1395 template <class ELEMENT>
1396 MoFEMErrorCode HookeElement::OpPostProcHookeElement<ELEMENT>::doWork(
1397  int side, EntityType type, EntitiesFieldData::EntData &data) {
1399 
1400  if (type != MBVERTEX) {
1402  }
1403 
1404  auto tensor_to_tensor = [](const auto &t1, auto &t2) {
1405  t2(0, 0) = t1(0, 0);
1406  t2(1, 1) = t1(1, 1);
1407  t2(2, 2) = t1(2, 2);
1408  t2(0, 1) = t2(1, 0) = t1(1, 0);
1409  t2(0, 2) = t2(2, 0) = t1(2, 0);
1410  t2(1, 2) = t2(2, 1) = t1(2, 1);
1411  };
1412 
1413  std::array<double, 9> def_val;
1414  def_val.fill(0);
1415 
1416  auto make_tag = [&](auto name, auto size) {
1417  Tag th;
1418  CHKERR postProcMesh.tag_get_handle(name, size, MB_TYPE_DOUBLE, th,
1419  MB_TAG_CREAT | MB_TAG_SPARSE,
1420  def_val.data());
1421  return th;
1422  };
1423 
1424  auto th_stress = make_tag("STRESS", 9);
1425  auto th_psi = make_tag("ENERGY", 1);
1426 
1427  const int nb_integration_pts = mapGaussPts.size();
1428 
1433 
1434  auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1435  auto t_H = getFTensor2FromMat<3, 3>(*dataAtPts->HMat);
1436 
1437  dataAtPts->stiffnessMat->resize(36, 1, false);
1439  MAT_TO_DDG(dataAtPts->stiffnessMat));
1440 
1441  EntityHandle ent = this->getFEEntityHandle();
1442  auto type = type_from_handle(ent);
1443  EntityHandle ent_3d = ent;
1444  if (type == MBTRI || type == MBQUAD) {
1445  Range ents;
1446  auto &m_field = this->getPtrFE()->mField;
1447  CHKERR m_field.get_moab().get_adjacencies(&ent, 1, 3, false, ents,
1448  moab::Interface::UNION);
1449 #ifndef NDEBUG
1450  if (ents.empty())
1451  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1452  "Could not find a 3D element adjacent to a given face element");
1453 #endif
1454  ent_3d = ents.front();
1455  }
1456 
1457  bool found_block = false;
1458  for (auto &m : (blockSetsPtr)) {
1459  if (m.second.tEts.find(ent_3d) != m.second.tEts.end()) {
1460  const double young = m.second.E;
1461  const double poisson = m.second.PoissonRatio;
1462  const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
1463 
1464  t_D(i, j, k, l) = 0.;
1465  t_D(0, 0, 0, 0) = t_D(1, 1, 1, 1) = t_D(2, 2, 2, 2) = 1 - poisson;
1466  t_D(0, 1, 0, 1) = t_D(0, 2, 0, 2) = t_D(1, 2, 1, 2) =
1467  0.5 * (1 - 2 * poisson);
1468  t_D(0, 0, 1, 1) = t_D(1, 1, 0, 0) = t_D(0, 0, 2, 2) = t_D(2, 2, 0, 0) =
1469  t_D(1, 1, 2, 2) = t_D(2, 2, 1, 1) = poisson;
1470  t_D(i, j, k, l) *= coefficient;
1471 
1472  found_block = true;
1473  break;
1474  }
1475  }
1476  if (!found_block)
1477  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1478  "Element not found in any of material blocksets");
1479 
1480  double detH = 0.;
1484  FTensor::Tensor2<double, 3, 3> t_small_strain;
1486  FTensor::Tensor2_symmetric<double, 3> t_small_strain_symm;
1487 
1488  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1489 
1490  if (isFieldDisp) {
1491  t_h(0, 0) += 1;
1492  t_h(1, 1) += 1;
1493  t_h(2, 2) += 1;
1494  }
1495 
1496  if (!isALE) {
1497  t_small_strain_symm(i, j) = (t_h(i, j) || t_h(j, i)) / 2.;
1498  } else {
1499  CHKERR determinantTensor3by3(t_H, detH);
1500  CHKERR invertTensor3by3(t_H, detH, t_invH);
1501  t_F(i, j) = t_h(i, k) * t_invH(k, j);
1502  t_small_strain_symm(i, j) = (t_F(i, j) || t_F(j, i)) / 2.;
1503  ++t_H;
1504  }
1505 
1506  t_small_strain_symm(0, 0) -= 1;
1507  t_small_strain_symm(1, 1) -= 1;
1508  t_small_strain_symm(2, 2) -= 1;
1509 
1510  // symmetric tensors need improvement
1511  t_stress_symm(i, j) = t_D(i, j, k, l) * t_small_strain_symm(k, l);
1512  tensor_to_tensor(t_stress_symm, t_stress);
1513 
1514  const double psi = 0.5 * t_stress_symm(i, j) * t_small_strain_symm(i, j);
1515 
1516  CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &psi);
1517  CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
1518  &t_stress(0, 0));
1519 
1520  ++t_h;
1521  }
1522 
1524 }
1525 
1526 template <int S>
1527 HookeElement::OpAnalyticalInternalStrain_dx<S>::OpAnalyticalInternalStrain_dx(
1528  const std::string row_field,
1529  boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1530  StrainFunction strain_fun)
1531  : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1532  strainFun(strain_fun) {}
1533 
1534 template <int S>
1536 HookeElement::OpAnalyticalInternalStrain_dx<S>::iNtegrate(EntData &row_data) {
1542 
1543  auto get_tensor1 = [](VectorDouble &v, const int r) {
1545  &v(r + 0), &v(r + 1), &v(r + 2));
1546  };
1547 
1548  const int nb_integration_pts = getGaussPts().size2();
1549  auto t_coords = getFTensor1CoordsAtGaussPts();
1550 
1551  // get element volume
1552  double vol = getVolume();
1553  auto t_w = getFTensor0IntegrationWeight();
1554 
1555  nF.resize(nbRows, false);
1556  nF.clear();
1557 
1558  // elastic stiffness tensor (4th rank tensor with minor and major
1559  // symmetry)
1561  MAT_TO_DDG(dataAtPts->stiffnessMat));
1562 
1563  // get derivatives of base functions on rows
1564  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1565  const int row_nb_base_fun = row_data.getN().size2();
1566 
1567  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1568 
1569  auto t_fun_strain = strainFun(t_coords);
1571  t_stress(i, j) = -t_D(i, j, k, l) * t_fun_strain(k, l);
1572 
1573  // calculate scalar weight times element volume
1574  double a = t_w * vol;
1575 
1576  auto t_nf = get_tensor1(nF, 0);
1577 
1578  int rr = 0;
1579  for (; rr != nbRows / 3; ++rr) {
1580  t_nf(i) += a * t_row_diff_base(j) * t_stress(i, j);
1581  ++t_row_diff_base;
1582  ++t_nf;
1583  }
1584 
1585  for (; rr != row_nb_base_fun; ++rr)
1586  ++t_row_diff_base;
1587 
1588  ++t_w;
1589  ++t_coords;
1590  ++t_D;
1591  }
1592 
1594 }
1595 
1596 template <int S>
1597 HookeElement::OpAnalyticalInternalAleStrain_dX<S>::
1598  OpAnalyticalInternalAleStrain_dX(
1599  const std::string row_field,
1600  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1601  StrainFunction strain_fun,
1602  boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1603  : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1604  strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1605 
1606 template <int S>
1607 MoFEMErrorCode HookeElement::OpAnalyticalInternalAleStrain_dX<S>::iNtegrate(
1608  EntData &row_data) {
1614 
1615  auto get_tensor1 = [](VectorDouble &v, const int r) {
1617  &v(r + 0), &v(r + 1), &v(r + 2));
1618  };
1619 
1620  const int nb_integration_pts = getGaussPts().size2();
1621 
1622  auto get_coords = [&]() { return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
1623  auto t_coords = get_coords();
1624 
1625  // get element volume
1626  double vol = getVolume();
1627  auto t_w = getFTensor0IntegrationWeight();
1628 
1629  nF.resize(nbRows, false);
1630  nF.clear();
1631 
1632  // elastic stiffness tensor (4th rank tensor with minor and major
1633  // symmetry)
1635  MAT_TO_DDG(dataAtPts->stiffnessMat));
1636  auto t_F = getFTensor2FromMat<3, 3>(*(dataAtPts->FMat));
1637  auto &det_H = *dataAtPts->detHVec;
1638  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1639 
1640  // get derivatives of base functions on rows
1641  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1642  const int row_nb_base_fun = row_data.getN().size2();
1643 
1644  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1645 
1646  auto t_fun_strain = strainFun(t_coords);
1648  t_stress(i, j) = -t_D(i, j, k, l) * t_fun_strain(k, l);
1649  FTensor::Tensor2<double, 3, 3> t_eshelby_stress;
1650  t_eshelby_stress(i, j) = -t_F(k, i) * t_stress(k, j);
1651 
1652  // calculate scalar weight times element volume
1653  double a = t_w * vol * det_H[gg];
1654 
1655  auto t_nf = get_tensor1(nF, 0);
1656 
1657  int rr = 0;
1658  for (; rr != nbRows / 3; ++rr) {
1659  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1660  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1661  t_nf(i) += a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
1662  ++t_row_diff_base;
1663  ++t_nf;
1664  }
1665 
1666  for (; rr != row_nb_base_fun; ++rr)
1667  ++t_row_diff_base;
1668 
1669  ++t_w;
1670  ++t_coords;
1671  ++t_F;
1672  ++t_invH;
1673  ++t_D;
1674  }
1675 
1677 }
1678 
1679 template <int S>
1680 HookeElement::OpAnalyticalInternalAleStrain_dx<S>::
1681  OpAnalyticalInternalAleStrain_dx(
1682  const std::string row_field,
1683  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1684  StrainFunction strain_fun,
1685  boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1686  : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1687  strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1688 
1689 template <int S>
1690 MoFEMErrorCode HookeElement::OpAnalyticalInternalAleStrain_dx<S>::iNtegrate(
1691  EntData &row_data) {
1697 
1698  auto get_tensor1 = [](VectorDouble &v, const int r) {
1700  &v(r + 0), &v(r + 1), &v(r + 2));
1701  };
1702 
1703  const int nb_integration_pts = getGaussPts().size2();
1704 
1705  auto get_coords = [&]() { return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
1706  auto t_coords = get_coords();
1707 
1708  // get element volume
1709  double vol = getVolume();
1710  auto t_w = getFTensor0IntegrationWeight();
1711 
1712  nF.resize(nbRows, false);
1713  nF.clear();
1714 
1715  // elastic stiffness tensor (4th rank tensor with minor and major
1716  // symmetry)
1718  MAT_TO_DDG(dataAtPts->stiffnessMat));
1719  auto &det_H = *dataAtPts->detHVec;
1720  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1721 
1722  // get derivatives of base functions on rows
1723  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1724  const int row_nb_base_fun = row_data.getN().size2();
1725 
1726  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1727 
1728  auto t_fun_strain = strainFun(t_coords);
1730  t_stress(i, j) = -t_D(i, j, k, l) * t_fun_strain(k, l);
1731 
1732  // calculate scalar weight times element volume
1733  double a = t_w * vol * det_H[gg];
1734 
1735  auto t_nf = get_tensor1(nF, 0);
1736 
1737  int rr = 0;
1738  for (; rr != nbRows / 3; ++rr) {
1739  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1740  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1741  t_nf(i) += a * t_row_diff_base_pulled(j) * t_stress(i, j);
1742  ++t_row_diff_base;
1743  ++t_nf;
1744  }
1745 
1746  for (; rr != row_nb_base_fun; ++rr)
1747  ++t_row_diff_base;
1748 
1749  ++t_w;
1750  ++t_coords;
1751  ++t_invH;
1752  ++t_D;
1753  }
1754 
1756 }
1757 
1758 #endif // __HOOKE_ELEMENT_HPP
OpPostProcHookeElement::isFieldDisp
bool isFieldDisp
Definition: HookeElement.hpp:644
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
OpAssemble::rowIndices
VectorInt rowIndices
Definition: HookeElement.hpp:370
MoFEM::K
VectorDouble K
Definition: Projection10NodeCoordsOnField.cpp:125
setBlocks
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr)
FTensor::Christof
Definition: Christof_value.hpp:9
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
DataAtIntegrationPts::energyVec
boost::shared_ptr< VectorDouble > energyVec
Definition: HookeElement.hpp:91
OpAssemble::isDiag
bool isDiag
true if this block is on diagonal
Definition: HookeElement.hpp:376
OpAleLhsWithDensity_dX_dX::rhoGradAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
Definition: HookeElement.hpp:491
DataAtIntegrationPts::forcesOnlyOnEntitiesRow
Range forcesOnlyOnEntitiesRow
Definition: HookeElement.hpp:114
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1631
PlasticOps::w
double w(double eqiv, double dot_tau, double f, double sigma_y, double sigma_Y)
Definition: PlasticOpsGeneric.hpp:276
OpAnalyticalInternalStrain_dx::strainFun
StrainFunction strainFun
Definition: HookeElement.hpp:583
OpPostProcHookeElement::isALE
bool isALE
Definition: HookeElement.hpp:643
FTensor::Tensor1
Definition: Tensor1_value.hpp:8
EntityHandle
DataAtIntegrationPts::FMat
boost::shared_ptr< MatrixDouble > FMat
Definition: HookeElement.hpp:83
OpCalculateStiffnessScaledByDensityField::rHo0
const double rHo0
p_0 reference density in E(p) = E * (p / p_0)^n
Definition: HookeElement.hpp:327
OpLhs_dx_dx
Definition: HookeElement.hpp:410
OpCalculateStiffnessScaledByDensityField::rhoAtGaussPtsPtr
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
Definition: HookeElement.hpp:325
OpAleLhsPre_dX_dx
Definition: HookeElement.hpp:536
OpCalculateStrainAle::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:139
BlockData
struct ConvectiveMassElement BlockData
OpCalculateMassMatrix::massData
MassBlockData & massData
Definition: HookeElement.hpp:218
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
OpAssemble::transK
MatrixDouble transK
Definition: HookeElement.hpp:365
OpAssemble::nbCols
int nbCols
number if dof on column
Definition: HookeElement.hpp:374
OpPostProcHookeElement::mapGaussPts
std::vector< EntityHandle > & mapGaussPts
Definition: HookeElement.hpp:642
ELEMENT
OpCalculateEnergy::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:173
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
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:549
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
OpCalculateStiffnessScaledByDensityField::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:323
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:222
OpAleLhsWithDensity_dX_dX
Definition: HookeElement.hpp:488
OpCalculateHomogeneousStiffness::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:203
MoFEM::th
Tag th
Definition: Projection10NodeCoordsOnField.cpp:122
OpAnalyticalInternalAleStrain_dX::matPosAtPtsPtr
boost::shared_ptr< MatrixDouble > matPosAtPtsPtr
Definition: HookeElement.hpp:608
BasicFiniteElements.hpp
ConvectiveMassElement::BlockData::tEts
Range tEts
elements in block set
Definition: ConvectiveMassElement.hpp:119
OpAleLhsWithDensity_dx_dX::rHo0
const double rHo0
Definition: HookeElement.hpp:469
DataAtIntegrationPts::eshelbyStressMat
boost::shared_ptr< MatrixDouble > eshelbyStressMat
Definition: HookeElement.hpp:92
OpCalculateEnergy::ghostVec
SmartPetscObj< Vec > ghostVec
Definition: HookeElement.hpp:174
FTensor::Tensor2_symmetric
Definition: Tensor2_symmetric_value.hpp:13
OpAleLhsWithDensity_dx_dX
Definition: HookeElement.hpp:464
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, SmartPetscObj< Vec > &v_energy_ptr)
sdf.r
int r
Definition: sdf.py:8
FTensor::Tensor2< double *, 3, 3 >
DataAtIntegrationPts::hMat
boost::shared_ptr< MatrixDouble > hMat
Definition: HookeElement.hpp:82
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
OpAleLhsWithDensity_dX_dX::rhoAtGaussPtsPtr
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
Definition: HookeElement.hpp:490
NonlinearElasticElement::BlockData::forcesOnlyOnEntitiesRow
Range forcesOnlyOnEntitiesRow
Definition: HookeElement.hpp:37
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:1566
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
OpCalculateMassMatrix
Assemble mass matrix for elastic element TODO: CHANGE FORMULA *.
Definition: HookeElement.hpp:212
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1489
ConvectiveMassElement
structure grouping operators and data used for calculation of mass (convective) element
Definition: ConvectiveMassElement.hpp:28
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Tensor3
Definition: Tensor3_value.hpp:12
OpAleLhsPre_dX_dx::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:544
OpCalculateStrainAle
Definition: HookeElement.hpp:130
OpCalculateMassMatrix::doWork
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Definition: HookeElement.hpp:231
OpAleRhs_dx
Definition: HookeElement.hpp:425
OpCalculateHomogeneousStiffness
Definition: HookeElement.hpp:190
a
constexpr double a
Definition: approx_sphere.cpp:30
OpCalculateStrain::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:127
OpAssemble::nbIntegrationPts
int nbIntegrationPts
number of integration points
Definition: HookeElement.hpp:375
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
NonlinearElasticElement::BlockData::E
double E
Definition: HookeElement.hpp:34
OpAnalyticalInternalAleStrain_dX
Definition: HookeElement.hpp:586
convert.type
type
Definition: convert.py:64
OpAssemble::K
MatrixDouble K
Definition: HookeElement.hpp:364
OpAnalyticalInternalAleStrain_dX::strainFun
StrainFunction strainFun
Definition: HookeElement.hpp:607
OpAssemble::colIndices
VectorInt colIndices
Definition: HookeElement.hpp:371
OpAnalyticalInternalAleStrain_dX::StrainFunction
boost::function< FTensor::Tensor2_symmetric< double, 3 > FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords) > StrainFunction
Definition: HookeElement.hpp:597
OpAssemble::nbRows
int nbRows
number of dofs on rows
Definition: HookeElement.hpp:373
DataAtIntegrationPts::DataAtIntegrationPts
DataAtIntegrationPts()
Definition: HookeElement.hpp:96
OpAleRhs_dX
Definition: HookeElement.hpp:512
OpPostProcHookeElement::blockSetsPtr
map< int, BlockData > & blockSetsPtr
Definition: HookeElement.hpp:640
OpCalculateEshelbyStress
Definition: HookeElement.hpp:177
HookeElement
OpCalculateEnergy
Definition: HookeElement.hpp:164
OpPostProcHookeElement::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:638
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1201
OpAnalyticalInternalStrain_dx
Definition: HookeElement.hpp:563
OpCalculateMassMatrix::locK
MatrixDouble locK
Definition: HookeElement.hpp:215
OpCalculateStiffnessScaledByDensityField::rhoN
const double rhoN
exponent n in E(p) = E * (p / p_0)^n
Definition: HookeElement.hpp:326
FTensor::PackPtr
Definition: FTensor.hpp:54
OpAnalyticalInternalAleStrain_dx
Definition: HookeElement.hpp:611
ConvectiveMassElement::BlockData
data for calculation inertia forces
Definition: ConvectiveMassElement.hpp:116
OpAnalyticalInternalAleStrain_dx::matPosAtPtsPtr
boost::shared_ptr< MatrixDouble > matPosAtPtsPtr
Definition: HookeElement.hpp:633
OpAleLhsWithDensity_dX_dX::rhoN
const double rhoN
Definition: HookeElement.hpp:492
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1876
NonlinearElasticElement::BlockData
data for calculation heat conductivity and heat capacity elements
Definition: HookeElement.hpp:32
FTensor::Tensor4
Definition: Tensor4_value.hpp:18
OpCalculateEshelbyStress::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:186
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
OpAnalyticalInternalAleStrain_dx::StrainFunction
boost::function< FTensor::Tensor2_symmetric< double, 3 > FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords) > StrainFunction
Definition: HookeElement.hpp:622
NonlinearElasticElement::BlockData::PoissonRatio
double PoissonRatio
Definition: HookeElement.hpp:35
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
EntData
EntitiesFieldData::EntData EntData
Definition: HookeElement.hpp:74
FTensor::Index< 'i', 3 >
DataAtIntegrationPts::invHMat
boost::shared_ptr< MatrixDouble > invHMat
Definition: HookeElement.hpp:87
ConvectiveMassElement::BlockData::rho0
double rho0
reference density
Definition: ConvectiveMassElement.hpp:117
convert.n
n
Definition: convert.py:82
OpAssemble::nF
VectorDouble nF
Definition: HookeElement.hpp:366
MoFEM::determinantTensor3by3
static auto determinantTensor3by3(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:1518
DataAtIntegrationPts::HMat
boost::shared_ptr< MatrixDouble > HMat
Definition: HookeElement.hpp:85
OpAnalyticalInternalStrain_dx::StrainFunction
boost::function< FTensor::Tensor2_symmetric< double, 3 > FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > &t_coords) > StrainFunction
Definition: HookeElement.hpp:574
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
Range
NonlinearElasticElement::BlockData::forcesOnlyOnEntitiesCol
Range forcesOnlyOnEntitiesCol
Definition: HookeElement.hpp:38
NonlinearElasticElement::BlockData::iD
int iD
Definition: HookeElement.hpp:33
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:526
MoFEM::EntitiesFieldData::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: EntitiesFieldData.hpp:1305
OpCalculateStrain
Definition: HookeElement.hpp:119
OpCalculateStiffnessScaledByDensityField::blockSetsPtr
boost::shared_ptr< map< int, BlockData > > blockSetsPtr
Definition: HookeElement.hpp:321
MAT_TO_DDG
#define MAT_TO_DDG(SM)
Definition: HookeElement.hpp:142
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
MoFEM::Types::VectorInt
UBlasVector< int > VectorInt
Definition: Types.hpp:67
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
DataAtIntegrationPts::eshelbyStress_dx
boost::shared_ptr< MatrixDouble > eshelbyStress_dx
Definition: HookeElement.hpp:94
FTensor::Ddg
Definition: Ddg_value.hpp:7
OpAleLhsWithDensity_dX_dX::rHo0
const double rHo0
Definition: HookeElement.hpp:493
NonlinearElasticElement::BlockData::tEts
Range tEts
constrains elements in block set
Definition: HookeElement.hpp:36
OpCalculateMassMatrix::dAta
BlockData & dAta
Definition: HookeElement.hpp:217
OpCalculateStress::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:161
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
OpAleLhs_dx_dX
Definition: HookeElement.hpp:449
MoFEM::BlockData
Definition: MeshsetsManager.cpp:752
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)
Definition: HookeElement.cpp:533
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
OpAnalyticalInternalAleStrain_dx::strainFun
StrainFunction strainFun
Definition: HookeElement.hpp:632
OpCalculateStress
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, boost::shared_ptr< DataAtIntegrationPts > data_at_pts=nullptr)
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
DataAtIntegrationPts::forcesOnlyOnEntitiesCol
Range forcesOnlyOnEntitiesCol
Definition: HookeElement.hpp:115
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
OpAssemble
Definition: HookeElement.hpp:340
OpPostProcHookeElement::postProcMesh
moab::Interface & postProcMesh
Definition: HookeElement.hpp:641
OpAleLhsWithDensity_dx_dX::rhoAtGaussPtsPtr
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
Definition: HookeElement.hpp:466
OpPostProcHookeElement
Definition: HookeElement.hpp:637
OpAleLhs_dx_dx
Definition: HookeElement.hpp:434
MoFEM::SmartPetscObj< Vec >
OpAleLhs_dX_dX
Definition: HookeElement.hpp:521
OpCalculateHomogeneousStiffness::blockSetsPtr
boost::shared_ptr< map< int, BlockData > > blockSetsPtr
Definition: HookeElement.hpp:201
invJac
MatrixDouble invJac
Definition: HookeElement.hpp:683
OpAssemble::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:368
k
FTensor::Index< 'k', 3 > k
Definition: matrix_function.cpp:20
OpAleLhsWithDensity_dx_dX::rhoGradAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
Definition: HookeElement.hpp:467
OpRhs_dx
Definition: HookeElement.hpp:401
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
OpAleLhsWithDensity_dx_dX::rhoN
const double rhoN
Definition: HookeElement.hpp:468
DataAtIntegrationPts::smallStrainMat
boost::shared_ptr< MatrixDouble > smallStrainMat
Definition: HookeElement.hpp:81
DataAtIntegrationPts::stiffnessMat
boost::shared_ptr< MatrixDouble > stiffnessMat
Definition: HookeElement.hpp:90
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
OpCalculateMassMatrix::translocK
MatrixDouble translocK
Definition: HookeElement.hpp:216
OpCalculateMassMatrix::commonData
boost::shared_ptr< DataAtIntegrationPts > commonData
Definition: HookeElement.hpp:220
OpCalculateStiffnessScaledByDensityField
Definition: HookeElement.hpp:318
ConvectiveMassElement::BlockData::a0
VectorDouble a0
constant acceleration
Definition: ConvectiveMassElement.hpp:118
DataAtIntegrationPts::detHVec
boost::shared_ptr< VectorDouble > detHVec
Definition: HookeElement.hpp:86
OpAleLhs_dX_dx
Definition: HookeElement.hpp:547
DataAtIntegrationPts::cauchyStressMat
boost::shared_ptr< MatrixDouble > cauchyStressMat
Definition: HookeElement.hpp:89