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