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 
41 
42  /** \brief data for calculation heat conductivity and heat capacity elements
43  * \ingroup nonlinear_elastic_elem
44  */
45  struct BlockData {
46  int iD;
47  double E;
48  double PoissonRatio;
49  Range tEts; ///< constrains elements in block set
52  };
53 };
54 
55 #endif // __NONLINEAR_ELASTIC_HPP
56 
57 /** \brief structure grouping operators and data used for calculation of
58  * nonlinear elastic element \ingroup nonlinear_elastic_elem
59  *
60  * In order to assemble matrices and right hand vectors, the loops over
61  * elements, entities over that elements and finally loop over integration
62  * points are executed.
63  *
64  * Following implementation separate those three categories of loops and to each
65  * loop attach operator.
66  *
67  */
68 
69 #ifndef __CONVECTIVE_MASS_ELEMENT_HPP
70 struct ConvectiveMassElement {
71  /** \brief data for calculation inertia forces
72  * \ingroup user_modules
73  */
74  struct BlockData {
75  double rho0; ///< reference density
76  VectorDouble a0; ///< constant acceleration
77  Range tEts; ///< elements in block set
78  };
79 }
80 
81 #endif //__CONVECTIVE_MASS_ELEMENT_HPP
82 struct HookeElement {
83 
86 
89  using VolUserDataOperator =
91 
92  struct DataAtIntegrationPts {
93 
94  boost::shared_ptr<MatrixDouble> smallStrainMat;
95  boost::shared_ptr<MatrixDouble> hMat;
96  boost::shared_ptr<MatrixDouble> FMat;
97 
98  boost::shared_ptr<MatrixDouble> HMat;
99  boost::shared_ptr<VectorDouble> detHVec;
100  boost::shared_ptr<MatrixDouble> invHMat;
101 
102  boost::shared_ptr<MatrixDouble> cauchyStressMat;
103  boost::shared_ptr<MatrixDouble> stiffnessMat;
104  boost::shared_ptr<VectorDouble> energyVec;
105  boost::shared_ptr<MatrixDouble> eshelbyStressMat;
106 
107  boost::shared_ptr<MatrixDouble> eshelbyStress_dx;
108 
110 
111  smallStrainMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
112  hMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
113  FMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
114 
115  HMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
116  detHVec = boost::shared_ptr<VectorDouble>(new VectorDouble());
117  invHMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
118 
119  cauchyStressMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
120  stiffnessMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
121  energyVec = boost::shared_ptr<VectorDouble>(new VectorDouble());
122  eshelbyStressMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
123  stiffnessMat = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
124 
125  eshelbyStress_dx = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
126  }
127 
130  };
131 
132  template <bool D = false>
134 
135  OpCalculateStrain(const std::string row_field, const std::string col_field,
136  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
137 
138  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
139 
140  private:
141  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
142  };
143 
145 
146  OpCalculateStrainAle(const std::string row_field,
147  const std::string col_field,
148  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
149 
150  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
151 
152  private:
153  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
154  };
155 
156 #define MAT_TO_DDG(SM) \
157  &(*SM)(0, 0), &(*SM)(1, 0), &(*SM)(2, 0), &(*SM)(3, 0), &(*SM)(4, 0), \
158  &(*SM)(5, 0), &(*SM)(6, 0), &(*SM)(7, 0), &(*SM)(8, 0), &(*SM)(9, 0), \
159  &(*SM)(10, 0), &(*SM)(11, 0), &(*SM)(12, 0), &(*SM)(13, 0), \
160  &(*SM)(14, 0), &(*SM)(15, 0), &(*SM)(16, 0), &(*SM)(17, 0), \
161  &(*SM)(18, 0), &(*SM)(19, 0), &(*SM)(20, 0), &(*SM)(21, 0), \
162  &(*SM)(22, 0), &(*SM)(23, 0), &(*SM)(24, 0), &(*SM)(25, 0), \
163  &(*SM)(26, 0), &(*SM)(27, 0), &(*SM)(28, 0), &(*SM)(29, 0), \
164  &(*SM)(30, 0), &(*SM)(31, 0), &(*SM)(32, 0), &(*SM)(33, 0), \
165  &(*SM)(34, 0), &(*SM)(35, 0)
166 
167  template <int S = 0> struct OpCalculateStress : public VolUserDataOperator {
168 
169  OpCalculateStress(const std::string row_field, const std::string col_field,
170  boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
171 
172  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
173 
174  protected:
175  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
176  };
177 
179 
180  OpCalculateEnergy(const std::string row_field, const std::string col_field,
181  boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
182  Vec ghost_vec = PETSC_NULL);
183 
185 
186  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
187 
188  protected:
189  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
190  Vec ghostVec;
191  };
192 
194 
196  const std::string row_field, const std::string col_field,
197  boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
198 
199  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
200 
201  protected:
202  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
203  };
204 
205  template <int S = 0>
207 
209  const std::string row_field, const std::string col_field,
210  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
211  boost::shared_ptr<DataAtIntegrationPts> data_at_pts);
212 
213  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
214 
215  protected:
216  boost::shared_ptr<map<int, BlockData>>
217  blockSetsPtr; ///< Structure keeping data about problem, like
218  ///< material parameters
219  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
220  };
221 
222  /** * @brief Assemble mass matrix for elastic element TODO: CHANGE FORMULA *
223  * \f[
224  * {\bf{M}} = \int\limits_\Omega
225  * \f]
226  *
227  */
230 
235 
236  boost::shared_ptr<DataAtIntegrationPts> commonData;
237 
238  OpCalculateMassMatrix(const std::string row_field,
239  const std::string col_field, BlockData &data,
240  MassBlockData &mass_data,
241  boost::shared_ptr<DataAtIntegrationPts> &common_data,
242  bool symm = true)
244  row_field, col_field, OPROWCOL, symm),
245  commonData(common_data), dAta(data), massData(mass_data) {}
246 
247  PetscErrorCode doWork(int row_side, int col_side, EntityType row_type,
248  EntityType col_type,
252 
253  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
255  &m(3 * r + 0, 3 * c + 0), &m(3 * r + 0, 3 * c + 1),
256  &m(3 * r + 0, 3 * c + 2), &m(3 * r + 1, 3 * c + 0),
257  &m(3 * r + 1, 3 * c + 1), &m(3 * r + 1, 3 * c + 2),
258  &m(3 * r + 2, 3 * c + 0), &m(3 * r + 2, 3 * c + 1),
259  &m(3 * r + 2, 3 * c + 2));
260  };
261 
262  const int row_nb_dofs = row_data.getIndices().size();
263  if (!row_nb_dofs)
265  const int col_nb_dofs = col_data.getIndices().size();
266  if (!col_nb_dofs)
268  if (dAta.tEts.find(getFEEntityHandle()) == dAta.tEts.end()) {
270  }
271  if (massData.tEts.find(getFEEntityHandle()) == massData.tEts.end()) {
273  }
274 
275  const bool diagonal_block =
276  (row_type == col_type) && (row_side == col_side);
277  // get number of integration points
278  // Set size can clear local tangent matrix
279  locK.resize(row_nb_dofs, col_nb_dofs, false);
280  locK.clear();
281 
282  const int row_nb_gauss_pts = row_data.getN().size1();
283  const int row_nb_base_functions = row_data.getN().size2();
284 
289 
290  double density = massData.rho0;
291 
292  // get integration weights
293  auto t_w = getFTensor0IntegrationWeight();
294 
295  // integrate local matrix for entity block
296  for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
297 
298  auto t_row_base_func = row_data.getFTensor0N(gg, 0);
299 
300  // Get volume and integration weight
301  double w = getVolume() * t_w;
302 
303  for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
304  auto t_col_base_func = col_data.getFTensor0N(gg, 0);
305  for (int col_bb = 0; col_bb != col_nb_dofs / 3; col_bb++) {
306  auto t_assemble = get_tensor2(locK, row_bb, col_bb);
307  t_assemble(i, j) += density * t_row_base_func * t_col_base_func * w;
308  // Next base function for column
309  ++t_col_base_func;
310  }
311  // Next base function for row
312  ++t_row_base_func;
313  }
314  // Next integration point for getting weight
315  ++t_w;
316  }
317 
318  CHKERR MatSetValues(getKSPB(), row_data, col_data, &locK(0, 0),
319  ADD_VALUES);
320 
321  // is symmetric
322  if (row_type != col_type || row_side != col_side) {
323  translocK.resize(col_nb_dofs, row_nb_dofs, false);
324  noalias(translocK) = trans(locK);
325 
326  CHKERR MatSetValues(getKSPB(), col_data, row_data, &translocK(0, 0),
327  ADD_VALUES);
328  }
329 
331  }
332  };
333 
335  protected:
336  boost::shared_ptr<map<int, BlockData>>
337  blockSetsPtr; ///< Structure keeping data about problem, like
338  ///< material parameters
339  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
340 
341  boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
342  const double rhoN; ///< exponent n in E(p) = E * (p / p_0)^n
343  const double rHo0; ///< p_0 reference density in E(p) = E * (p / p_0)^n
344  // // where p is density, E - youngs modulus
345  public:
347  const std::string row_field, const std::string col_field,
348  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
349  boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
350  boost::shared_ptr<VectorDouble> rho_at_gauss_pts, const double rho_n,
351  const double rho_0);
352 
353  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
354  };
355 
356  struct OpAssemble : public VolUserDataOperator {
357 
358  OpAssemble(const std::string row_field, const std::string col_field,
359  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
360  const char type, bool symm = false);
361 
362  /**
363  * \brief Do calculations for give operator
364  * @param row_side row side number (local number) of entity on element
365  * @param col_side column side number (local number) of entity on element
366  * @param row_type type of row entity MBVERTEX, MBEDGE, MBTRI or MBTET
367  * @param col_type type of column entity MBVERTEX, MBEDGE, MBTRI or MBTET
368  * @param row_data data for row
369  * @param col_data data for column
370  * @return error code
371  */
372  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
373  EntityType col_type, EntData &row_data,
374  EntData &col_data);
375 
376  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
377 
378  protected:
379  // Finite element stiffness sub-matrix K_ij
383 
384  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
385 
388 
389  int nbRows; ///< number of dofs on rows
390  int nbCols; ///< number if dof on column
391  int nbIntegrationPts; ///< number of integration points
392  bool isDiag; ///< true if this block is on diagonal
393 
394  virtual MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
395 
396  virtual MoFEMErrorCode iNtegrate(EntData &row_data);
397 
398  /**
399  * \brief Assemble local entity block matrix
400  * @param row_data row data (consist base functions on row entity)
401  * @param col_data column data (consist base functions on column
402  * entity)
403  * @return error code
404  */
405  MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data);
406 
407  /**
408  * \brief Assemble local entity right-hand vector
409  * @param row_data row data (consist base functions on row entity)
410  * @param col_data column data (consist base functions on column
411  * entity)
412  * @return error code
413  */
414  MoFEMErrorCode aSsemble(EntData &row_data);
415  };
416 
417  struct OpRhs_dx : public OpAssemble {
418 
419  OpRhs_dx(const std::string row_field, const std::string col_field,
420  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
421 
422  protected:
423  MoFEMErrorCode iNtegrate(EntData &row_data);
424  };
425 
426  template <int S = 0> struct OpLhs_dx_dx : public OpAssemble {
427 
428  OpLhs_dx_dx(const std::string row_field, const std::string col_field,
429  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
430 
431  protected:
432  /**
433  * \brief Integrate B^T D B operator
434  * @param row_data row data (consist base functions on row entity)
435  * @param col_data column data (consist base functions on column entity)
436  * @return error code
437  */
438  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
439  };
440 
441  struct OpAleRhs_dx : public OpAssemble {
442 
443  OpAleRhs_dx(const std::string row_field, const std::string col_field,
444  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
445 
446  protected:
447  MoFEMErrorCode iNtegrate(EntData &row_data);
448  };
449 
450  template <int S = 0> struct OpAleLhs_dx_dx : public OpAssemble {
451 
452  OpAleLhs_dx_dx(const std::string row_field, const std::string col_field,
453  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
454 
455  protected:
456  /**
457  * \brief Integrate B^T D B operator
458  * @param row_data row data (consist base functions on row entity)
459  * @param col_data column data (consist base functions on column entity)
460  * @return error code
461  */
462  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
463  };
464 
465  template <int S = 0> struct OpAleLhs_dx_dX : public OpAssemble {
466 
467  OpAleLhs_dx_dX(const std::string row_field, const std::string col_field,
468  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
469 
470  protected:
471  /**
472  * \brief Integrate tangent stiffness for spatial momentum
473  * @param row_data row data (consist base functions on row entity)
474  * @param col_data column data (consist base functions on column entity)
475  * @return error code
476  */
477  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
478  };
479 
481 
482  boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
483  boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
484  const double rhoN;
485  const double rHo0;
486 
488  const std::string row_field, const std::string col_field,
489  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
490  boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
491  boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
492  const double rho_n, const double rho_0);
493 
494  protected:
495  /**
496  * \brief Integrate tangent stiffness for spatial momentum
497  * @param row_data row data (consist base functions on row entity)
498  * @param col_data column data (consist base functions on column entity)
499  * @return error code
500  */
501  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
502  };
503 
505 
506  boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
507  boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
508  const double rhoN;
509  const double rHo0;
510 
512  const std::string row_field, const std::string col_field,
513  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
514  boost::shared_ptr<VectorDouble> rho_at_gauss_pts,
515  boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts,
516  const double rho_n, const double rho_0);
517 
518  protected:
519  /**
520  * \brief Integrate tangent stiffness for material momentum
521  * @param row_data row data (consist base functions on row entity)
522  * @param col_data column data (consist base functions on column entity)
523  * @return error code
524  */
525  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
526  };
527 
528  struct OpAleRhs_dX : public OpAssemble {
529 
530  OpAleRhs_dX(const std::string row_field, const std::string col_field,
531  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
532 
533  protected:
534  MoFEMErrorCode iNtegrate(EntData &row_data);
535  };
536 
537  template <int S = 0> struct OpAleLhs_dX_dX : public OpAssemble {
538 
539  OpAleLhs_dX_dX(const std::string row_field, const std::string col_field,
540  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
541 
542  protected:
543  /**
544  * \brief Integrate tangent stiffness for material momentum
545  * @param row_data row data (consist base functions on row entity)
546  * @param col_data column data (consist base functions on column entity)
547  * @return error code
548  */
549  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
550  };
551 
552  template <int S = 0> struct OpAleLhsPre_dX_dx : public VolUserDataOperator {
553 
554  OpAleLhsPre_dX_dx(const std::string row_field, const std::string col_field,
555  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts);
556 
557  MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data);
558 
559  private:
560  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
561  };
562 
563  struct OpAleLhs_dX_dx : public OpAssemble {
564 
565  OpAleLhs_dX_dx(const std::string row_field, const std::string col_field,
566  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
567  : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
568 
569  protected:
570  /**
571  * \brief Integrate tangent stiffness for material momentum
572  * @param row_data row data (consist base functions on row entity)
573  * @param col_data column data (consist base functions on column entity)
574  * @return error code
575  */
576  MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data);
577  };
578 
579  template <int S> struct OpAnalyticalInternalStain_dx : public OpAssemble {
580 
581  typedef boost::function<
582 
584 
586 
587  )
588 
589  >
591 
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);
600  };
601 
602  template <int S> struct OpAnalyticalInternalAleStain_dX : public OpAssemble {
603 
604  typedef boost::function<
605 
607 
609 
610  )
611 
612  >
614 
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);
624  boost::shared_ptr<MatrixDouble> matPosAtPtsPtr;
625  };
626 
627  template <int S> struct OpAnalyticalInternalAleStain_dx : public OpAssemble {
628 
629  typedef boost::function<
630 
632 
634 
635  )
636 
637  >
639 
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);
649  boost::shared_ptr<MatrixDouble> matPosAtPtsPtr;
650  };
651 
652  template <class ELEMENT>
654  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
655  map<int, BlockData>
656  &blockSetsPtr; // FIXME: (works only with the first block)
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,
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
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:
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) {
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  }
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>
746  EntData &col_data) {
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) {
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 
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)
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 
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 
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) {
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);
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 
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 
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) {
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 
922 
923  // elastic stiffness tensor (4th rank tensor with minor and major
924  // symmetry)
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  }
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>
944  EntData &col_data) {
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) {
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 
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)
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 
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 
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>
1037  EntData &col_data) {
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) {
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 
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)
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 
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 
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 
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>
1166  EntData &col_data) {
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) {
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 
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)
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 
1216  t_F_dX(i, j, k, l) = -(t_h(i, m) * t_invH(m, k)) * t_invH(l, j);
1217 
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 
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) {
1334 
1335  const int nb_integration_pts = row_data.getN().size1();
1336 
1337  auto get_eshelby_stress_dx = [this, nb_integration_pts]() {
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 
1359 
1360  // Elastic stiffness tensor (4th rank tensor with minor and major
1361  // symmetry)
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 
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) {
1415 
1416  if (type != MBVERTEX) {
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 
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);
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.;
1473  FTensor::Tensor2<double, 3, 3> t_small_strain;
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 
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>
1525 
1526  auto get_tensor1 = [](VectorDouble &v, const int r) {
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)
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);
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 
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>
1609 
1610  auto get_tensor1 = [](VectorDouble &v, const int r) {
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)
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);
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 
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>
1692 
1693  auto get_tensor1 = [](VectorDouble &v, const int r) {
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)
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);
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 
1751 }
1752 
1753 #endif // __HOOKE_ELEMENT_HPP
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > matPosAtPtsPtr
MoFEMErrorCode aSsemble(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Assemble local entity block matrix.
bool isALE
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Integrate B^T D B operator.
Deprecated interface functions.
boost::shared_ptr< MatrixDouble > smallStrainMat
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
static MoFEMErrorCode setOperators(boost::shared_ptr< ForcesAndSourcesCore > fe_lhs_ptr, boost::shared_ptr< ForcesAndSourcesCore > fe_rhs_ptr, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, const EntityType type=MBTET, boost::shared_ptr< DataAtIntegrationPts > data_at_pts=nullptr)
MoFEMErrorCode determinantTensor3by3(T1 &t, T2 &det)
Calculate determinant 3 by 3.
Definition: Templates.hpp:482
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< MatrixDouble > eshelbyStress_dx
DataForcesAndSourcesCore::EntData EntData
int nbCols
number if dof on column
boost::shared_ptr< map< int, BlockData > > blockSetsPtr
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
boost::shared_ptr< MatrixDouble > hMat
MatrixDouble invJac
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
structure grouping operators and data used for calculation of mass (convective) element \ nonlinear_e...
MatrixDouble K
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
FTensor::Tensor1< double *, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > stiffnessMat
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
bool isDiag
true if this block is on diagonal
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< MatrixDouble > invHMat
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:463
const double rhoN
exponent n in E(p) = E * (p / p_0)^n
std::vector< EntityHandle > & mapGaussPts
boost::shared_ptr< MatrixDouble > HMat
OpAleLhs_dX_dx(const std::string row_field, const std::string col_field, boost::shared_ptr< DataAtIntegrationPts > &data_at_pts)
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Range tEts
constrains elements in block set
boost::function< FTensor::Tensor2_symmetric< double, 3 > FTensor::Tensor1< FTensor::PackPtr< double *, 3 >, 3 > &t_coords) > StrainFunctions
MatrixDouble transK
boost::shared_ptr< MatrixDouble > FMat
map< int, BlockData > & blockSetsPtr
ForcesAndSourcesCore::UserDataOperator UserDataOperator
boost::shared_ptr< MatrixDouble > matPosAtPtsPtr
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
OpCalculateMassMatrix(const std::string row_field, const std::string col_field, BlockData &data, MassBlockData &mass_data, boost::shared_ptr< DataAtIntegrationPts > &common_data, bool symm=true)
int nbIntegrationPts
number of integration points
boost::shared_ptr< map< int, BlockData > > blockSetsPtr
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
VectorDouble a0
constant acceleration
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< MatrixDouble > eshelbyStressMat
boost::shared_ptr< VectorDouble > detHVec
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
int nbRows
number of dofs on rows
boost::shared_ptr< DataAtIntegrationPts > commonData
#define MAT_TO_DDG(SM)
VectorInt colIndices
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
static MoFEMErrorCode setBlocks(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr)
static MoFEMErrorCode addElasticElement(MoFEM::Interface &m_field, boost::shared_ptr< map< int, BlockData >> &block_sets_ptr, const std::string element_name, const std::string x_field, const std::string X_field, const bool ale)
const double rHo0
p_0 reference density in E(p) = E * (p / p_0)^n
#define CHKERR
Inline error check.
Definition: definitions.h:602
boost::function< FTensor::Tensor2_symmetric< double, 3 > FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords) > StrainFunctions
FTensor::Index< 'k', 2 > k
Definition: ContactOps.hpp:28
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:72
VectorDouble nF
static MoFEMErrorCode calculateEnergy(DM dm, boost::shared_ptr< map< int, BlockData >> block_sets_ptr, const std::string x_field, const std::string X_field, const bool ale, const bool field_disp, Vec *v_energy_ptr)
boost::shared_ptr< VectorDouble > energyVec
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
struct ConvectiveMassElement BlockData
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
VectorInt rowIndices
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
PetscErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
moab::Interface & postProcMesh
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:73
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
data for calculation heat conductivity and heat capacity elements
Data operator to do calculations at integration points.Is inherited and implemented by user to do cal...
data for calculation inertia forces
MassBlockData & massData
FTensor::Index< 'l', 2 > l
Definition: ContactOps.hpp:29
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< MatrixDouble > cauchyStressMat
Assemble mass matrix for elastic element TODO: CHANGE FORMULA * .
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
double w(const double g, const double t)
Definition: ContactOps.hpp:160
structure grouping operators and data used for calculation of nonlinear elastic elementIn order to as...
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26
boost::function< FTensor::Tensor2_symmetric< double, 3 > FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords) > StrainFunctions