v0.10.0
HookeElement.hpp
Go to the documentation of this file.
1 /**
2  * \file HookeElement.hpp
3  * \example HookeElement.hpp
4  *
5  * \brief Operators and data structures for linear elastic
6  * analysis
7  *
8  * Implemention of operators for Hooke material. Implementation is extended to
9  * the case when the mesh is moving as results of topological changes, also the
10  * calculation of material forces and associated tangent matrices are added to
11  * implementation.
12  *
13  * In other words, spatial deformation is small but topological changes large.
14  */
15 
16 /*
17  * This file is part of MoFEM.
18  * MoFEM is free software: you can redistribute it and/or modify it under
19  * the terms of the GNU Lesser General Public License as published by the
20  * Free Software Foundation, either version 3 of the License, or (at your
21  * option) any later version.
22  *
23  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
24  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
25  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
26  * License for more details.
27  *
28  * You should have received a copy of the GNU Lesser General Public
29  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
30 
31 #ifndef __HOOKE_ELEMENT_HPP
32 #define __HOOKE_ELEMENT_HPP
33 
34 #ifndef __BASICFINITEELEMENTS_HPP__
35 #include <BasicFiniteElements.hpp>
36 #endif // __BASICFINITEELEMENTS_HPP__
37 
38 #ifndef __NONLINEAR_ELASTIC_HPP
39 
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 
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;
661 
662  OpPostProcHookeElement(const string row_field,
663  boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
664  map<int, BlockData> &block_sets_ptr,
665  moab::Interface &post_proc_mesh,
666  std::vector<EntityHandle> &map_gauss_pts,
667  bool is_ale = false, bool is_field_disp = true);
668 
669  MoFEMErrorCode doWork(int side, EntityType type,
671  };
672 
673  static MoFEMErrorCode
674  setBlocks(MoFEM::Interface &m_field,
675  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr);
676 
677  static MoFEMErrorCode
679  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
680  const std::string element_name, const std::string x_field,
681  const std::string X_field, const bool ale);
682 
683  static MoFEMErrorCode
684  setOperators(boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr,
685  boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr,
686  boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
687  const std::string x_field, const std::string X_field,
688  const bool ale, const bool field_disp,
689  const EntityType type = MBTET,
690  boost::shared_ptr<DataAtIntegrationPts> data_at_pts = nullptr);
691 
692  static MoFEMErrorCode
693  calculateEnergy(DM dm, boost::shared_ptr<map<int, BlockData>> block_sets_ptr,
694  const std::string x_field, const std::string X_field,
695  const bool ale, const bool field_disp, Vec *v_energy_ptr);
696 
697 private:
699 };
700 
701 template <bool D>
702 HookeElement::OpCalculateStrain<D>::OpCalculateStrain(
703  const std::string row_field, const std::string col_field,
704  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
705  : VolUserDataOperator(row_field, col_field, OPROW, true),
706  dataAtPts(data_at_pts) {
707  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
708 }
709 
710 template <bool D>
711 MoFEMErrorCode HookeElement::OpCalculateStrain<D>::doWork(int row_side,
712  EntityType row_type,
713  EntData &row_data) {
717  // get number of integration points
718  const int nb_integration_pts = getGaussPts().size2();
719  dataAtPts->smallStrainMat->resize(6, nb_integration_pts, false);
720  auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
721  auto t_h = getFTensor2FromMat<3, 3>(*(dataAtPts->hMat));
722 
723  for (int gg = 0; gg != nb_integration_pts; ++gg) {
724  t_strain(i, j) = (t_h(i, j) || t_h(j, i)) / 2.;
725 
726  // If displacement field, not field o spatial positons is given
727  if (!D) {
728  t_strain(0, 0) -= 1;
729  t_strain(1, 1) -= 1;
730  t_strain(2, 2) -= 1;
731  }
732 
733  ++t_strain;
734  ++t_h;
735  }
737 }
738 
739 template <int S>
740 HookeElement::OpAleLhs_dx_dx<S>::OpAleLhs_dx_dx(
741  const std::string row_field, const std::string col_field,
742  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
743  : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
744 
745 template <int S>
747  EntData &col_data) {
749 
750  // get sub-block (3x3) of local stiffens matrix, here represented by
751  // second order tensor
752  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
754  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
755  &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
756  &m(r + 2, c + 2));
757  };
758 
763 
764  // get element volume
765  double vol = getVolume();
766 
767  // get intergrayion weights
768  auto t_w = getFTensor0IntegrationWeight();
769 
770  // get derivatives of base functions on rows
771  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
772  const int row_nb_base_fun = row_data.getN().size2();
773 
774  // Elastic stiffness tensor (4th rank tensor with minor and major
775  // symmetry)
777  MAT_TO_DDG(dataAtPts->stiffnessMat));
778 
779  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
780  auto &det_H = *dataAtPts->detHVec;
781 
782  // iterate over integration points
783  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
784 
785  // calculate scalar weight times element volume
786  double a = t_w * vol * det_H[gg];
787 
788  // iterate over row base functions
789  int rr = 0;
790  for (; rr != nbRows / 3; ++rr) {
791 
792  // get sub matrix for the row
793  auto t_m = get_tensor2(K, 3 * rr, 0);
794 
795  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
796  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
797 
799  // I mix up the indices here so that it behaves like a
800  // Dg. That way I don't have to have a separate wrapper
801  // class Christof_Expr, which simplifies things.
802  t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base_pulled(i));
803 
804  // get derivatives of base functions for columns
805  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
806 
807  // iterate column base functions
808  for (int cc = 0; cc != nbCols / 3; ++cc) {
809 
810  FTensor::Tensor1<double, 3> t_col_diff_base_pulled;
811  t_col_diff_base_pulled(j) = t_col_diff_base(i) * t_invH(i, j);
812 
813  // integrate block local stiffens matrix
814  t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base_pulled(k);
815 
816  // move to next column base function
817  ++t_col_diff_base;
818 
819  // move to next block of local stiffens matrix
820  ++t_m;
821  }
822 
823  // move to next row base function
824  ++t_row_diff_base;
825  }
826 
827  for (; rr != row_nb_base_fun; ++rr)
828  ++t_row_diff_base;
829 
830  // move to next integration weight
831  ++t_w;
832  ++t_D;
833  ++t_invH;
834  }
835 
837 }
838 
839 template <int S>
840 HookeElement::OpCalculateHomogeneousStiffness<S>::
841  OpCalculateHomogeneousStiffness(
842  const std::string row_field, const std::string col_field,
843  boost::shared_ptr<map<int, BlockData>> &block_sets_ptr,
844  boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
845  : VolUserDataOperator(row_field, col_field, OPROW, true),
846  blockSetsPtr(block_sets_ptr), dataAtPts(data_at_pts) {
847  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
848 }
849 
850 template <int S>
851 MoFEMErrorCode HookeElement::OpCalculateHomogeneousStiffness<S>::doWork(
852  int row_side, EntityType row_type, EntData &row_data) {
854 
855  for (auto &m : (*blockSetsPtr)) {
856  if (m.second.tEts.find(getFEEntityHandle()) != m.second.tEts.end()) {
857 
858  dataAtPts->stiffnessMat->resize(36, 1, false);
860  MAT_TO_DDG(dataAtPts->stiffnessMat));
861  const double young = m.second.E;
862  const double poisson = m.second.PoissonRatio;
863 
864  // coefficient used in intermediate calculation
865  const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
866 
871 
872  t_D(i, j, k, l) = 0.;
873 
874  t_D(0, 0, 0, 0) = 1 - poisson;
875  t_D(1, 1, 1, 1) = 1 - poisson;
876  t_D(2, 2, 2, 2) = 1 - poisson;
877 
878  t_D(0, 1, 0, 1) = 0.5 * (1 - 2 * poisson);
879  t_D(0, 2, 0, 2) = 0.5 * (1 - 2 * poisson);
880  t_D(1, 2, 1, 2) = 0.5 * (1 - 2 * poisson);
881 
882  t_D(0, 0, 1, 1) = poisson;
883  t_D(1, 1, 0, 0) = poisson;
884  t_D(0, 0, 2, 2) = poisson;
885  t_D(2, 2, 0, 0) = poisson;
886  t_D(1, 1, 2, 2) = poisson;
887  t_D(2, 2, 1, 1) = poisson;
888 
889  t_D(i, j, k, l) *= coefficient;
890 
891  break;
892  }
893  }
894 
896 }
897 
898 template <int S>
899 HookeElement::OpCalculateStress<S>::OpCalculateStress(
900  const std::string row_field, const std::string col_field,
901  boost::shared_ptr<DataAtIntegrationPts> data_at_pts)
902  : VolUserDataOperator(row_field, col_field, OPROW, true),
903  dataAtPts(data_at_pts) {
904  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
905 }
906 
907 template <int S>
908 MoFEMErrorCode HookeElement::OpCalculateStress<S>::doWork(int row_side,
909  EntityType row_type,
910  EntData &row_data) {
912  // get number of integration points
913  const int nb_integration_pts = getGaussPts().size2();
914  auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
915  dataAtPts->cauchyStressMat->resize(6, nb_integration_pts, false);
916  auto t_cauchy_stress =
917  getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
918 
923 
924  // elastic stiffness tensor (4th rank tensor with minor and major
925  // symmetry)
927  MAT_TO_DDG(dataAtPts->stiffnessMat));
928  for (int gg = 0; gg != nb_integration_pts; ++gg) {
929  t_cauchy_stress(i, j) = t_D(i, j, k, l) * t_strain(k, l);
930  ++t_strain;
931  ++t_cauchy_stress;
932  ++t_D;
933  }
935 }
936 
937 template <int S>
938 HookeElement::OpLhs_dx_dx<S>::OpLhs_dx_dx(
939  const std::string row_field, const std::string col_field,
940  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
941  : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
942 
943 template <int S>
945  EntData &col_data) {
947 
948  // get sub-block (3x3) of local stiffens matrix, here represented by
949  // second order tensor
950  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
952  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
953  &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
954  &m(r + 2, c + 2));
955  };
956 
961 
962  // get element volume
963  double vol = getVolume();
964 
965  // get intergrayion weights
966  auto t_w = getFTensor0IntegrationWeight();
967 
968  // get derivatives of base functions on rows
969  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
970  const int row_nb_base_fun = row_data.getN().size2();
971 
972  // Elastic stiffness tensor (4th rank tensor with minor and major
973  // symmetry)
975  MAT_TO_DDG(dataAtPts->stiffnessMat));
976 
977  // iterate over integration points
978  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
979 
980  // calculate scalar weight times element volume
981  double a = t_w * vol;
982  if (getHoGaussPtsDetJac().size()) {
983  a *= getHoGaussPtsDetJac()[gg];
984  }
985 
986  // iterate over row base functions
987  int rr = 0;
988  for (; rr != nbRows / 3; ++rr) {
989 
990  // get sub matrix for the row
991  auto t_m = get_tensor2(K, 3 * rr, 0);
992 
993  // get derivatives of base functions for columns
994  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
995 
997  // I mix up the indices here so that it behaves like a
998  // Dg. That way I don't have to have a separate wrapper
999  // class Christof_Expr, which simplifies things.
1000  t_rowD(l, j, k) = t_D(i, j, k, l) * (a * t_row_diff_base(i));
1001 
1002  // iterate column base functions
1003  for (int cc = 0; cc != nbCols / 3; ++cc) {
1004 
1005  // integrate block local stiffens matrix
1006  t_m(i, j) += t_rowD(i, j, k) * t_col_diff_base(k);
1007 
1008  // move to next column base function
1009  ++t_col_diff_base;
1010 
1011  // move to next block of local stiffens matrix
1012  ++t_m;
1013  }
1014 
1015  // move to next row base function
1016  ++t_row_diff_base;
1017  }
1018 
1019  for (; rr != row_nb_base_fun; ++rr)
1020  ++t_row_diff_base;
1021 
1022  // move to next integration weight
1023  ++t_w;
1024  ++t_D;
1025  }
1026 
1028 }
1029 
1030 template <int S>
1031 HookeElement::OpAleLhs_dx_dX<S>::OpAleLhs_dx_dX(
1032  const std::string row_field, const std::string col_field,
1033  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1034  : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, false) {}
1035 
1036 template <int S>
1038  EntData &col_data) {
1040 
1041  // get sub-block (3x3) of local stiffens matrix, here represented by
1042  // second order tensor
1043  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1045  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1046  &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1047  &m(r + 2, c + 2));
1048  };
1049 
1056 
1057  // get element volume
1058  double vol = getVolume();
1059 
1060  // get intergrayion weights
1061  auto t_w = getFTensor0IntegrationWeight();
1062 
1063  // get derivatives of base functions on rows
1064  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1065  const int row_nb_base_fun = row_data.getN().size2();
1066 
1067  // Elastic stiffness tensor (4th rank tensor with minor and major
1068  // symmetry)
1070  MAT_TO_DDG(dataAtPts->stiffnessMat));
1071 
1072  auto t_cauchy_stress =
1073  getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1074  auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1075  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1076  auto &det_H = *dataAtPts->detHVec;
1077 
1078  // iterate over integration points
1079  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1080 
1081  // calculate scalar weight times element volume
1082  double a = t_w * vol * det_H[gg];
1083 
1085  t_F_dX(i, j, k, l) = -(t_h(i, m) * t_invH(m, k)) * t_invH(l, j);
1086 
1087  // iterate over row base functions
1088  int rr = 0;
1089  for (; rr != nbRows / 3; ++rr) {
1090 
1091  // get sub matrix for the row
1092  auto t_m = get_tensor2(K, 3 * rr, 0);
1093 
1094  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1095  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1096 
1097  FTensor::Tensor1<double, 3> t_row_stress;
1098  t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_cauchy_stress(i, j);
1099 
1100  FTensor::Tensor3<double, 3, 3, 3> t_row_diff_base_pulled_dX;
1101  t_row_diff_base_pulled_dX(j, k, l) =
1102  -(t_invH(i, k) * t_row_diff_base(i)) * t_invH(l, j);
1103 
1104  FTensor::Tensor3<double, 3, 3, 3> t_row_dX_stress;
1105  t_row_dX_stress(i, k, l) =
1106  a * (t_row_diff_base_pulled_dX(j, k, l) * t_cauchy_stress(j, i));
1107 
1109  t_row_D(l, j, k) = (a * t_row_diff_base_pulled(i)) * t_D(i, j, k, l);
1110 
1111  FTensor::Tensor3<double, 3, 3, 3> t_row_stress_dX;
1112  // FIXME: This operator is not implemented, doing operation by hand
1113  // t_row_stress_dX(i, m, n) = t_row_D(i, k, l) * t_F_dX(k, l, m, n);
1114  t_row_stress_dX(i, j, k) = 0;
1115  for (int ii = 0; ii != 3; ++ii)
1116  for (int mm = 0; mm != 3; ++mm)
1117  for (int nn = 0; nn != 3; ++nn) {
1118  auto &v = t_row_stress_dX(ii, mm, nn);
1119  for (int kk = 0; kk != 3; ++kk)
1120  for (int ll = 0; ll != 3; ++ll)
1121  v += t_row_D(ii, kk, ll) * t_F_dX(kk, ll, mm, nn);
1122  }
1123 
1124  // get derivatives of base functions for columns
1125  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1126 
1127  // iterate column base functions
1128  for (int cc = 0; cc != nbCols / 3; ++cc) {
1129 
1130  t_m(i, k) += t_row_stress(i) * (t_invH(j, k) * t_col_diff_base(j));
1131  t_m(i, k) += t_row_dX_stress(i, k, l) * t_col_diff_base(l);
1132  t_m(i, k) += t_row_stress_dX(i, k, l) * t_col_diff_base(l);
1133 
1134  // move to next column base function
1135  ++t_col_diff_base;
1136 
1137  // move to next block of local stiffens matrix
1138  ++t_m;
1139  }
1140 
1141  // move to next row base function
1142  ++t_row_diff_base;
1143  }
1144 
1145  for (; rr != row_nb_base_fun; ++rr)
1146  ++t_row_diff_base;
1147 
1148  // move to next integration weight
1149  ++t_w;
1150  ++t_D;
1151  ++t_cauchy_stress;
1152  ++t_invH;
1153  ++t_h;
1154  }
1155 
1157 }
1158 
1159 template <int S>
1160 HookeElement::OpAleLhs_dX_dX<S>::OpAleLhs_dX_dX(
1161  const std::string row_field, const std::string col_field,
1162  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1163  : OpAssemble(row_field, col_field, data_at_pts, OPROWCOL, true) {}
1164 
1165 template <int S>
1167  EntData &col_data) {
1169 
1170  // get sub-block (3x3) of local stiffens matrix, here represented by
1171  // second order tensor
1172  auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
1174  &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
1175  &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
1176  &m(r + 2, c + 2));
1177  };
1178 
1185 
1186  // get element volume
1187  double vol = getVolume();
1188 
1189  // get intergrayion weights
1190  auto t_w = getFTensor0IntegrationWeight();
1191 
1192  // get derivatives of base functions on rows
1193  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1194  const int row_nb_base_fun = row_data.getN().size2();
1195 
1196  // Elastic stiffness tensor (4th rank tensor with minor and major
1197  // symmetry)
1199  MAT_TO_DDG(dataAtPts->stiffnessMat));
1200  auto t_cauchy_stress =
1201  getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1202  auto t_strain = getFTensor2SymmetricFromMat<3>(*(dataAtPts->smallStrainMat));
1203  auto t_eshelby_stress =
1204  getFTensor2FromMat<3, 3>(*dataAtPts->eshelbyStressMat);
1205  auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1206  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1207  auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
1208  auto &det_H = *dataAtPts->detHVec;
1209 
1210  // iterate over integration points
1211  for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1212 
1213  // calculate scalar weight times element volume
1214  double a = t_w * vol * det_H[gg];
1215 
1217  t_F_dX(i, j, k, l) = -(t_h(i, m) * t_invH(m, k)) * t_invH(l, j);
1218 
1220  t_D_strain_dX(i, j, m, n) = 0.;
1221  for (int ii = 0; ii != 3; ++ii)
1222  for (int jj = 0; jj != 3; ++jj)
1223  for (int ll = 0; ll != 3; ++ll)
1224  for (int kk = 0; kk != 3; ++kk) {
1225  auto &v = t_D_strain_dX(ii, jj, kk, ll);
1226  for (int mm = 0; mm != 3; ++mm)
1227  for (int nn = 0; nn != 3; ++nn)
1228  v += t_D(ii, jj, mm, nn) * t_F_dX(mm, nn, kk, ll);
1229  }
1230 
1231  FTensor::Tensor4<double, 3, 3, 3, 3> t_eshelby_stress_dX;
1232  t_eshelby_stress_dX(i, j, m, n) = t_F(k, i) * t_D_strain_dX(k, j, m, n);
1233 
1234  for (int ii = 0; ii != 3; ++ii)
1235  for (int jj = 0; jj != 3; ++jj)
1236  for (int mm = 0; mm != 3; ++mm)
1237  for (int nn = 0; nn != 3; ++nn) {
1238  auto &v = t_eshelby_stress_dX(ii, jj, mm, nn);
1239  for (int kk = 0; kk != 3; ++kk)
1240  v += t_F_dX(kk, ii, mm, nn) * t_cauchy_stress(kk, jj);
1241  }
1242 
1243  t_eshelby_stress_dX(i, j, k, l) *= -1;
1244 
1245  FTensor::Tensor2<double, 3, 3> t_energy_dX;
1246  t_energy_dX(k, l) = t_F_dX(i, j, k, l) * t_cauchy_stress(i, j);
1247  t_energy_dX(k, l) +=
1248  (t_strain(m, n) * t_D(m, n, i, j)) * t_F_dX(i, j, k, l);
1249  t_energy_dX(k, l) /= 2.;
1250 
1251  for (int kk = 0; kk != 3; ++kk)
1252  for (int ll = 0; ll != 3; ++ll) {
1253  auto v = t_energy_dX(kk, ll);
1254  for (int ii = 0; ii != 3; ++ii)
1255  t_eshelby_stress_dX(ii, ii, kk, ll) += v;
1256  }
1257 
1258  // iterate over row base functions
1259  int rr = 0;
1260  for (; rr != nbRows / 3; ++rr) {
1261 
1262  // get sub matrix for the row
1263  auto t_m = get_tensor2(K, 3 * rr, 0);
1264 
1265  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1266  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1267 
1268  FTensor::Tensor1<double, 3> t_row_stress;
1269  t_row_stress(i) = a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
1270 
1271  FTensor::Tensor3<double, 3, 3, 3> t_row_diff_base_pulled_dX;
1272  t_row_diff_base_pulled_dX(j, k, l) =
1273  -(t_row_diff_base(i) * t_invH(i, k)) * t_invH(l, j);
1274 
1275  FTensor::Tensor3<double, 3, 3, 3> t_row_dX_stress;
1276  t_row_dX_stress(i, k, l) =
1277  a * (t_row_diff_base_pulled_dX(j, k, l) * t_eshelby_stress(i, j));
1278 
1279  FTensor::Tensor3<double, 3, 3, 3> t_row_stress_dX;
1280  t_row_stress_dX(i, m, n) =
1281  a * t_row_diff_base_pulled(j) * t_eshelby_stress_dX(i, j, m, n);
1282 
1283  // get derivatives of base functions for columns
1284  auto t_col_diff_base = col_data.getFTensor1DiffN<3>(gg, 0);
1285 
1286  // iterate column base functions
1287  for (int cc = 0; cc != nbCols / 3; ++cc) {
1288 
1289  t_m(i, k) += t_row_stress(i) * (t_invH(j, k) * t_col_diff_base(j));
1290  t_m(i, k) += t_row_dX_stress(i, k, l) * t_col_diff_base(l);
1291  t_m(i, k) += t_row_stress_dX(i, k, l) * t_col_diff_base(l);
1292 
1293  // move to next column base function
1294  ++t_col_diff_base;
1295 
1296  // move to next block of local stiffens matrix
1297  ++t_m;
1298  }
1299 
1300  // move to next row base function
1301  ++t_row_diff_base;
1302  }
1303 
1304  for (; rr != row_nb_base_fun; ++rr)
1305  ++t_row_diff_base;
1306 
1307  // move to next integration weight
1308  ++t_w;
1309  ++t_D;
1310  ++t_cauchy_stress;
1311  ++t_strain;
1312  ++t_eshelby_stress;
1313  ++t_h;
1314  ++t_invH;
1315  ++t_F;
1316  }
1317 
1319 }
1320 
1321 template <int S>
1322 HookeElement::OpAleLhsPre_dX_dx<S>::OpAleLhsPre_dX_dx(
1323  const std::string row_field, const std::string col_field,
1324  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts)
1325  : VolUserDataOperator(row_field, col_field, OPROW, true),
1326  dataAtPts(data_at_pts) {
1327  std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE], false);
1328 }
1329 
1330 template <int S>
1331 MoFEMErrorCode HookeElement::OpAleLhsPre_dX_dx<S>::doWork(int row_side,
1332  EntityType row_type,
1333  EntData &row_data) {
1335 
1336  const int nb_integration_pts = row_data.getN().size1();
1337 
1338  auto get_eshelby_stress_dx = [this, nb_integration_pts]() {
1340  t_eshelby_stress_dx;
1341  dataAtPts->eshelbyStress_dx->resize(81, nb_integration_pts, false);
1342  int mm = 0;
1343  for (int ii = 0; ii != 3; ++ii)
1344  for (int jj = 0; jj != 3; ++jj)
1345  for (int kk = 0; kk != 3; ++kk)
1346  for (int ll = 0; ll != 3; ++ll)
1347  t_eshelby_stress_dx.ptr(ii, jj, kk, ll) =
1348  &(*dataAtPts->eshelbyStress_dx)(mm++, 0);
1349  return t_eshelby_stress_dx;
1350  };
1351 
1352  auto t_eshelby_stress_dx = get_eshelby_stress_dx();
1353 
1360 
1361  // Elastic stiffness tensor (4th rank tensor with minor and major
1362  // symmetry)
1364  MAT_TO_DDG(dataAtPts->stiffnessMat));
1365  auto t_cauchy_stress =
1366  getFTensor2SymmetricFromMat<3>(*(dataAtPts->cauchyStressMat));
1367  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1368  auto t_F = getFTensor2FromMat<3, 3>(*dataAtPts->FMat);
1369 
1370  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1371 
1372  t_eshelby_stress_dx(i, j, m, n) =
1373  (t_F(k, i) * t_D(k, j, m, l)) * t_invH(n, l);
1374  for (int ii = 0; ii != 3; ++ii)
1375  for (int jj = 0; jj != 3; ++jj)
1376  for (int mm = 0; mm != 3; ++mm)
1377  for (int nn = 0; nn != 3; ++nn) {
1378  auto &v = t_eshelby_stress_dx(ii, jj, mm, nn);
1379  v += t_invH(nn, ii) * t_cauchy_stress(mm, jj);
1380  }
1381  t_eshelby_stress_dx(i, j, k, l) *= -1;
1382 
1383  FTensor::Tensor2<double, 3, 3> t_energy_dx;
1384  t_energy_dx(m, n) = t_invH(n, j) * t_cauchy_stress(m, j);
1385 
1386  for (int mm = 0; mm != 3; ++mm)
1387  for (int nn = 0; nn != 3; ++nn) {
1388  auto v = t_energy_dx(mm, nn);
1389  for (int ii = 0; ii != 3; ++ii)
1390  t_eshelby_stress_dx(ii, ii, mm, nn) += v;
1391  }
1392 
1393  ++t_D;
1394  ++t_invH;
1395  ++t_cauchy_stress;
1396  ++t_eshelby_stress_dx;
1397  ++t_F;
1398  }
1399 
1401 }
1402 
1403 template <class ELEMENT>
1404 HookeElement::OpPostProcHookeElement<ELEMENT>::OpPostProcHookeElement(
1405  const string row_field, boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1406  map<int, BlockData> &block_sets_ptr, moab::Interface &post_proc_mesh,
1407  std::vector<EntityHandle> &map_gauss_pts, bool is_ale, bool is_field_disp)
1408  : ELEMENT::UserDataOperator(row_field, UserDataOperator::OPROW),
1409  dataAtPts(data_at_pts), blockSetsPtr(block_sets_ptr),
1410  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts), isALE(is_ale),
1411  isFieldDisp(is_field_disp) {}
1412 
1413 template <class ELEMENT>
1414 MoFEMErrorCode HookeElement::OpPostProcHookeElement<ELEMENT>::doWork(
1415  int side, EntityType type, DataForcesAndSourcesCore::EntData &data) {
1417 
1418  if (type != MBVERTEX) {
1420  }
1421 
1422  auto tensor_to_tensor = [](const auto &t1, auto &t2) {
1423  t2(0, 0) = t1(0, 0);
1424  t2(1, 1) = t1(1, 1);
1425  t2(2, 2) = t1(2, 2);
1426  t2(0, 1) = t2(1, 0) = t1(1, 0);
1427  t2(0, 2) = t2(2, 0) = t1(2, 0);
1428  t2(1, 2) = t2(2, 1) = t1(2, 1);
1429  };
1430 
1431  std::array<double,9> def_val;
1432  def_val.fill(0);
1433 
1434  auto make_tag = [&](auto name, auto size) {
1435  Tag th;
1436  CHKERR postProcMesh.tag_get_handle(name, size, MB_TYPE_DOUBLE, th,
1437  MB_TAG_CREAT | MB_TAG_SPARSE,
1438  def_val.data());
1439  return th;
1440  };
1441 
1442  auto th_stress = make_tag("STRESS", 9);
1443  auto th_psi = make_tag("ENERGY", 1);
1444 
1445  const int nb_integration_pts = mapGaussPts.size();
1446 
1451 
1452  auto t_h = getFTensor2FromMat<3, 3>(*dataAtPts->hMat);
1453  auto t_H = getFTensor2FromMat<3, 3>(*dataAtPts->HMat);
1454 
1455  dataAtPts->stiffnessMat->resize(36, 1, false);
1457  MAT_TO_DDG(dataAtPts->stiffnessMat));
1458  for (auto &m : (blockSetsPtr)) {
1459  const double young = m.second.E;
1460  const double poisson = m.second.PoissonRatio;
1461 
1462  const double coefficient = young / ((1 + poisson) * (1 - 2 * poisson));
1463 
1464  t_D(i, j, k, l) = 0.;
1465  t_D(0, 0, 0, 0) = t_D(1, 1, 1, 1) = t_D(2, 2, 2, 2) = 1 - poisson;
1466  t_D(0, 1, 0, 1) = t_D(0, 2, 0, 2) = t_D(1, 2, 1, 2) =
1467  0.5 * (1 - 2 * poisson);
1468  t_D(0, 0, 1, 1) = t_D(1, 1, 0, 0) = t_D(0, 0, 2, 2) = t_D(2, 2, 0, 0) =
1469  t_D(1, 1, 2, 2) = t_D(2, 2, 1, 1) = poisson;
1470  t_D(i, j, k, l) *= coefficient;
1471 
1472  break; // FIXME: calculates only first block
1473  }
1474 
1475  double detH = 0.;
1479  FTensor::Tensor2<double, 3, 3> t_small_strain;
1481  FTensor::Tensor2_symmetric<double, 3> t_small_strain_symm;
1482 
1483  for (int gg = 0; gg != nb_integration_pts; ++gg) {
1484 
1485  if (isFieldDisp) {
1486  t_h(0, 0) += 1;
1487  t_h(1, 1) += 1;
1488  t_h(2, 2) += 1;
1489  }
1490 
1491  if (!isALE) {
1492  t_small_strain_symm(i, j) = (t_h(i, j) || t_h(j, i)) / 2.;
1493  } else {
1494  CHKERR determinantTensor3by3(t_H, detH);
1495  CHKERR invertTensor3by3(t_H, detH, t_invH);
1496  t_F(i, j) = t_h(i, k) * t_invH(k, j);
1497  t_small_strain_symm(i, j) = (t_F(i, j) || t_F(j, i)) / 2.;
1498  ++t_H;
1499  }
1500 
1501  t_small_strain_symm(0, 0) -= 1;
1502  t_small_strain_symm(1, 1) -= 1;
1503  t_small_strain_symm(2, 2) -= 1;
1504 
1505  // symmetric tensors need improvement
1506  t_stress_symm(i, j) = t_D(i, j, k, l) * t_small_strain_symm(k, l);
1507  tensor_to_tensor(t_stress_symm, t_stress);
1508 
1509  const double psi = 0.5 * t_stress_symm(i, j) * t_small_strain_symm(i, j);
1510 
1511  CHKERR postProcMesh.tag_set_data(th_psi, &mapGaussPts[gg], 1, &psi);
1512  CHKERR postProcMesh.tag_set_data(th_stress, &mapGaussPts[gg], 1,
1513  &t_stress(0, 0));
1514 
1515  ++t_h;
1516  }
1517 
1519 }
1520 
1521 template <int S>
1522 HookeElement::OpAnalyticalInternalStain_dx<S>::OpAnalyticalInternalStain_dx(
1523  const std::string row_field,
1524  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1525  StrainFunctions strain_fun)
1526  : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1527  strainFun(strain_fun) {}
1528 
1529 template <int S>
1537 
1538  auto get_tensor1 = [](VectorDouble &v, const int r) {
1540  &v(r + 0), &v(r + 1), &v(r + 2));
1541  };
1542 
1543  const int nb_integration_pts = getGaussPts().size2();
1544 
1545  auto get_coords = [&]() {
1546  if (getHoCoordsAtGaussPts().size1() == nb_integration_pts)
1547  return getFTensor1HoCoordsAtGaussPts();
1548  else
1549  return getFTensor1CoordsAtGaussPts();
1550  };
1551  auto t_coords = get_coords();
1552 
1553  // get element volume
1554  double vol = getVolume();
1555  auto t_w = getFTensor0IntegrationWeight();
1556 
1557  nF.resize(nbRows, false);
1558  nF.clear();
1559 
1560  // elastic stiffness tensor (4th rank tensor with minor and major
1561  // symmetry)
1563  MAT_TO_DDG(dataAtPts->stiffnessMat));
1564 
1565  // get derivatives of base functions on rows
1566  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1567  const int row_nb_base_fun = row_data.getN().size2();
1568 
1569  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1570 
1571  auto t_fun_strain = strainFun(t_coords);
1573  t_stress(i, j) = t_D(i, j, k, l) * t_fun_strain(k, l);
1574 
1575  // calculate scalar weight times element volume
1576  double a = t_w * vol;
1577 
1578  if (getHoGaussPtsDetJac().size()) {
1579  // If HO geometry
1580  a *= getHoGaussPtsDetJac()[gg];
1581  }
1582 
1583  auto t_nf = get_tensor1(nF, 0);
1584 
1585  int rr = 0;
1586  for (; rr != nbRows / 3; ++rr) {
1587  t_nf(i) += a * t_row_diff_base(j) * t_stress(i, j);
1588  ++t_row_diff_base;
1589  ++t_nf;
1590  }
1591 
1592  for (; rr != row_nb_base_fun; ++rr)
1593  ++t_row_diff_base;
1594 
1595  ++t_w;
1596  ++t_coords;
1597  ++t_D;
1598  }
1599 
1601 }
1602 
1603 template <int S>
1604 HookeElement::OpAnalyticalInternalAleStain_dX<S>::
1605  OpAnalyticalInternalAleStain_dX(
1606  const std::string row_field,
1607  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1608  StrainFunctions strain_fun,
1609  boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1610  : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1611  strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1612 
1613 template <int S>
1621 
1622  auto get_tensor1 = [](VectorDouble &v, const int r) {
1624  &v(r + 0), &v(r + 1), &v(r + 2));
1625  };
1626 
1627  const int nb_integration_pts = getGaussPts().size2();
1628 
1629  auto get_coords = [&]() { return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
1630  auto t_coords = get_coords();
1631 
1632  // get element volume
1633  double vol = getVolume();
1634  auto t_w = getFTensor0IntegrationWeight();
1635 
1636  nF.resize(nbRows, false);
1637  nF.clear();
1638 
1639  // elastic stiffness tensor (4th rank tensor with minor and major
1640  // symmetry)
1642  MAT_TO_DDG(dataAtPts->stiffnessMat));
1643  auto t_F = getFTensor2FromMat<3, 3>(*(dataAtPts->FMat));
1644  auto &det_H = *dataAtPts->detHVec;
1645  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1646 
1647  // get derivatives of base functions on rows
1648  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1649  const int row_nb_base_fun = row_data.getN().size2();
1650 
1651  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1652 
1653  auto t_fun_strain = strainFun(t_coords);
1655  t_stress(i, j) = t_D(i, j, k, l) * t_fun_strain(k, l);
1656  FTensor::Tensor2<double, 3, 3> t_eshelby_stress;
1657  t_eshelby_stress(i, j) = -t_F(k, i) * t_stress(k, j);
1658 
1659  // calculate scalar weight times element volume
1660  double a = t_w * vol * det_H[gg];
1661 
1662  auto t_nf = get_tensor1(nF, 0);
1663 
1664  int rr = 0;
1665  for (; rr != nbRows / 3; ++rr) {
1666  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1667  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1668  t_nf(i) += a * t_row_diff_base_pulled(j) * t_eshelby_stress(i, j);
1669  ++t_row_diff_base;
1670  ++t_nf;
1671  }
1672 
1673  for (; rr != row_nb_base_fun; ++rr)
1674  ++t_row_diff_base;
1675 
1676  ++t_w;
1677  ++t_coords;
1678  ++t_F;
1679  ++t_invH;
1680  ++t_D;
1681  }
1682 
1684 }
1685 
1686 template <int S>
1687 HookeElement::OpAnalyticalInternalAleStain_dx<S>::
1688  OpAnalyticalInternalAleStain_dx(
1689  const std::string row_field,
1690  boost::shared_ptr<DataAtIntegrationPts> &data_at_pts,
1691  StrainFunctions strain_fun,
1692  boost::shared_ptr<MatrixDouble> mat_pos_at_pts_ptr)
1693  : OpAssemble(row_field, row_field, data_at_pts, OPROW, true),
1694  strainFun(strain_fun), matPosAtPtsPtr(mat_pos_at_pts_ptr) {}
1695 
1696 template <int S>
1704 
1705  auto get_tensor1 = [](VectorDouble &v, const int r) {
1707  &v(r + 0), &v(r + 1), &v(r + 2));
1708  };
1709 
1710  const int nb_integration_pts = getGaussPts().size2();
1711 
1712  auto get_coords = [&]() { return getFTensor1FromMat<3>(*matPosAtPtsPtr); };
1713  auto t_coords = get_coords();
1714 
1715  // get element volume
1716  double vol = getVolume();
1717  auto t_w = getFTensor0IntegrationWeight();
1718 
1719  nF.resize(nbRows, false);
1720  nF.clear();
1721 
1722  // elastic stiffness tensor (4th rank tensor with minor and major
1723  // symmetry)
1725  MAT_TO_DDG(dataAtPts->stiffnessMat));
1726  auto &det_H = *dataAtPts->detHVec;
1727  auto t_invH = getFTensor2FromMat<3, 3>(*dataAtPts->invHMat);
1728 
1729  // get derivatives of base functions on rows
1730  auto t_row_diff_base = row_data.getFTensor1DiffN<3>();
1731  const int row_nb_base_fun = row_data.getN().size2();
1732 
1733  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1734 
1735  auto t_fun_strain = strainFun(t_coords);
1737  t_stress(i, j) = t_D(i, j, k, l) * t_fun_strain(k, l);
1738 
1739  // calculate scalar weight times element volume
1740  double a = t_w * vol * det_H[gg];
1741 
1742  auto t_nf = get_tensor1(nF, 0);
1743 
1744  int rr = 0;
1745  for (; rr != nbRows / 3; ++rr) {
1746  FTensor::Tensor1<double, 3> t_row_diff_base_pulled;
1747  t_row_diff_base_pulled(i) = t_row_diff_base(j) * t_invH(j, i);
1748  t_nf(i) += a * t_row_diff_base_pulled(j) * t_stress(i, j);
1749  ++t_row_diff_base;
1750  ++t_nf;
1751  }
1752 
1753  for (; rr != row_nb_base_fun; ++rr)
1754  ++t_row_diff_base;
1755 
1756  ++t_w;
1757  ++t_coords;
1758  ++t_invH;
1759  ++t_D;
1760  }
1761 
1763 }
1764 
1765 #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:571
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< MatrixDouble > eshelbyStress_dx
DataForcesAndSourcesCore::EntData EntData
int nbCols
number if dof on column
const double D
diffusivity
boost::shared_ptr< map< int, BlockData > > blockSetsPtr
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
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:485
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:516
static Index< 'l', 3 > l
bool isDiag
true if this block is on diagonal
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< MatrixDouble > invHMat
static Index< 'n', 3 > n
MoFEMErrorCode invertTensor3by3(ublas::matrix< T, L, A > &jac_data, ublas::vector< T, A > &det_data, ublas::matrix< T, L, A > &inv_jac_data)
Calculate inverse of tensor rank 2 at integration points.
Definition: Templates.hpp:552
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
static Index< 'm', 3 > m
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
static Index< 'i', 3 > i
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
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
VectorDouble a0
constant acceleration
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
static Index< 'j', 3 > j
boost::shared_ptr< MatrixDouble > eshelbyStressMat
boost::shared_ptr< VectorDouble > detHVec
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:604
const double r
rate factor
boost::function< FTensor::Tensor2_symmetric< double, 3 > FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords) > StrainFunctions
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:73
bool isFieldDisp
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)
static Index< 'k', 3 > k
boost::shared_ptr< VectorDouble > energyVec
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
struct ConvectiveMassElement BlockData
Data on single entity (This is passed as argument to DataOperator::doWork)
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:415
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:74
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
Index< 'K', 3 > K
MassBlockData & massData
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:169
structure grouping operators and data used for calculation of nonlinear elastic elementIn order to as...
boost::function< FTensor::Tensor2_symmetric< double, 3 > FTensor::Tensor1< FTensor::PackPtr< double *, 1 >, 3 > &t_coords) > StrainFunctions