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