v0.9.2
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 
661  OpPostProcHookeElement(const string row_field,
662  boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
663  map<int, BlockData> &block_sets_ptr,
664  moab::Interface &post_proc_mesh,
665  std::vector<EntityHandle> &map_gauss_pts,
666  bool is_ale = false);
667 
668  MoFEMErrorCode doWork(int side, EntityType type,
669  DataForcesAndSourcesCore::EntData &data);
670  };
671 
672  static MoFEMErrorCode
673  setBlocks(MoFEM::Interface &m_field,
674  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr);
675 
676  static MoFEMErrorCode
677  addElasticElement(MoFEM::Interface &m_field,
678  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
679  const std::string element_name, const std::string x_field,
680  const std::string X_field, const bool ale);
681 
682  static MoFEMErrorCode
683  setOperators(boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr,
684  boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr,
685  boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
686  const std::string x_field, const std::string X_field,
687  const bool ale, const bool field_disp,
688  const EntityType type = MBTET,
689  boost::shared_ptr<DataAtIntegrationPts> data_at_pts = nullptr);
690 
691  static MoFEMErrorCode
692  calculateEnergy(DM dm, boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
693  const std::string x_field, const std::string X_field,
694  const bool ale, const bool field_disp, Vec *v_energy_ptr);
695 
696 private:
697  MatrixDouble invJac;
698 };
699 
700 template <bool D>
701 HookeElement::OpCalculateStrain<D>::OpCalculateStrain(
702  const std::string row_field, const std::string col_field,
703  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
704  : VolUserDataOperator(row_field, col_field, OPROW, true),
705  dataAtPts(data_at_pts) {
706  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
707 }
708 
709 template <bool D>
710 MoFEMErrorCode HookeElement::OpCalculateStrain<D>::doWork(int row_side,
711  EntityType row_type,
712  EntData &row_data) {
713  MoFEMFunctionBegin;
714  FTensor::Index<'i', 3> i;
715  FTensor::Index<'j', 3> j;
716  // get number of integration points
717  const int nb_integration_pts = getGaussPts().size2();
718  dataAtPts->smallStrainMat->resize(6, nb_integration_pts, false);
719  auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
720  auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hMat));
721 
722  for (int gg = 0; gg != nb_integration_pts; ++gg) {
723  t_strain(i, j) = (t_h(i, j) || t_h(j, i)) / 2.;
724 
725  // If displacement field, not field o spatial positons is given
726  if (!D) {
727  t_strain(0, 0) -= 1;
728  t_strain(1, 1) -= 1;
729  t_strain(2, 2) -= 1;
730  }
731 
732  ++t_strain;
733  ++t_h;
734  }
735  MoFEMFunctionReturn(0);
736 }
737 
738 template <int S>
739 HookeElement::OpAleLhs_dx_dx<S>::OpAleLhs_dx_dx(
740  const std::string row_field, const std::string col_field,
741  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
742  : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
743 
744 template <int S>
745 MoFEMErrorCode HookeElement::OpAleLhs_dx_dx<S>::iNtegrate(EntData &row_data,
746  EntData &col_data) {
747  MoFEMFunctionBegin;
748 
749  // get sub-block (3x3) of local stiffens matrix, here represented by
750  // second order tensor
751  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
752  return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
753  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
754  &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
755  &m(r + 2, c + 2));
756  };
757 
758  FTensor::Index<'i', 3> i;
759  FTensor::Index<'j', 3> j;
760  FTensor::Index<'k', 3> k;
761  FTensor::Index<'l', 3> l;
762 
763  // get element volume
764  double vol = getVolume();
765 
766  // get intergrayion weights
767  auto t_w = getFTensor0IntegrationWeight();
768 
769  // get derivatives of base functions on rows
770  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
771  const int row_nb_base_fun = row_data.getN().size2();
772 
773  // Elastic stiffness tensor (4th rank tensor with minor and major
774  // symmetry)
775  FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
776  MAT_TO_DDG(dataAtPts->stiffnessMat));
777 
778  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
779  auto &det_H = *dataAtPts->detHVec;
780 
781  // iterate over integration points
782  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
783 
784  // calculate scalar weight times element volume
785  double a = t_w * vol * det_H[gg];
786 
787  // iterate over row base functions
788  int rr = 0;
789  for (; rr != nbRows / 3; ++rr) {
790 
791  // get sub matrix for the row
792  auto t_m = get_tensor2(K, 3 * rr, 0);
793 
794  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
795  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
796 
797  FTensor::Christof<double, 3, 3> t_rowD;
798  // I mix up the indices here so that it behaves like a
799  // Dg. That way I don't have to have a separate wrapper
800  // class Christof_Expr, which simplifies things.
801  t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base_pulled(i));
802 
803  // get derivatives of base functions for columns
804  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
805 
806  // iterate column base functions
807  for (int cc = 0; cc != nbCols / 3; ++cc) {
808 
809  FTensor::Tensor1<double, 3> t_col_diff_base_pulled;
810  t_col_diff_base_pulled(j) = t_col_diff_base(i) * t_invH(i, j);
811 
812  // integrate block local stiffens matrix
813  t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base_pulled(k);
814 
815  // move to next column base function
816  ++t_col_diff_base;
817 
818  // move to next block of local stiffens matrix
819  ++t_m;
820  }
821 
822  // move to next row base function
823  ++t_row_diff_base;
824  }
825 
826  for (; rr != row_nb_base_fun; ++rr)
827  ++t_row_diff_base;
828 
829  // move to next integration weight
830  ++t_w;
831  ++t_D;
832  ++t_invH;
833  }
834 
835  MoFEMFunctionReturn(0);
836 }
837 
838 template <int S>
839 HookeElement::OpCalculateHomogeneousStiffness<S>::
840  OpCalculateHomogeneousStiffness(
841  const std::string row_field, const std::string col_field,
842  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
843  boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
844  : VolUserDataOperator(row_field, col_field, OPROW, true),
845  blockSetsPtr(block_sets_ptr), dataAtPts(data_at_pts) {
846  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
847 }
848 
849 template <int S>
850 MoFEMErrorCode HookeElement::OpCalculateHomogeneousStiffness<S>::doWork(
851  int row_side, EntityType row_type, EntData &row_data) {
852  MoFEMFunctionBegin;
853 
854  for (auto &m : (*blockSetsPtr)) {
855  if (m.second.tEts.find(getFEEntityHandle()) != m.second.tEts.end()) {
856 
857  dataAtPts->stiffnessMat->resize(36, 1, false);
858  FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
859  MAT_TO_DDG(dataAtPts->stiffnessMat));
860  const double young = m.second.E;
861  const double poisson = m.second.PoissonRatio;
862 
863  // coefficient used in intermediate calculation
864  const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
865 
866  FTensor::Index<'i', 3> i;
867  FTensor::Index<'j', 3> j;
868  FTensor::Index<'k', 3> k;
869  FTensor::Index<'l', 3> l;
870 
871  t_D(i, j, k, l) = 0.;
872 
873  t_D(0, 0, 0, 0) = 1 - poisson;
874  t_D(1, 1, 1, 1) = 1 - poisson;
875  t_D(2, 2, 2, 2) = 1 - poisson;
876 
877  t_D(0, 1, 0, 1) = 0.5 * (1 - 2 * poisson);
878  t_D(0, 2, 0, 2) = 0.5 * (1 - 2 * poisson);
879  t_D(1, 2, 1, 2) = 0.5 * (1 - 2 * poisson);
880 
881  t_D(0, 0, 1, 1) = poisson;
882  t_D(1, 1, 0, 0) = poisson;
883  t_D(0, 0, 2, 2) = poisson;
884  t_D(2, 2, 0, 0) = poisson;
885  t_D(1, 1, 2, 2) = poisson;
886  t_D(2, 2, 1, 1) = poisson;
887 
888  t_D(i, j, k, l) *= coefficient;
889 
890  break;
891  }
892  }
893 
894  MoFEMFunctionReturn(0);
895 }
896 
897 template <int S>
898 HookeElement::OpCalculateStress<S>::OpCalculateStress(
899  const std::string row_field, const std::string col_field,
900  boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
901  : VolUserDataOperator(row_field, col_field, OPROW, true),
902  dataAtPts(data_at_pts) {
903  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
904 }
905 
906 template <int S>
907 MoFEMErrorCode HookeElement::OpCalculateStress<S>::doWork(int row_side,
908  EntityType row_type,
909  EntData &row_data) {
910  MoFEMFunctionBegin;
911  // get number of integration points
912  const int nb_integration_pts = getGaussPts().size2();
913  auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
914  dataAtPts->cauchyStressMat->resize(6, nb_integration_pts, false);
915  auto t_cauchy_stress =
916  getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
917 
918  FTensor::Index<'i', 3> i;
919  FTensor::Index<'j', 3> j;
920  FTensor::Index<'k', 3> k;
921  FTensor::Index<'l', 3> l;
922 
923  // elastic stiffness tensor (4th rank tensor with minor and major
924  // symmetry)
925  FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
926  MAT_TO_DDG(dataAtPts->stiffnessMat));
927  for (int gg = 0; gg != nb_integration_pts; ++gg) {
928  t_cauchy_stress(i, j) = t_D(i, j, k, l) * t_strain(k, l);
929  ++t_strain;
930  ++t_cauchy_stress;
931  ++t_D;
932  }
933  MoFEMFunctionReturn(0);
934 }
935 
936 template <int S>
937 HookeElement::OpLhs_dx_dx<S>::OpLhs_dx_dx(
938  const std::string row_field, const std::string col_field,
939  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
940  : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
941 
942 template <int S>
943 MoFEMErrorCode HookeElement::OpLhs_dx_dx<S>::iNtegrate(EntData &row_data,
944  EntData &col_data) {
945  MoFEMFunctionBegin;
946 
947  // get sub-block (3x3) of local stiffens matrix, here represented by
948  // second order tensor
949  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
950  return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
951  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
952  &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
953  &m(r + 2, c + 2));
954  };
955 
956  FTensor::Index<'i', 3> i;
957  FTensor::Index<'j', 3> j;
958  FTensor::Index<'k', 3> k;
959  FTensor::Index<'l', 3> l;
960 
961  // get element volume
962  double vol = getVolume();
963 
964  // get intergrayion weights
965  auto t_w = getFTensor0IntegrationWeight();
966 
967  // get derivatives of base functions on rows
968  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
969  const int row_nb_base_fun = row_data.getN().size2();
970 
971  // Elastic stiffness tensor (4th rank tensor with minor and major
972  // symmetry)
973  FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
974  MAT_TO_DDG(dataAtPts->stiffnessMat));
975 
976  // iterate over integration points
977  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
978 
979  // calculate scalar weight times element volume
980  double a = t_w * vol;
981  if (getHoGaussPtsDetJac().size()) {
982  a *= getHoGaussPtsDetJac()[gg];
983  }
984 
985  // iterate over row base functions
986  int rr = 0;
987  for (; rr != nbRows / 3; ++rr) {
988 
989  // get sub matrix for the row
990  auto t_m = get_tensor2(K, 3 * rr, 0);
991 
992  // get derivatives of base functions for columns
993  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
994 
995  FTensor::Christof<double, 3, 3> t_rowD;
996  // I mix up the indices here so that it behaves like a
997  // Dg. That way I don't have to have a separate wrapper
998  // class Christof_Expr, which simplifies things.
999  t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base(i));
1000 
1001  // iterate column base functions
1002  for (int cc = 0; cc != nbCols / 3; ++cc) {
1003 
1004  // integrate block local stiffens matrix
1005  t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base(k);
1006 
1007  // move to next column base function
1008  ++t_col_diff_base;
1009 
1010  // move to next block of local stiffens matrix
1011  ++t_m;
1012  }
1013 
1014  // move to next row base function
1015  ++t_row_diff_base;
1016  }
1017 
1018  for (; rr != row_nb_base_fun; ++rr)
1019  ++t_row_diff_base;
1020 
1021  // move to next integration weight
1022  ++t_w;
1023  ++t_D;
1024  }
1025 
1026  MoFEMFunctionReturn(0);
1027 }
1028 
1029 template <int S>
1030 HookeElement::OpAleLhs_dx_dX<S>::OpAleLhs_dx_dX(
1031  const std::string row_field, const std::string col_field,
1032  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1033  : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
1034 
1035 template <int S>
1036 MoFEMErrorCode HookeElement::OpAleLhs_dx_dX<S>::iNtegrate(EntData &row_data,
1037  EntData &col_data) {
1038  MoFEMFunctionBegin;
1039 
1040  // get sub-block (3x3) of local stiffens matrix, here represented by
1041  // second order tensor
1042  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1043  return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
1044  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1045  &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1046  &m(r + 2, c + 2));
1047  };
1048 
1049  FTensor::Index<'i', 3> i;
1050  FTensor::Index<'j', 3> j;
1051  FTensor::Index<'k', 3> k;
1052  FTensor::Index<'l', 3> l;
1053  FTensor::Index<'m', 3> m;
1054  FTensor::Index<'n', 3> n;
1055 
1056  // get element volume
1057  double vol = getVolume();
1058 
1059  // get intergrayion weights
1060  auto t_w = getFTensor0IntegrationWeight();
1061 
1062  // get derivatives of base functions on rows
1063  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1064  const int row_nb_base_fun = row_data.getN().size2();
1065 
1066  // Elastic stiffness tensor (4th rank tensor with minor and major
1067  // symmetry)
1068  FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
1069  MAT_TO_DDG(dataAtPts->stiffnessMat));
1070 
1071  auto t_cauchy_stress =
1072  getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1073  auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1074  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1075  auto &det_H = *dataAtPts->detHVec;
1076 
1077  // iterate over integration points
1078  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1079 
1080  // calculate scalar weight times element volume
1081  double a = t_w * vol * det_H[gg];
1082 
1083  FTensor::Tensor4<double, 3, 3, 3, 3> t_F_dX;
1084  t_F_dX(i, j, k, l) = -(t_h(i, m) * t_invH(m, k)) * t_invH(l, j);
1085 
1086  // iterate over row base functions
1087  int rr = 0;
1088  for (; rr != nbRows / 3; ++rr) {
1089 
1090  // get sub matrix for the row
1091  auto t_m = get_tensor2(K, 3 * rr, 0);
1092 
1093  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1094  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1095 
1096  FTensor::Tensor1<double, 3> t_row_stress;
1097  t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_cauchy_stress(i, j);
1098 
1099  FTensor::Tensor3<double, 3, 3, 3> t_row_diff_base_pulled_dX;
1100  t_row_diff_base_pulled_dX(j, k, l) =
1101  -(t_invH(i, k) * t_row_diff_base(i)) * t_invH(l, j);
1102 
1103  FTensor::Tensor3<double, 3, 3, 3> t_row_dX_stress;
1104  t_row_dX_stress(i, k, l) =
1105  a * (t_row_diff_base_pulled_dX(j, k, l) * t_cauchy_stress(j, i));
1106 
1107  FTensor::Christof<double, 3, 3> t_row_D;
1108  t_row_D(l, j, k) = (a * t_row_diff_base_pulled(i)) * t_D(i, j, k, l);
1109 
1110  FTensor::Tensor3<double, 3, 3, 3> t_row_stress_dX;
1111  // FIXME: This operator is not implemented, doing operation by hand
1112  // t_row_stress_dX(i, m, n) = t_row_D(i, k, l) * t_F_dX(k, l, m, n);
1113  t_row_stress_dX(i, j, k) = 0;
1114  for (int ii = 0; ii != 3; ++ii)
1115  for (int mm = 0; mm != 3; ++mm)
1116  for (int nn = 0; nn != 3; ++nn) {
1117  auto &v = t_row_stress_dX(ii, mm, nn);
1118  for (int kk = 0; kk != 3; ++kk)
1119  for (int ll = 0; ll != 3; ++ll)
1120  v += t_row_D(ii, kk, ll) * t_F_dX(kk, ll, mm, nn);
1121  }
1122 
1123  // get derivatives of base functions for columns
1124  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1125 
1126  // iterate column base functions
1127  for (int cc = 0; cc != nbCols / 3; ++cc) {
1128 
1129  t_m(i, k) += t_row_stress(i) * (t_invH(j, k) * t_col_diff_base(j));
1130  t_m(i, k) += t_row_dX_stress(i, k, l) * t_col_diff_base(l);
1131  t_m(i, k) += t_row_stress_dX(i, k, l) * t_col_diff_base(l);
1132 
1133  // move to next column base function
1134  ++t_col_diff_base;
1135 
1136  // move to next block of local stiffens matrix
1137  ++t_m;
1138  }
1139 
1140  // move to next row base function
1141  ++t_row_diff_base;
1142  }
1143 
1144  for (; rr != row_nb_base_fun; ++rr)
1145  ++t_row_diff_base;
1146 
1147  // move to next integration weight
1148  ++t_w;
1149  ++t_D;
1150  ++t_cauchy_stress;
1151  ++t_invH;
1152  ++t_h;
1153  }
1154 
1155  MoFEMFunctionReturn(0);
1156 }
1157 
1158 template <int S>
1159 HookeElement::OpAleLhs_dX_dX<S>::OpAleLhs_dX_dX(
1160  const std::string row_field, const std::string col_field,
1161  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1162  : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
1163 
1164 template <int S>
1165 MoFEMErrorCode HookeElement::OpAleLhs_dX_dX<S>::iNtegrate(EntData &row_data,
1166  EntData &col_data) {
1167  MoFEMFunctionBegin;
1168 
1169  // get sub-block (3x3) of local stiffens matrix, here represented by
1170  // second order tensor
1171  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1172  return FTensor::Tensor2<FTensor::PackPtr<double *, 3>, 3, 3>(
1173  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1174  &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1175  &m(r + 2, c + 2));
1176  };
1177 
1178  FTensor::Index<'i', 3> i;
1179  FTensor::Index<'j', 3> j;
1180  FTensor::Index<'k', 3> k;
1181  FTensor::Index<'l', 3> l;
1182  FTensor::Index<'m', 3> m;
1183  FTensor::Index<'n', 3> n;
1184 
1185  // get element volume
1186  double vol = getVolume();
1187 
1188  // get intergrayion weights
1189  auto t_w = getFTensor0IntegrationWeight();
1190 
1191  // get derivatives of base functions on rows
1192  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1193  const int row_nb_base_fun = row_data.getN().size2();
1194 
1195  // Elastic stiffness tensor (4th rank tensor with minor and major
1196  // symmetry)
1197  FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
1198  MAT_TO_DDG(dataAtPts->stiffnessMat));
1199  auto t_cauchy_stress =
1200  getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1201  auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
1202  auto t_eshelby_stress =
1203  getFTensor2FromMat<3, 3>(*dataAtPts->eshelbyStressMat);
1204  auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1205  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1206  auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
1207  auto &det_H = *dataAtPts->detHVec;
1208 
1209  // iterate over integration points
1210  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1211 
1212  // calculate scalar weight times element volume
1213  double a = t_w * vol * det_H[gg];
1214 
1215  FTensor::Tensor4<double, 3, 3, 3, 3> t_F_dX;
1216  t_F_dX(i, j, k, l) = -(t_h(i, m) * t_invH(m, k)) * t_invH(l, j);
1217 
1218  FTensor::Tensor4<double, 3, 3, 3, 3> t_D_strain_dX;
1219  t_D_strain_dX(i, j, m, n) = 0.;
1220  for (int ii = 0; ii != 3; ++ii)
1221  for (int jj = 0; jj != 3; ++jj)
1222  for (int ll = 0; ll != 3; ++ll)
1223  for (int kk = 0; kk != 3; ++kk) {
1224  auto &v = t_D_strain_dX(ii, jj, kk, ll);
1225  for (int mm = 0; mm != 3; ++mm)
1226  for (int nn = 0; nn != 3; ++nn)
1227  v += t_D(ii, jj, mm, nn) * t_F_dX(mm, nn, kk, ll);
1228  }
1229 
1230  FTensor::Tensor4<double, 3, 3, 3, 3> t_eshelby_stress_dX;
1231  t_eshelby_stress_dX(i, j, m, n) = t_F(k, i) * t_D_strain_dX(k, j, m, n);
1232 
1233  for (int ii = 0; ii != 3; ++ii)
1234  for (int jj = 0; jj != 3; ++jj)
1235  for (int mm = 0; mm != 3; ++mm)
1236  for (int nn = 0; nn != 3; ++nn) {
1237  auto &v = t_eshelby_stress_dX(ii, jj, mm, nn);
1238  for (int kk = 0; kk != 3; ++kk)
1239  v += t_F_dX(kk, ii, mm, nn) * t_cauchy_stress(kk, jj);
1240  }
1241 
1242  t_eshelby_stress_dX(i, j, k, l) *= -1;
1243 
1244  FTensor::Tensor2<double, 3, 3> t_energy_dX;
1245  t_energy_dX(k, l) = t_F_dX(i, j, k, l) * t_cauchy_stress(i, j);
1246  t_energy_dX(k, l) +=
1247  (t_strain(m, n) * t_D(m, n, i, j)) * t_F_dX(i, j, k, l);
1248  t_energy_dX(k, l) /= 2.;
1249 
1250  for (int kk = 0; kk != 3; ++kk)
1251  for (int ll = 0; ll != 3; ++ll) {
1252  auto v = t_energy_dX(kk, ll);
1253  for (int ii = 0; ii != 3; ++ii)
1254  t_eshelby_stress_dX(ii, ii, kk, ll) += v;
1255  }
1256 
1257  // iterate over row base functions
1258  int rr = 0;
1259  for (; rr != nbRows / 3; ++rr) {
1260 
1261  // get sub matrix for the row
1262  auto t_m = get_tensor2(K, 3 * rr, 0);
1263 
1264  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1265  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1266 
1267  FTensor::Tensor1<double, 3> t_row_stress;
1268  t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
1269 
1270  FTensor::Tensor3<double, 3, 3, 3> t_row_diff_base_pulled_dX;
1271  t_row_diff_base_pulled_dX(j, k, l) =
1272  -(t_row_diff_base(i) * t_invH(i, k)) * t_invH(l, j);
1273 
1274  FTensor::Tensor3<double, 3, 3, 3> t_row_dX_stress;
1275  t_row_dX_stress(i, k, l) =
1276  a * (t_row_diff_base_pulled_dX(j, k, l) * t_eshelby_stress(i, j));
1277 
1278  FTensor::Tensor3<double, 3, 3, 3> t_row_stress_dX;
1279  t_row_stress_dX(i, m, n) =
1280  a * t_row_diff_base_pulled(j) * t_eshelby_stress_dX(i, j, m, n);
1281 
1282  // get derivatives of base functions for columns
1283  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1284 
1285  // iterate column base functions
1286  for (int cc = 0; cc != nbCols / 3; ++cc) {
1287 
1288  t_m(i, k) += t_row_stress(i) * (t_invH(j, k) * t_col_diff_base(j));
1289  t_m(i, k) += t_row_dX_stress(i, k, l) * t_col_diff_base(l);
1290  t_m(i, k) += t_row_stress_dX(i, k, l) * t_col_diff_base(l);
1291 
1292  // move to next column base function
1293  ++t_col_diff_base;
1294 
1295  // move to next block of local stiffens matrix
1296  ++t_m;
1297  }
1298 
1299  // move to next row base function
1300  ++t_row_diff_base;
1301  }
1302 
1303  for (; rr != row_nb_base_fun; ++rr)
1304  ++t_row_diff_base;
1305 
1306  // move to next integration weight
1307  ++t_w;
1308  ++t_D;
1309  ++t_cauchy_stress;
1310  ++t_strain;
1311  ++t_eshelby_stress;
1312  ++t_h;
1313  ++t_invH;
1314  ++t_F;
1315  }
1316 
1317  MoFEMFunctionReturn(0);
1318 }
1319 
1320 template <int S>
1321 HookeElement::OpAleLhsPre_dX_dx<S>::OpAleLhsPre_dX_dx(
1322  const std::string row_field, const std::string col_field,
1323  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1324  : VolUserDataOperator(row_field, col_field, OPROW, true),
1325  dataAtPts(data_at_pts) {
1326  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1327 }
1328 
1329 template <int S>
1330 MoFEMErrorCode HookeElement::OpAleLhsPre_dX_dx<S>::doWork(int row_side,
1331  EntityType row_type,
1332  EntData &row_data) {
1333  MoFEMFunctionBegin;
1334 
1335  const int nb_integration_pts = row_data.getN().size1();
1336 
1337  auto get_eshelby_stress_dx = [this, nb_integration_pts]() {
1338  FTensor::Tensor4<FTensor::PackPtr<double *, 1>, 3, 3, 3, 3>
1339  t_eshelby_stress_dx;
1340  dataAtPts->eshelbyStress_dx->resize(81, nb_integration_pts, false);
1341  int mm = 0;
1342  for (int ii = 0; ii != 3; ++ii)
1343  for (int jj = 0; jj != 3; ++jj)
1344  for (int kk = 0; kk != 3; ++kk)
1345  for (int ll = 0; ll != 3; ++ll)
1346  t_eshelby_stress_dx.ptr(ii, jj, kk, ll) =
1347  &(*dataAtPts->eshelbyStress_dx)(mm++, 0);
1348  return t_eshelby_stress_dx;
1349  };
1350 
1351  auto t_eshelby_stress_dx = get_eshelby_stress_dx();
1352 
1353  FTensor::Index<'i', 3> i;
1354  FTensor::Index<'j', 3> j;
1355  FTensor::Index<'k', 3> k;
1356  FTensor::Index<'l', 3> l;
1357  FTensor::Index<'m', 3> m;
1358  FTensor::Index<'n', 3> n;
1359 
1360  // Elastic stiffness tensor (4th rank tensor with minor and major
1361  // symmetry)
1362  FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
1363  MAT_TO_DDG(dataAtPts->stiffnessMat));
1364  auto t_cauchy_stress =
1365  getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1366  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1367  auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
1368 
1369  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1370 
1371  t_eshelby_stress_dx(i, j, m, n) =
1372  (t_F(k, i) * t_D(k, j, m, l)) * t_invH(n, l);
1373  for (int ii = 0; ii != 3; ++ii)
1374  for (int jj = 0; jj != 3; ++jj)
1375  for (int mm = 0; mm != 3; ++mm)
1376  for (int nn = 0; nn != 3; ++nn) {
1377  auto &v = t_eshelby_stress_dx(ii, jj, mm, nn);
1378  v += t_invH(nn, ii) * t_cauchy_stress(mm, jj);
1379  }
1380  t_eshelby_stress_dx(i, j, k, l) *= -1;
1381 
1382  FTensor::Tensor2<double, 3, 3> t_energy_dx;
1383  t_energy_dx(m, n) = t_invH(n, j) * t_cauchy_stress(m, j);
1384 
1385  for (int mm = 0; mm != 3; ++mm)
1386  for (int nn = 0; nn != 3; ++nn) {
1387  auto v = t_energy_dx(mm, nn);
1388  for (int ii = 0; ii != 3; ++ii)
1389  t_eshelby_stress_dx(ii, ii, mm, nn) += v;
1390  }
1391 
1392  ++t_D;
1393  ++t_invH;
1394  ++t_cauchy_stress;
1395  ++t_eshelby_stress_dx;
1396  ++t_F;
1397  }
1398 
1399  MoFEMFunctionReturn(0);
1400 }
1401 
1402 template <class ELEMENT>
1403 HookeElement::OpPostProcHookeElement<ELEMENT>::OpPostProcHookeElement(
1404  const string row_field, boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1405  map<int, BlockData> &block_sets_ptr, moab::Interface &post_proc_mesh,
1406  std::vector<EntityHandle> &map_gauss_pts, bool is_ale)
1407  : ELEMENT::UserDataOperator(row_field, UserDataOperator::OPROW),
1408  dataAtPts(data_at_pts), blockSetsPtr(block_sets_ptr),
1409  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts), isALE(is_ale) {}
1410 
1411 template <class ELEMENT>
1412 MoFEMErrorCode HookeElement::OpPostProcHookeElement<ELEMENT>::doWork(
1413  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1414  MoFEMFunctionBegin;
1415 
1416  if (type != MBVERTEX) {
1417  MoFEMFunctionReturnHot(0);
1418  }
1419 
1420  auto tensor_to_tensor = [](const auto &t1, auto &t2) {
1421  t2(0, 0) = t1(0, 0);
1422  t2(1, 1) = t1(1, 1);
1423  t2(2, 2) = t1(2, 2);
1424  t2(0, 1) = t2(1, 0) = t1(1, 0);
1425  t2(0, 2) = t2(2, 0) = t1(2, 0);
1426  t2(1, 2) = t2(2, 1) = t1(2, 1);
1427  };
1428 
1429  double def_VAL[9];
1430  bzero(def_VAL, 9 * sizeof(double));
1431 
1432  Tag th_stress;
1433  CHKERR postProcMesh.tag_get_handle("STRESS", 9, MB_TYPE_DOUBLE, th_stress,
1434  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
1435  Tag th_psi;
1436  CHKERR postProcMesh.tag_get_handle("ENERGY", 1, MB_TYPE_DOUBLE, th_psi,
1437  MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
1438 
1439  const int nb_integration_pts = mapGaussPts.size();
1440 
1441  FTensor::Index<'i', 3> i;
1442  FTensor::Index<'j', 3> j;
1443  FTensor::Index<'k', 3> k;
1444  FTensor::Index<'l', 3> l;
1445 
1446  auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hMat));
1447  auto t_H = getFTensor2FromMat<3, 3>(*dataAtPts->HMat);
1448 
1449  dataAtPts->stiffnessMat->resize(36, 1, false);
1450  FTensor::Ddg<FTensor::PackPtr<double *, 1>, 3, 3> t_D(
1451  MAT_TO_DDG(dataAtPts->stiffnessMat));
1452  for (auto &m : (blockSetsPtr)) {
1453  const double young = m.second.E;
1454  const double poisson = m.second.PoissonRatio;
1455 
1456  const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
1457 
1458  t_D(i, j, k, l) = 0.;
1459  t_D(0, 0, 0, 0) = t_D(1, 1, 1, 1) = t_D(2, 2, 2, 2) = 1 - poisson;
1460  t_D(0, 1, 0, 1) = t_D(0, 2, 0, 2) = t_D(1, 2, 1, 2) =
1461  0.5 * (1 - 2 * poisson);
1462  t_D(0, 0, 1, 1) = t_D(1, 1, 0, 0) = t_D(0, 0, 2, 2) = t_D(2, 2, 0, 0) =
1463  t_D(1, 1, 2, 2) = t_D(2, 2, 1, 1) = poisson;
1464  t_D(i, j, k, l) *= coefficient;
1465 
1466  break; // FIXME: calculates only first block
1467  }
1468 
1469  double detH = 0.;
1470  FTensor::Tensor2<double, 3, 3> t_invH;
1471  FTensor::Tensor2<double, 3, 3> t_F;
1472  FTensor::Tensor2<double, 3, 3> t_stress;
1473  FTensor::Tensor2<double, 3, 3> t_small_strain;
1474  FTensor::Tensor2_symmetric<double, 3> t_stress_symm;
1475  FTensor::Tensor2_symmetric<double, 3> t_small_strain_symm;
1476 
1477  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1478 
1479  if (!isALE) {
1480  t_small_strain_symm(i, j) = (t_h(i, j) || t_h(j, i)) / 2.;
1481  } else {
1482  CHKERR determinantTensor3by3(t_H, detH);
1483  CHKERR invertTensor3by3(t_H, detH, t_invH);
1484  t_F(i, j) = t_h(i, k) * t_invH(k, j);
1485  t_small_strain_symm(i, j) = (t_F(i, j) || t_F(j, i)) / 2.;
1486 
1487  t_small_strain_symm(0, 0) -= 1;
1488  t_small_strain_symm(1, 1) -= 1;
1489  t_small_strain_symm(2, 2) -= 1;
1490  ++t_H;
1491  }
1492 
1493  // symmetric tensors need improvement
1494  t_stress_symm(i, j) = t_D(i, j, k, l) * t_small_strain_symm(k, l);
1495  tensor_to_tensor(t_stress_symm, t_stress);
1496 
1497  const double psi = 0.5 * t_stress_symm(i, j) * t_small_strain_symm(i, j);
1498 
1499  CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &psi);
1500  CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
1501  &t_stress(0, 0));
1502 
1503  ++t_h;
1504  }
1505 
1506  MoFEMFunctionReturn(0);
1507 }
1508 
1509 template <int S>
1510 HookeElement::OpAnalyticalInternalStain_dx<S>::OpAnalyticalInternalStain_dx(
1511  const std::string row_field,
1512  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1513  StrainFunctions strain_fun)
1514  : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1515  strainFun(strain_fun) {}
1516 
1517 template <int S>
1518 MoFEMErrorCode
1519 HookeElement::OpAnalyticalInternalStain_dx<S>::iNtegrate(EntData &row_data) {
1520  FTensor::Index<'i', 3> i;
1521  FTensor::Index<'j', 3> j;
1522  FTensor::Index<'k', 3> k;
1523  FTensor::Index<'l', 3> l;
1524  MoFEMFunctionBegin;
1525 
1526  auto get_tensor1 = [](VectorDouble &v, const int r) {
1527  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
1528  &v(r + 0), &v(r + 1), &v(r + 2));
1529  };
1530 
1531  const int nb_integration_pts = getGaussPts().size2();
1532 
1533  auto get_coords = [&]() {
1534  if (getHoCoordsAtGaussPts().size1() == nb_integration_pts)
1535  return getFTensor1HoCoordsAtGaussPts();
1536  else
1537  return getFTensor1CoordsAtGaussPts();
1538  };
1539  auto t_coords = get_coords();
1540 
1541  // get element volume
1542  double vol = getVolume();
1543  auto t_w = getFTensor0IntegrationWeight();
1544 
1545  nF.resize(nbRows, false);
1546  nF.clear();
1547 
1548  // elastic stiffness tensor (4th rank tensor with minor and major
1549  // symmetry)
1550  FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
1551  MAT_TO_DDG(dataAtPts->stiffnessMat));
1552 
1553  // get derivatives of base functions on rows
1554  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1555  const int row_nb_base_fun = row_data.getN().size2();
1556 
1557  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1558 
1559  auto t_fun_strain = strainFun(t_coords);
1560  FTensor::Tensor2_symmetric<double, 3> t_stress;
1561  t_stress(i, j) = t_D(i, j, k, l) * t_fun_strain(k, l);
1562 
1563  // calculate scalar weight times element volume
1564  double a = t_w * vol;
1565 
1566  if (getHoGaussPtsDetJac().size()) {
1567  // If HO geometry
1568  a *= getHoGaussPtsDetJac()[gg];
1569  }
1570 
1571  auto t_nf = get_tensor1(nF, 0);
1572 
1573  int rr = 0;
1574  for (; rr != nbRows / 3; ++rr) {
1575  t_nf(i) += a * t_row_diff_base(j) * t_stress(i, j);
1576  ++t_row_diff_base;
1577  ++t_nf;
1578  }
1579 
1580  for (; rr != row_nb_base_fun; ++rr)
1581  ++t_row_diff_base;
1582 
1583  ++t_w;
1584  ++t_coords;
1585  ++t_D;
1586  }
1587 
1588  MoFEMFunctionReturn(0);
1589 }
1590 
1591 template <int S>
1592 HookeElement::OpAnalyticalInternalAleStain_dX<S>::
1593  OpAnalyticalInternalAleStain_dX(
1594  const std::string row_field,
1595  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1596  StrainFunctions strain_fun,
1597  boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1598  : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1599  strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1600 
1601 template <int S>
1602 MoFEMErrorCode
1603 HookeElement::OpAnalyticalInternalAleStain_dX<S>::iNtegrate(EntData &row_data) {
1604  FTensor::Index<'i', 3> i;
1605  FTensor::Index<'j', 3> j;
1606  FTensor::Index<'k', 3> k;
1607  FTensor::Index<'l', 3> l;
1608  MoFEMFunctionBegin;
1609 
1610  auto get_tensor1 = [](VectorDouble &v, const int r) {
1611  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
1612  &v(r + 0), &v(r + 1), &v(r + 2));
1613  };
1614 
1615  const int nb_integration_pts = getGaussPts().size2();
1616 
1617  auto get_coords = [&]() { return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
1618  auto t_coords = get_coords();
1619 
1620  // get element volume
1621  double vol = getVolume();
1622  auto t_w = getFTensor0IntegrationWeight();
1623 
1624  nF.resize(nbRows, false);
1625  nF.clear();
1626 
1627  // elastic stiffness tensor (4th rank tensor with minor and major
1628  // symmetry)
1629  FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
1630  MAT_TO_DDG(dataAtPts->stiffnessMat));
1631  auto t_F = getFTensor2FromMat<3, 3>(*(dataAtPts->FMat));
1632  auto &det_H = *dataAtPts->detHVec;
1633  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1634 
1635  // get derivatives of base functions on rows
1636  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1637  const int row_nb_base_fun = row_data.getN().size2();
1638 
1639  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1640 
1641  auto t_fun_strain = strainFun(t_coords);
1642  FTensor::Tensor2_symmetric<double, 3> t_stress;
1643  t_stress(i, j) = t_D(i, j, k, l) * t_fun_strain(k, l);
1644  FTensor::Tensor2<double, 3, 3> t_eshelby_stress;
1645  t_eshelby_stress(i, j) = -t_F(k, i) * t_stress(k, j);
1646 
1647  // calculate scalar weight times element volume
1648  double a = t_w * vol * det_H[gg];
1649 
1650  auto t_nf = get_tensor1(nF, 0);
1651 
1652  int rr = 0;
1653  for (; rr != nbRows / 3; ++rr) {
1654  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1655  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1656  t_nf(i) += a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
1657  ++t_row_diff_base;
1658  ++t_nf;
1659  }
1660 
1661  for (; rr != row_nb_base_fun; ++rr)
1662  ++t_row_diff_base;
1663 
1664  ++t_w;
1665  ++t_coords;
1666  ++t_F;
1667  ++t_invH;
1668  ++t_D;
1669  }
1670 
1671  MoFEMFunctionReturn(0);
1672 }
1673 
1674 template <int S>
1675 HookeElement::OpAnalyticalInternalAleStain_dx<S>::
1676  OpAnalyticalInternalAleStain_dx(
1677  const std::string row_field,
1678  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1679  StrainFunctions strain_fun,
1680  boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1681  : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1682  strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1683 
1684 template <int S>
1685 MoFEMErrorCode
1686 HookeElement::OpAnalyticalInternalAleStain_dx<S>::iNtegrate(EntData &row_data) {
1687  FTensor::Index<'i', 3> i;
1688  FTensor::Index<'j', 3> j;
1689  FTensor::Index<'k', 3> k;
1690  FTensor::Index<'l', 3> l;
1691  MoFEMFunctionBegin;
1692 
1693  auto get_tensor1 = [](VectorDouble &v, const int r) {
1694  return FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3>(
1695  &v(r + 0), &v(r + 1), &v(r + 2));
1696  };
1697 
1698  const int nb_integration_pts = getGaussPts().size2();
1699 
1700  auto get_coords = [&]() { return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
1701  auto t_coords = get_coords();
1702 
1703  // get element volume
1704  double vol = getVolume();
1705  auto t_w = getFTensor0IntegrationWeight();
1706 
1707  nF.resize(nbRows, false);
1708  nF.clear();
1709 
1710  // elastic stiffness tensor (4th rank tensor with minor and major
1711  // symmetry)
1712  FTensor::Ddg<FTensor::PackPtr<double *, S>, 3, 3> t_D(
1713  MAT_TO_DDG(dataAtPts->stiffnessMat));
1714  auto &det_H = *dataAtPts->detHVec;
1715  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1716 
1717  // get derivatives of base functions on rows
1718  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1719  const int row_nb_base_fun = row_data.getN().size2();
1720 
1721  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1722 
1723  auto t_fun_strain = strainFun(t_coords);
1724  FTensor::Tensor2_symmetric<double, 3> t_stress;
1725  t_stress(i, j) = t_D(i, j, k, l) * t_fun_strain(k, l);
1726 
1727  // calculate scalar weight times element volume
1728  double a = t_w * vol * det_H[gg];
1729 
1730  auto t_nf = get_tensor1(nF, 0);
1731 
1732  int rr = 0;
1733  for (; rr != nbRows / 3; ++rr) {
1734  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1735  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1736  t_nf(i) += a * t_row_diff_base_pulled(j) * t_stress(i, j);
1737  ++t_row_diff_base;
1738  ++t_nf;
1739  }
1740 
1741  for (; rr != row_nb_base_fun; ++rr)
1742  ++t_row_diff_base;
1743 
1744  ++t_w;
1745  ++t_coords;
1746  ++t_invH;
1747  ++t_D;
1748  }
1749 
1750  MoFEMFunctionReturn(0);
1751 }
1752 
1753 #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:482
OpLhs_dx_dx
Definition: HookeElement.hpp:426
OpCalculateStiffnessScaledByDensityField::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:339
DataAtIntegrationPts::eshelbyStress_dx
boost::shared_ptr< MatrixDouble > eshelbyStress_dx
Definition: HookeElement.hpp:107
EntData
DataForcesAndSourcesCore::EntData EntData
Definition: HookeElement.hpp:87
OpAssemble::nbCols
int nbCols
number if dof on column
Definition: HookeElement.hpp:390
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:75
OpPlasticTools::n
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
DataAtIntegrationPts::hMat
boost::shared_ptr< MatrixDouble > hMat
Definition: HookeElement.hpp:95
invJac
MatrixDouble invJac
Definition: HookeElement.hpp:697
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
ConvectiveMassElement
structure grouping operators and data used for calculation of mass (convective) element \ nonlinear_e...
Definition: ConvectiveMassElement.hpp:42
OpAssemble::K
MatrixDouble K
Definition: HookeElement.hpp:380
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:514
OpAleRhs_dx
Definition: HookeElement.hpp:441
OpAssemble::isDiag
bool isDiag
true if this block is on diagonal
Definition: HookeElement.hpp:392
OpCalculateStress::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:175
DataAtIntegrationPts::invHMat
boost::shared_ptr< MatrixDouble > invHMat
Definition: HookeElement.hpp:100
OpAnalyticalInternalAleStain_dx
Definition: HookeElement.hpp:627
OpAleLhsWithDensity_dx_dX::rHo0
const double rHo0
Definition: HookeElement.hpp:485
MoFEM::invertTensor3by3
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:463
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:1150
NonlinearElasticElement::BlockData::tEts
Range tEts
constrains elements in block set
Definition: HookeElement.hpp:49
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
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
OpPlasticTools::m
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
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
DataAtIntegrationPts::eshelbyStressMat
boost::shared_ptr< MatrixDouble > eshelbyStressMat
Definition: HookeElement.hpp:105
DataAtIntegrationPts::detHVec
boost::shared_ptr< VectorDouble > detHVec
Definition: HookeElement.hpp:99
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: DataStructures.hpp:1552
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:602
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
OpContactTools::k
FTensor::Index< 'k', 2 > k
Definition: ContactOps.hpp:28
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:72
OpAssemble::nF
VectorDouble nF
Definition: HookeElement.hpp:382
calculateEnergy
static MoFEMErrorCode calculateEnergy(DM dm, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, Vec *v_energy_ptr)
DataAtIntegrationPts::energyVec
boost::shared_ptr< VectorDouble > energyVec
Definition: HookeElement.hpp:104
OpAleLhs_dx_dx
Definition: HookeElement.hpp:450
OpAnalyticalInternalAleStain_dx::strainFun
StrainFunctions strainFun
Definition: HookeElement.hpp:648
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
BlockData
struct ConvectiveMassElement BlockData
MoFEM::DataForcesAndSourcesCore::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: DataStructures.hpp:53
OpContactTools::j
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
OpCalculateEnergy::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:189
OpAleLhsPre_dX_dx
Definition: HookeElement.hpp:552
OpAssemble::rowIndices
VectorInt rowIndices
Definition: HookeElement.hpp:386
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
OpCalculateMassMatrix::doWork
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Definition: HookeElement.hpp:247
OpPostProcHookeElement::postProcMesh
moab::Interface & postProcMesh
Definition: HookeElement.hpp:657
MoFEM::Types::VectorDouble
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:73
MoFEM::DataForcesAndSourcesCore::EntData::getN
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
Definition: DataStructures.hpp:1251
OpAleRhs_dX
Definition: HookeElement.hpp:528
NonlinearElasticElement::BlockData
data for calculation heat conductivity and heat capacity elements
Definition: HookeElement.hpp:45
MoFEM::ForcesAndSourcesCore::UserDataOperator
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
Definition: ForcesAndSourcesCore.hpp:84
ConvectiveMassElement::BlockData
data for calculation inertia forces
Definition: ConvectiveMassElement.hpp:116
OpCalculateMassMatrix::massData
MassBlockData & massData
Definition: HookeElement.hpp:234
OpContactTools::l
FTensor::Index< 'l', 2 > l
Definition: ContactOps.hpp:29
OpAleLhsPre_dX_dx::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: HookeElement.hpp:560
DataAtIntegrationPts::cauchyStressMat
boost::shared_ptr< MatrixDouble > cauchyStressMat
Definition: HookeElement.hpp:102
OpCalculateMassMatrix
Assemble mass matrix for elastic element TODO: CHANGE FORMULA * .
Definition: HookeElement.hpp:228
MoFEM::DataForcesAndSourcesCore::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: DataStructures.hpp:1422
OpContactTools::w
double w(const double g, const double t)
Definition: ContactOps.hpp:160
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic elementIn order to as...
Definition: HookeElement.hpp:40
OpContactTools::i
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26
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 Sun Aug 2 2020 20:49:22 for MoFEM by Doxygen 1.8.15 and hosted at University of Glasgow