v0.8.23
EshelbianPlasticity.hpp
Go to the documentation of this file.
1 /**
2  * \file EshelbianPlasticity.hpp
3  * \brief Eshelbian plasticity interface
4  *
5  * \brief Problem implementation for mix element for large-strain elasticity
6  *
7  * \todo Implementation of plasticity
8  */
9 
10 #ifndef __ESHELBIAN_PLASTICITY_HPP__
11 #define __ESHELBIAN_PLASTICITY_HPP__
12 
13 namespace EshelbianPlasticity {
14 
16 
17 static const MOFEMuuid IDD_MOFEMEshelbianCrackInterface =
19 
20 typedef boost::shared_ptr<MatrixDouble> MatrixPtr;
21 typedef boost::shared_ptr<VectorDouble> VectorPtr;
22 
26  &m(0, 0), &m(1, 0), &m(2, 0), &m(3, 0), &m(4, 0), &m(5, 0), &m(6, 0),
27  &m(7, 0), &m(8, 0), &m(9, 0), &m(10, 0), &m(11, 0), &m(12, 0), &m(13, 0),
28  &m(14, 0), &m(15, 0), &m(16, 0), &m(17, 0), &m(18, 0), &m(19, 0),
29  &m(20, 0), &m(21, 0), &m(22, 0), &m(23, 0), &m(24, 0), &m(25, 0),
30  &m(26, 0));
31 }
32 
33 template <typename T>
40  t_R(i, j) = 0;
41  t_R(0, 0) = 1;
42  t_R(1, 1) = 1;
43  t_R(2, 2) = 1;
44 
45  const double angle = sqrt(t_omega(i) * t_omega(i));
46 
47  const double tol = 1e-12;
48  if (fabs(angle) < tol) {
49  return t_R;
50  }
51 
53  t_Omega(i, j) = levi_civita(i, j, k) * t_omega(k);
54  const double a = sin(angle) / angle;
55  const double ss_2 = sin(angle / 2.);
56  const double b = 2 * ss_2 * ss_2 / (angle * angle);
57  t_R(i, j) += a * t_Omega(i, j);
58  t_R(i, j) += b * t_Omega(i, k) * t_Omega(k, j);
59 
60  return t_R;
61 }
62 
63 template <typename T>
66 
71 
72  double angle = sqrt(t_omega(i) * t_omega(i));
73  const double tol = 1e-12;
74  if (fabs(angle) < tol) {
76  auto t_R = getRotationFormVector(t_omega);
77  t_diff_R(i, j, k) = levi_civita(i, l, k) * t_R(l, j);
78  return t_diff_R;
79  }
80 
83  t_Omega(i, j) = levi_civita(i, j, k) * t_omega(k);
84 
85  const double angle2 = angle * angle;
86  const double ss = sin(angle);
87  const double a = ss / angle;
88  const double ss_2 = sin(angle / 2.);
89  const double b = 2 * ss_2 * ss_2 / angle2;
90  t_diff_R(i, j, k) = 0;
91  t_diff_R(i, j, k) = a * levi_civita(i, j, k);
92  t_diff_R(i, j, k) += b * (t_Omega(i, l) * levi_civita(l, j, k) +
93  levi_civita(i, l, k) * t_Omega(l, j));
94 
95  const double cc = cos(angle);
96  const double cc_2 = cos(angle / 2.);
97  const double diff_a = (angle * cc - ss) / angle2;
98 
99  const double diff_b =
100  (2. * angle * ss_2 * cc_2 - 4. * ss_2 * ss_2) / (angle2 * angle);
101  t_diff_R(i, j, k) += diff_a * t_Omega(i, j) * (t_omega(k) / angle);
102  t_diff_R(i, j, k) +=
103  diff_b * t_Omega(i, l) * t_Omega(l, j) * (t_omega(k) / angle);
104 
105  return t_diff_R;
106 }
107 
112 
114  Mat Suu;
115  AO aoSuu;
116  Mat SBubble;
118  Mat SOmega;
120  Mat Sw;
121  AO aoSw;
123  : Suu(PETSC_NULL), aoSuu(PETSC_NULL), SBubble(PETSC_NULL),
124  aoSBubble(PETSC_NULL), Sw(PETSC_NULL), aoSw(PETSC_NULL),
125  SOmega(PETSC_NULL), aoSOmega(PETSC_NULL) {}
127  if (Suu)
128  MatDestroy(&Suu);
129  if (aoSuu)
130  AODestroy(&aoSuu);
131  if (SBubble)
132  MatDestroy(&SBubble);
133  if (aoSBubble)
134  AODestroy(&aoSBubble);
135  if (Sw)
136  MatDestroy(&Sw);
137  if (aoSw)
138  AODestroy(&aoSw);
139  if (SOmega)
140  MatDestroy(&SOmega);
141  if (aoSOmega)
142  AODestroy(&aoSOmega);
143  }
144 
147  if (this->Suu)
148  MatDestroy(&this->Suu);
149  if (this->aoSuu)
150  AODestroy(&this->aoSuu);
151  this->Suu = Suu;
152  this->aoSuu = aoSuu;
153  PetscObjectReference((PetscObject)this->Suu);
154  PetscObjectReference((PetscObject)this->aoSuu);
156  }
157 
160  if (this->SBubble)
161  MatDestroy(&this->SBubble);
162  if (this->aoSBubble)
163  AODestroy(&this->aoSBubble);
164  this->SBubble = SBubble;
165  this->aoSBubble = aoSBubble;
166  PetscObjectReference((PetscObject)this->SBubble);
167  PetscObjectReference((PetscObject)this->aoSBubble);
169  }
170 
173  if (this->Sw)
174  MatDestroy(&this->Sw);
175  if (this->aoSw)
176  AODestroy(&this->aoSw);
177  this->Sw = Sw;
178  this->aoSw = aoSw;
179  PetscObjectReference((PetscObject)this->Sw);
180  PetscObjectReference((PetscObject)this->aoSw);
182  }
183 
186  if (this->SOmega)
187  MatDestroy(&this->SOmega);
188  if (this->aoSOmega)
189  AODestroy(&this->aoSOmega);
190  this->SOmega = SOmega;
191  this->aoSOmega = aoSOmega;
192  PetscObjectReference((PetscObject)this->SOmega);
193  PetscObjectReference((PetscObject)this->aoSOmega);
195  }
196 };
197 
198 template <typename E> struct EpElement : public E, EpElementBase {
199  EpElement(MoFEM::Interface &m_field) : E(m_field), EpElementBase() {}
200 };
201 
202 template <> struct EpElement<BasicMethod> : public FEMethod, EpElementBase {
204 };
205 
206 template <> struct EpElement<FEMethod> : public FEMethod, EpElementBase {
208 };
209 
210 struct EpFEMethod : EpElement<FEMethod> {
214  if (Suu)
215  CHKERR MatZeroEntries(Suu);
216  if (SBubble)
217  CHKERR MatZeroEntries(SBubble);
218  if (Sw)
219  CHKERR MatZeroEntries(Sw);
220  if (SOmega)
221  CHKERR MatZeroEntries(SOmega);
223  }
224 
227  auto assemble = [](Mat &a) {
229  if (a) {
230  CHKERR MatAssemblyBegin(a, MAT_FINAL_ASSEMBLY);
231  CHKERR MatAssemblyEnd(a, MAT_FINAL_ASSEMBLY);
232  }
234  };
235  CHKERR assemble(Suu);
236  CHKERR assemble(SBubble);
237  CHKERR assemble(SOmega);
238  CHKERR assemble(Sw);
239  // std::string wait;
240  // CHKERR MatView(SOmega, PETSC_VIEWER_DRAW_WORLD);
241  // std::cin >> wait;
242  // CHKERR MatView(Sw, PETSC_VIEWER_DRAW_WORLD);
243  // std::cin >> wait;
244  // CHKERR MatView(SBubble, PETSC_VIEWER_DRAW_WORLD);
245  // std::cin >> wait;
246  // CHKERR MatView(Suu, PETSC_VIEWER_DRAW_WORLD);
247  // std::cin >> wait;
248  // CHKERR MatView(ts_B, PETSC_VIEWER_DRAW_WORLD);
249  // std::cin >> wait;
251  }
252 };
253 
255 
276 
286 
307 
308  // Not really data at integration points, used to calculate Schur complement
309 
310  struct BlockMatData {
311 
312  std::string rowField;
313  std::string colField;
314  EntityType rowType;
315  EntityType colType;
316  int rowSide;
317  int colSide;
321 
323 
324  BlockMatData(const std::string row_field, const std::string col_field,
325  EntityType row_type, EntityType col_type, int row_side,
326  int col_side, const MatrixDouble &m, const VectorInt row_ind,
327  VectorInt col_ind)
328  : rowField(row_field), colField(col_field), rowType(row_type),
329  colType(col_type), rowSide(row_side), colSide(col_side),
330  setAtElement(true) {
331 
332  M.resize(m.size1(), m.size2(), false);
333  noalias(M) = m;
334  rowInd.resize(row_ind.size(), false);
335  noalias(rowInd) = row_ind;
336  colInd.resize(col_ind.size(), false);
337  noalias(colInd) = col_ind;
338 
339  }
340 
341  void setInd(const VectorInt &row_ind, const VectorInt &col_ind) const {
342  auto &const_row_ind = const_cast<VectorInt &>(rowInd);
343  auto &const_col_ind = const_cast<VectorInt &>(colInd);
344  const_row_ind.resize(row_ind.size(), false);
345  noalias(const_row_ind) = row_ind;
346  const_col_ind.resize(col_ind.size(), false);
347  noalias(const_col_ind) = col_ind;
348  }
349 
350  void setMat(const MatrixDouble &m) const {
351  auto &const_m = const_cast<MatrixDouble &>(M);
352  const_m.resize(m.size1(), m.size2(), false);
353  noalias(const_m) = m;
354  }
355 
356  void addMat(const MatrixDouble &m) const {
357  auto &const_m = const_cast<MatrixDouble &>(M);
358  const_m += m;
359  }
360 
361  void clearMat() const {
362  auto &const_m = const_cast<MatrixDouble &>(M);
363  const_m.clear();
364  }
365 
366  void setSetAtElement() const {
367  bool &set = const_cast<bool &>(setAtElement);
368  set = true;
369  }
370 
371  void unSetAtElement() const {
372  bool &set = const_cast<bool &>(setAtElement);
373  set = false;
374  }
375  };
376 
377  typedef multi_index_container<
378  BlockMatData,
379  indexed_by<
380 
381  ordered_unique<
382 
383  composite_key<
384  BlockMatData,
385 
386  member<BlockMatData, std::string, &BlockMatData::rowField>,
387  member<BlockMatData, std::string, &BlockMatData::colField>,
388  member<BlockMatData, EntityType, &BlockMatData::rowType>,
389  member<BlockMatData, EntityType, &BlockMatData::colType>,
390  member<BlockMatData, int, &BlockMatData::rowSide>,
391  member<BlockMatData, int, &BlockMatData::colSide>
392 
393  >>,
394  ordered_non_unique<
395 
396  composite_key<
397  BlockMatData,
398 
399  member<BlockMatData, std::string, &BlockMatData::rowField>,
400  member<BlockMatData, std::string, &BlockMatData::colField>,
401  member<BlockMatData, EntityType, &BlockMatData::rowType>,
402  member<BlockMatData, EntityType, &BlockMatData::colType>
403 
404  >>,
405  ordered_non_unique<
406 
407  composite_key<
408  BlockMatData,
409 
410  member<BlockMatData, std::string, &BlockMatData::rowField>,
411  member<BlockMatData, std::string, &BlockMatData::colField>
412 
413  >>,
414  ordered_non_unique<
415  member<BlockMatData, std::string, &BlockMatData::rowField>>,
416  ordered_non_unique<
417  member<BlockMatData, std::string, &BlockMatData::colField>>
418 
419  >>
421 
425 };
426 
427 struct OpJacobian;
428 
430 
435 
438 
439  PhysicalEquations() = delete;
440  PhysicalEquations(const int size_active, const int size_dependent)
441  : activeVariables(size_active, 0), dependentVariables(size_dependent, 0),
442  dependentVariablesDirevatives(size_dependent * size_active, 0) {}
443  virtual ~PhysicalEquations() {}
444 
445  virtual MoFEMErrorCode recordTape(const int tag) = 0;
446 
447  virtual OpJacobian *
448  returnOpJacobian(const std::string &field_name, const int tag,
449  const bool eval_rhs, const bool eval_lhs,
450  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
451  boost::shared_ptr<PhysicalEquations> &physics_ptr) = 0;
452 
453  std::vector<double> activeVariables;
454  std::vector<double> dependentVariables;
455  std::vector<double> dependentVariablesDirevatives;
456 
457  /** \name Active variables */
458 
459  /**@{*/
460 
461  template <int S>
462  inline static DTensor2Ptr get_VecTensor2(std::vector<double> &v) {
463  return DTensor2Ptr(&v[S + 0], &v[S + 1], &v[S + 2], &v[S + 3], &v[S + 4],
464  &v[S + 5], &v[S + 6], &v[S + 7], &v[S + 8]);
465  }
466 
467  template <int S>
468  inline static DTensor0Ptr get_VecTensor0(std::vector<double> &v) {
469  return DTensor0Ptr(&v[S + 0]);
470  }
471 
472  template <int S0, int S1, int S2>
473  inline static DTensor3Ptr get_vecTensor3(std::vector<double> &v) {
474  constexpr int A00 = 18 * S0 + 18 * 0 + 9 * S1 + 3 * S2;
475  constexpr int A01 = 18 * S0 + 18 * 1 + 9 * S1 + 3 * S2;
476  constexpr int A02 = 18 * S0 + 18 * 2 + 9 * S1 + 3 * S2;
477  constexpr int A10 = 18 * S0 + 18 * 3 + 9 * S1 + 3 * S2;
478  constexpr int A11 = 18 * S0 + 18 * 4 + 9 * S1 + 3 * S2;
479  constexpr int A12 = 18 * S0 + 18 * 5 + 9 * S1 + 3 * S2;
480  constexpr int A20 = 18 * S0 + 18 * 6 + 9 * S1 + 3 * S2;
481  constexpr int A21 = 18 * S0 + 18 * 7 + 9 * S1 + 3 * S2;
482  constexpr int A22 = 18 * S0 + 18 * 8 + 9 * S1 + 3 * S2;
483  return DTensor3Ptr(
484 
485  &v[A00 + 0], &v[A00 + 1], &v[A00 + 2], &v[A01 + 0], &v[A01 + 1],
486  &v[A01 + 2], &v[A02 + 0], &v[A02 + 1], &v[A02 + 2],
487 
488  &v[A10 + 0], &v[A10 + 1], &v[A10 + 2], &v[A11 + 0], &v[A11 + 1],
489  &v[A11 + 2], &v[A12 + 0], &v[A12 + 1], &v[A12 + 2],
490 
491  &v[A20 + 0], &v[A20 + 1], &v[A20 + 2], &v[A21 + 0], &v[A21 + 1],
492  &v[A21 + 2], &v[A22 + 0], &v[A22 + 1], &v[A22 + 2]
493 
494  );
495  }
496 
497  /**@}*/
498 
499  /** \name Active variables */
500 
501  /**@{*/
502 
503  inline DTensor2Ptr get_h() { return get_VecTensor2<0>(activeVariables); }
504  inline DTensor2Ptr get_H() { return get_VecTensor2<9>(activeVariables); }
505 
506  /**@}*/
507 
508  /** \name Dependent variables */
509 
510  /**@{*/
511 
512  inline DTensor2Ptr get_P() { return get_VecTensor2<0>(dependentVariables); }
514  return get_VecTensor2<9>(dependentVariables);
515  }
516  inline DTensor0Ptr get_Phi() {
517  return get_VecTensor0<18>(dependentVariables);
518  }
519  inline double &get_PhiRef() { return dependentVariables[18]; }
521  return get_VecTensor2<9 + 9 + 1>(dependentVariables);
522  }
523 
524  /**@}*/
525 
526  /** \name Derivatives of dependent variables */
527 
528  /**@{*/
529 
531  return get_vecTensor3<0, 0, 0>(dependentVariablesDirevatives);
532  }
534  return get_vecTensor3<0, 1, 0>(dependentVariablesDirevatives);
535  }
537  return get_vecTensor3<0, 0, 1>(dependentVariablesDirevatives);
538  }
540  return get_vecTensor3<0, 1, 1>(dependentVariablesDirevatives);
541  }
543  return get_vecTensor3<0, 0, 2>(dependentVariablesDirevatives);
544  }
546  return get_vecTensor3<0, 1, 2>(dependentVariablesDirevatives);
547  }
548 
550  return get_vecTensor3<9, 0, 0>(dependentVariablesDirevatives);
551  }
553  return get_vecTensor3<9, 1, 0>(dependentVariablesDirevatives);
554  }
556  return get_vecTensor3<9, 0, 1>(dependentVariablesDirevatives);
557  }
559  return get_vecTensor3<9, 1, 1>(dependentVariablesDirevatives);
560  }
562  return get_vecTensor3<9, 0, 2>(dependentVariablesDirevatives);
563  }
565  return get_vecTensor3<9, 1, 2>(dependentVariablesDirevatives);
566  }
567 
569  return get_VecTensor2<18 * (9 + 9) + 0>(dependentVariablesDirevatives);
570  }
572  return get_VecTensor2<18 * (9 + 9) + 9>(dependentVariablesDirevatives);
573  }
574 
576  return get_vecTensor3<9 + 9 + 1, 0, 0>(dependentVariablesDirevatives);
577  }
579  return get_vecTensor3<9 + 9 + 1, 1, 0>(dependentVariablesDirevatives);
580  }
582  return get_vecTensor3<9 + 9 + 1, 0, 1>(dependentVariablesDirevatives);
583  }
585  return get_vecTensor3<9 + 9 + 1, 1, 1>(dependentVariablesDirevatives);
586  }
588  return get_vecTensor3<9 + 9 + 1, 0, 2>(dependentVariablesDirevatives);
589  }
591  return get_vecTensor3<9 + 9 + 1, 1, 2>(dependentVariablesDirevatives);
592  }
593 
594  /**@}*/
595 };
596 
597 struct BcDisp {
598  BcDisp(std::string name, std::vector<double> &attr, Range &faces);
599  std::string blockName;
600  Range faces;
603 };
604 typedef std::vector<BcDisp> BcDispVec;
605 
606 struct BcRot {
607  BcRot(std::string name, std::vector<double> &attr, Range &faces);
608  std::string blockName;
609  Range faces;
611  double theta;
612 };
613 typedef std::vector<BcRot> BcRotVec;
614 
615 typedef std::vector<Range> TractionFreeBc;
616 
617 struct TractionBc {
618  TractionBc(std::string name, std::vector<double> &attr, Range &faces);
619  std::string blockName;
620  Range faces;
623 };
624 typedef std::vector<TractionBc> TractionBcVec;
625 
626 struct OpJacobian : public UserDataOperator {
627  const int tAg; ///< adol-c tape
628  const bool evalRhs;
629  const bool evalLhs;
630  boost::shared_ptr<DataAtIntegrationPts>
631  dataAtPts; ///< data at integration pts
632  boost::shared_ptr<PhysicalEquations>
633  physicsPtr; ///< material physical equations
634 
635  OpJacobian(const std::string &field_name, const int tag, const bool eval_rhs,
636  const bool eval_lhs,
637  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
638  boost::shared_ptr<PhysicalEquations> &physics_ptr)
639  : UserDataOperator(field_name, OPROW), tAg(tag), evalRhs(eval_rhs),
640  evalLhs(eval_lhs), dataAtPts(data_ptr), physicsPtr(physics_ptr) {}
641 
642  virtual MoFEMErrorCode evaluateRhs(EntData &data) = 0;
643  virtual MoFEMErrorCode evaluateLhs(EntData &data) = 0;
644 
645  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
646 };
647 
648 template <typename T> struct OpAssembleBasic : public T {
649 
650  const bool assembleSymmetry;
651 
652  boost::shared_ptr<DataAtIntegrationPts>
653  dataAtPts; ///< data at integration pts
654 
655  OpAssembleBasic(const std::string &field_name,
656  boost::shared_ptr<DataAtIntegrationPts> data_ptr,
657  const char type)
658  : T(field_name, type), dataAtPts(data_ptr), assembleSymmetry(false) {}
659 
660  OpAssembleBasic(const std::string &row_field, const std::string &col_field,
661  boost::shared_ptr<DataAtIntegrationPts> data_ptr,
662  const char type, const bool assemble_symmetry)
663  : T(row_field, col_field, type, false), dataAtPts(data_ptr),
664  assembleSymmetry(assemble_symmetry) {}
665 
666  VectorDouble nF; ///< local right hand side vector
667  MatrixDouble K; ///< local tangent matrix
669 
672  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
674  }
675 
676  virtual MoFEMErrorCode integrate(int row_side, EntityType row_type,
677  EntData &data) {
679  CHKERR integrate(data);
681  }
682 
683  virtual MoFEMErrorCode assemble(EntData &data) {
685  double *vec_ptr = &*nF.begin();
686  int nb_dofs = data.getIndices().size();
687  int *ind_ptr = &*data.getIndices().begin();
688  CHKERR VecSetValues(T::getFEMethod()->ts_F, nb_dofs, ind_ptr, vec_ptr,
689  ADD_VALUES);
691  }
692 
693  virtual MoFEMErrorCode assemble(int row_side, EntityType row_type,
694  EntData &data) {
696  CHKERR assemble(data);
698  }
699 
700  virtual MoFEMErrorCode integrate(EntData &row_data, EntData &col_data) {
702  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
704  }
705 
706  virtual MoFEMErrorCode assemble(int row_side, int col_side,
707  EntityType row_type, EntityType col_type,
708  EntData &row_data, EntData &col_data) {
710 
711  auto &bmc = dataAtPts->blockMatContainor;
712  int *row_ind_ptr = &*row_data.getIndices().begin();
713  int *col_ind_ptr = &*col_data.getIndices().begin();
714  int row_nb_dofs = row_data.getIndices().size();
715  int col_nb_dofs = col_data.getIndices().size();
716 
717  auto add_block = [&](auto &row_name, auto &col_name, auto row_side,
718  auto col_side, auto row_type, auto col_type,
719  const auto &m, const auto &row_ind,
720  const auto &col_ind) {
721  auto it = bmc.get<0>().find(boost::make_tuple(
722  row_name, col_name, row_type, col_type, row_side, col_side));
723  if (it != bmc.get<0>().end()) {
724  it->setMat(m);
725  it->setInd(row_ind, col_ind);
726  it->setSetAtElement();
727  } else
729  row_name, col_name, row_type, col_type, row_side, col_side, m,
730  row_ind, col_ind));
731  };
732 
733  add_block(this->rowFieldName, this->colFieldName, row_side, col_side,
734  row_type, col_type, K, row_data.getIndices(),
735  col_data.getIndices());
736  if (assembleSymmetry) {
737  transposeK.resize(col_nb_dofs, row_nb_dofs, false);
738  noalias(transposeK) = trans(K);
739  add_block(this->colFieldName, this->rowFieldName, col_side, row_side,
740  col_type, row_type, transposeK, col_data.getIndices(),
741  row_data.getIndices());
742  }
743 
744  double *mat_ptr = &*K.data().begin();
745  CHKERR MatSetValues(T::getFEMethod()->ts_B, row_nb_dofs, row_ind_ptr,
746  col_nb_dofs, col_ind_ptr, mat_ptr, ADD_VALUES);
747  if (assembleSymmetry) {
748  double *mat_ptr = &*transposeK.data().begin();
749  CHKERR MatSetValues(T::getFEMethod()->ts_B, col_nb_dofs, col_ind_ptr,
750  row_nb_dofs, row_ind_ptr, mat_ptr, ADD_VALUES);
751  }
752 
754  }
755 
756  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
758  if (data.getIndices().empty())
760  nF.resize(data.getIndices().size(), false);
761  nF.clear();
762  CHKERR integrate(side, type, data);
763  CHKERR assemble(side, type, data);
765  }
766 
767  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
768  EntityType col_type, EntData &row_data,
769  EntData &col_data) {
771  if (row_data.getIndices().empty())
773  if (col_data.getIndices().empty())
775  K.resize(row_data.getIndices().size(), col_data.getIndices().size(), false);
776  K.clear();
777  CHKERR integrate(row_data, col_data);
778  CHKERR assemble(row_side, col_side, row_type, col_type, row_data, col_data);
780  }
781 };
782 
783 struct OpAssembleVolume : public OpAssembleBasic<VolUserDataOperator> {
784  OpAssembleVolume(const std::string &field,
785  boost::shared_ptr<DataAtIntegrationPts> data_ptr,
786  const char type)
787  : OpAssembleBasic<VolUserDataOperator>(field, data_ptr, type) {}
788 
789  OpAssembleVolume(const std::string &row_field, const std::string &col_field,
790  boost::shared_ptr<DataAtIntegrationPts> data_ptr,
791  const char type, const bool assemble_symmetry)
793  row_field, col_field, data_ptr, type, assemble_symmetry) {}
794 };
795 
796 struct OpAssembleFace : public OpAssembleBasic<FaceUserDataOperator> {
797  OpAssembleFace(const std::string &field,
798  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
799  const char type)
800  : OpAssembleBasic<FaceUserDataOperator>(field, data_ptr, type) {}
801 
802  OpAssembleFace(const std::string &row_field, const std::string &col_field,
803  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
804  const char type, const bool assemble_symmetry)
806  row_field, col_field, data_ptr, type, assemble_symmetry) {}
807 };
808 
810  boost::shared_ptr<DataAtIntegrationPts>
811  dataAtPts; ///< data at integration pts
813  const std::string &field_name,
814  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
815  : VolUserDataOperator(field_name, OPROW), dataAtPts(data_ptr) {}
816  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
817 };
818 
820  OpSpatialEquilibrium(const std::string &field_name,
821  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
822  : OpAssembleVolume(field_name, data_ptr, OPROW) {}
824 };
825 
827  OpSpatialRotation(const std::string &field_name,
828  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
829  : OpAssembleVolume(field_name, data_ptr, OPROW) {}
831 };
832 
834  const double alpha_u;
835  OpSpatialPhysical(const std::string &field_name,
836  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
837  const double alpha)
838  : OpAssembleVolume(field_name, data_ptr, OPROW), alpha_u(alpha) {}
840 };
841 
843  OpSpatialConsistencyP(const std::string &field_name,
844  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
845  : OpAssembleVolume(field_name, data_ptr, OPROW) {}
847 };
848 
850  OpSpatialConsistencyBubble(const std::string &field_name,
851  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
852  : OpAssembleVolume(field_name, data_ptr, OPROW) {}
854 };
855 
857  OpSpatialConsistencyDivTerm(const std::string &field_name,
858  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
859  : OpAssembleVolume(field_name, data_ptr, OPROW) {}
861 };
862 
863 struct OpDispBc : public OpAssembleFace {
864  boost::shared_ptr<BcDispVec> bcDispPtr;
865  OpDispBc(const std::string &field_name,
866  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
867  boost::shared_ptr<BcDispVec> &bc_disp_ptr)
868  : OpAssembleFace(field_name, data_ptr, OPROW), bcDispPtr(bc_disp_ptr) {}
870 };
871 
872 struct OpDispBc_dx : public OpAssembleFace {
873  boost::shared_ptr<BcDispVec> bcDispPtr;
874  OpDispBc_dx(const std::string &row_field_name,
875  const std::string &col_field_name,
876  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
877  boost::shared_ptr<BcDispVec> &bc_disp_ptr)
878  : OpAssembleFace(row_field_name, col_field_name, data_ptr, OPROWCOL,
879  false),
880  bcDispPtr(bc_disp_ptr) {}
881  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
882 };
883 
884 struct OpRotationBc : public OpAssembleFace {
885  boost::shared_ptr<BcRotVec> bcRotPtr;
886  OpRotationBc(const std::string &field_name,
887  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
888  boost::shared_ptr<BcRotVec> &bc_rot_ptr)
889  : OpAssembleFace(field_name, data_ptr, OPROW), bcRotPtr(bc_rot_ptr) {}
891 };
892 
894  boost::shared_ptr<BcRotVec> bcRotPtr;
895  OpRotationBc_dx(const std::string &row_field_name,
896  const std::string &col_field_name,
897  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
898  boost::shared_ptr<BcRotVec> &bc_rot_ptr)
899  : OpAssembleFace(row_field_name, col_field_name, data_ptr, OPROWCOL,
900  false),
901  bcRotPtr(bc_rot_ptr) {}
902  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
903 };
904 
905 struct FeTractionBc;
906 
908  OpTractionBc(std::string row_field, FeTractionBc &bc_fe)
909  : FaceUserDataOperator(row_field, FaceUserDataOperator::OPROW),
910  bcFe(bc_fe) {}
911  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
912 
913 protected:
917 };
918 
919 struct FeTractionBc : public FEMethod {
921  FeTractionBc(MoFEM::Interface &m_field, const std::string field_name,
922  boost::shared_ptr<TractionBcVec> &bc)
923  : FEMethod(), mField(m_field), bcData(bc), fieldName(field_name) {}
925  MoFEMErrorCode postProcess() { return 0; }
926 
927  friend MoFEMErrorCode OpTractionBc::doWork(int side, EntityType type,
928  EntData &data);
929 
930 protected:
932  boost::shared_ptr<TractionBcVec> bcData;
933  std::string fieldName;
934 };
935 
936 template <>
938  EpElement(MoFEM::Interface &m_field, const std::string field_name,
939  boost::shared_ptr<TractionBcVec> &bc)
940  : FeTractionBc(m_field, field_name, bc), EpElementBase() {}
945  }
946 };
947 
949  OpSpatialEquilibrium_dw_dP(const std::string &row_field,
950  const std::string &col_field,
951  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
952  const bool assemble_off = false)
953  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
954  assemble_off) {
955  sYmm = false;
956  }
957  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
958 
959 };
960 
962  const double alpha_u;
963  OpSpatialPhysical_du_du(const std::string &row_field,
964  const std::string &col_field,
965  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
966  const double alpha)
967  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false),
968  alpha_u(alpha) {
969  sYmm = false;
970  }
971  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
972 };
973 
975  OpSpatialPhysical_du_dP(const std::string &row_field,
976  const std::string &col_field,
977  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
978  const bool assemble_off = false)
979  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
980  assemble_off) {
981  sYmm = false;
982  }
983 
984  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
985 };
986 
989  const std::string &row_field, const std::string &col_field,
990  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
991  const bool assemble_off = false)
992  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
993  assemble_off) {
994  sYmm = false;
995  }
996 
997  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
998 
999 };
1000 
1002  OpSpatialPhysical_du_domega(const std::string &row_field,
1003  const std::string &col_field,
1004  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1005  const bool assemble_off = false)
1006  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
1007  assemble_off) {
1008  sYmm = false;
1009  }
1010  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1011 };
1012 
1014  OpSpatialPhysical_du_dx(const std::string &row_field,
1015  const std::string &col_field,
1016  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1017  const bool assemble_off = false)
1018  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
1019  assemble_off) {
1020  sYmm = false;
1021  }
1022  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1023 };
1024 
1026  OpSpatialRotation_domega_dP(const std::string &row_field,
1027  const std::string &col_field,
1028  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1029  const bool assemble_off = false)
1030  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
1031  assemble_off) {
1032  sYmm = false;
1033  }
1034  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1035 
1036 };
1037 
1040  const std::string &row_field, const std::string &col_field,
1041  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1042  const bool assemble_off = false)
1043  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
1044  assemble_off) {
1045  sYmm = false;
1046  }
1047  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1048 
1049 };
1050 
1053  const std::string &row_field, const std::string &col_field,
1054  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1055  const bool assemble_off = false)
1056  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
1057  assemble_off) {
1058  sYmm = false;
1059  }
1060  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1061 
1062 };
1063 
1065  OpSpatialRotation_domega_du(const std::string &row_field,
1066  const std::string &col_field,
1067  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1068  const bool assemble_off = false)
1069  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
1070  assemble_off) {
1071  sYmm = false;
1072  }
1073  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1074 };
1075 
1077  OpSpatialRotation_domega_dx(const std::string &row_field,
1078  const std::string &col_field,
1079  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1080  const bool assemble_off = false)
1081  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
1082  assemble_off) {
1083  sYmm = false;
1084  }
1085  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1086 };
1087 
1090  const std::string &row_field, const std::string &col_field,
1091  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
1092  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false) {
1093  sYmm = false;
1094  }
1095  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1096 };
1097 
1100  const std::string &row_field, const std::string &col_field,
1101  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1102  const bool assemble_off = false)
1103  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
1104  assemble_off) {
1105  sYmm = false;
1106  }
1107  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1108 
1109 
1110 };
1111 
1113  OpSpatialConsistency_dP_dx(const std::string &row_field,
1114  const std::string &col_field,
1115  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1116  const bool assemble_off = false)
1117  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
1118  assemble_off) {
1119  sYmm = false;
1120  }
1121  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1122 };
1123 
1126  const std::string &row_field, const std::string &col_field,
1127  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1128  const bool assemble_off = false)
1129  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
1130  assemble_off) {
1131  sYmm = false;
1132  }
1133  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1134 };
1135 
1137 
1138  moab::Interface &postProcMesh;
1139  std::vector<EntityHandle> &mapGaussPts;
1140  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
1141 
1142  OpPostProcDataStructure(moab::Interface &post_proc_mesh,
1143  std::vector<EntityHandle> &map_gauss_pts,
1144  const std::string field_name,
1145  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
1146  : VolUserDataOperator(field_name, UserDataOperator::OPROW),
1147  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
1148  dataAtPts(data_ptr) {}
1149 
1150  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
1151 };
1152 
1154  OpSpatialPrj(const std::string &row_field,
1155  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
1156  : OpAssembleVolume(row_field, data_ptr, OPROW) {}
1157  MoFEMErrorCode integrate(EntData &row_data);
1158 };
1159 
1161  OpSpatialPrj_dx_dx(const std::string &row_field, const std::string &col_field,
1162  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
1163  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false) {
1164  // FIXME: That is symmetric
1165  sYmm = false;
1166  }
1167  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1168 };
1169 
1171  OpSpatialPrj_dx_dw(const std::string &row_field, const std::string &col_field,
1172  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
1173  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false) {
1174  sYmm = false;
1175  }
1176  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1177 };
1178 
1180  OpSpatialSchurBegin(const std::string &row_field,
1181  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
1182  : OpAssembleVolume(row_field, data_ptr, OPROW) {
1183  sYmm = false;
1184  }
1185  MoFEMErrorCode integrate(EntData &row_data) { return 0; }
1186  MoFEMErrorCode assemble(int row_side, EntityType row_type, EntData &data);
1187 
1188 };
1189 
1191  OpSpatialPreconditionMass(const std::string &row_field, MatrixPtr m_ptr)
1192  : OpAssembleVolume(row_field, nullptr, OPROW), mPtr(m_ptr) {
1193  sYmm = false;
1194  }
1195  MoFEMErrorCode integrate(EntData &row_data);
1196  MoFEMErrorCode assemble(int row_side, EntityType row_type, EntData &data) {
1197  return 0;
1198  }
1199 
1200 private:
1202 
1203 };
1204 
1206  OpSpatialSchurEnd(const std::string &row_field,
1207  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1208  const double eps)
1209  : OpAssembleVolume(row_field, data_ptr, OPROW), eps(eps) {
1210  sYmm = false;
1211 
1212 
1213 
1214  }
1215 
1216 
1217  MoFEMErrorCode integrate(EntData &row_data) { return 0; }
1218  MoFEMErrorCode assemble(int side, EntityType type, EntData &data);
1219  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1220  return assemble(side, type, data);
1221  }
1222 
1223 private:
1224  const double eps;
1225 
1229 
1230  map<std::string, MatrixDouble> invBlockMat;
1231  map<std::string, DataAtIntegrationPts::BlockMatContainor> blockMat;
1232 
1233 };
1234 
1236 
1237  OpCalculateStrainEnergy(const std::string field_name,
1238  boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1239  boost::shared_ptr<double> &e)
1240  : VolUserDataOperator(field_name, VolUserDataOperator::OPROW),
1241  dataAtPts(data_at_pts), energy(e) {}
1242 
1243  MoFEMErrorCode doWork(int side, EntityType type,
1245 
1246 private:
1247  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
1248  boost::shared_ptr<double> energy;
1249 };
1250 
1252 
1253  /**
1254  * \brief Getting interface of core database
1255  * @param uuid unique ID of interface
1256  * @param iface returned pointer to interface
1257  * @return error code
1258  */
1259  MoFEMErrorCode query_interface(const MOFEMuuid &uuid,
1260  UnknownInterface **iface) const;
1261 
1263 
1264  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
1265  boost::shared_ptr<PhysicalEquations> physicalEquations;
1266 
1267  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> elasticFeRhs;
1268  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> elasticFeLhs;
1269  boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> elasticBcLhs;
1270  boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> elasticBcRhs;
1271  boost::shared_ptr<EpFEMethod> schurAssembly;
1272 
1273  DM dM; ///< Coupled problem all fields
1274  DM dmElastic; ///< Elastic problem
1275  DM dmElasticSchurStreach; ///< Sub problem of dmElastic Schur
1276  DM dmElasticSchurBubble; ///< Sub problem of dmElastic Schur
1277  DM dmElasticSchurSpatialDisp; ///< Sub problem of dmElastic Schur
1278  DM dmElasticSchurOmega; ///< Sub problem of dmElastic Schur
1279  DM dmMaterial; ///< Material problem
1280 
1281  const std::string piolaStress;
1282  const std::string eshelbyStress;
1283  const std::string spatialDisp;
1284  const std::string spatialPosition;
1285  const std::string materialDisp;
1286  const std::string streachTensor;
1287  const std::string rotAxis;
1288  const std::string materialGradient;
1289  const std::string tauField;
1290  const std::string lambdaField;
1291  const std::string bubbleField;
1292 
1293  const std::string elementVolumeName;
1294  const std::string naturalBcElement;
1295  const std::string essentialBcElement;
1296 
1297  EshelbianCore(MoFEM::Interface &m_field);
1298  virtual ~EshelbianCore();
1299 
1301  PetscBool ghostFrameOn;
1304  double alpha_u;
1306 
1308 
1309  boost::shared_ptr<BcDispVec> bcSpatialDispVecPtr;
1310  boost::shared_ptr<BcRotVec> bcSpatialRotationVecPtr;
1311  boost::shared_ptr<TractionBcVec> bcSpatialTraction;
1312  boost::shared_ptr<TractionFreeBc> bcSpatialFreeTraction;
1313 
1314  template <typename BC>
1315  MoFEMErrorCode getBc(boost::shared_ptr<BC> &bc_vec_ptr,
1316  const std::string block_set_name,
1317  const int nb_attributes) {
1320  auto block_name = it->getName();
1321  if (block_name.compare(0, block_set_name.length(), block_set_name) == 0) {
1322  std::vector<double> block_attributes;
1323  CHKERR it->getAttributes(block_attributes);
1324  if (block_attributes.size() != nb_attributes) {
1325  SETERRQ3(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1326  "In block %s expected %d attributes, but given %d",
1327  it->getName().c_str(), nb_attributes,
1328  block_attributes.size());
1329  }
1330  Range faces;
1331  CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
1332  true);
1333  bc_vec_ptr->emplace_back(block_name, block_attributes, faces);
1334  }
1335  }
1337  }
1338 
1340  bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
1341  return getBc(bcSpatialDispVecPtr, "SPATIAL_DISP_BC", 6);
1342  }
1343 
1345  bcSpatialRotationVecPtr = boost::make_shared<BcRotVec>();
1346  return getBc(bcSpatialRotationVecPtr, "SPATIAL_ROTATION_BC", 4);
1347  }
1348 
1350  bcSpatialTraction = boost::make_shared<TractionBcVec>();
1351  return getBc(bcSpatialTraction, "SPATIAL_TRACTION_BC", 6);
1352  }
1353 
1355  boost::shared_ptr<TractionFreeBc> &bc_ptr,
1356  const std::string disp_block_set_name,
1357  const std::string rot_block_set_name);
1358  inline MoFEMErrorCode
1361  boost::shared_ptr<TractionFreeBc>(new TractionFreeBc());
1362  return getTractionFreeBc(meshset, bcSpatialFreeTraction, "SPATIAL_DISP_BC",
1363  "SPATIAL_ROTATION_BC");
1364  }
1365 
1366  MoFEMErrorCode addFields(const EntityHandle meshset = 0);
1369  MoFEMErrorCode addDMs(const BitRefLevel bit = BitRefLevel().set(0));
1370 
1372  const double lambda,
1373  const double mu,
1374  const double sigma_y);
1375 
1376  MoFEMErrorCode addMaterial_HMHMooneyRivlin(const int tape, const double alpha,
1377  const double beta,
1378  const double lambda,
1379  const double sigma_y);
1380 
1382  const int tag, const bool do_rhs, const bool do_lhs,
1383  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe);
1384 
1386  const int tag, const bool add_elastic, const bool add_material,
1387  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe_rhs,
1388  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe_lhs);
1389 
1391  const bool add_elastic, const bool add_material,
1392  boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> &fe_rhs,
1393  boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> &fe_lhs);
1394 
1395  MoFEMErrorCode setElasticElementOps(const int tag);
1397 
1398  MoFEMErrorCode setUpTSElastic(TS ts, Mat m, Vec f);
1399  MoFEMErrorCode solveElastic(TS ts, Vec x);
1400 
1401  MoFEMErrorCode postProcessResults(const int tag, const std::string file);
1402 };
1403 
1404 } // namespace EshelbianPlasticity
1405 
1406 #endif //__ESHELBIAN_PLASTICITY_HPP__
boost::shared_ptr< PhysicalEquations > physicsPtr
material physical equations
MoFEMErrorCode setGenericFaceElementOps(const bool add_elastic, const bool add_material, boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore >> &fe_rhs, boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore >> &fe_lhs)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
OpSpatialPhysical_du_domega(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const bool assemble_off=false)
OpSpatialPhysical_du_dBubble(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const bool assemble_off=false)
map< std::string, DataAtIntegrationPts::BlockMatContainor > blockMat
MoFEMErrorCode setBaseVolumeElementOps(const int tag, const bool do_rhs, const bool do_lhs, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore >> &fe)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
VectorDouble nF
local right hand side vector
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
virtual moab::Interface & get_moab()=0
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
OpJacobian(const std::string &field_name, const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, boost::shared_ptr< PhysicalEquations > &physics_ptr)
MoFEMErrorCode setElasticElementOps(const int tag)
OpSpatialPrj(const std::string &row_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr)
OpDispBc(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, boost::shared_ptr< BcDispVec > &bc_disp_ptr)
MoFEMErrorCode integrate(EntData &data)
boost::shared_ptr< EpFEMethod > schurAssembly
OpSpatialSchurEnd(const std::string &row_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const double eps)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
FTensor::Tensor3< adouble, 3, 3, 3 > ATensor3
OpRotationBc(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, boost::shared_ptr< BcRotVec > &bc_rot_ptr)
boost::shared_ptr< BcRotVec > bcRotPtr
multi_index_container< BlockMatData, indexed_by< ordered_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField >, member< BlockMatData, EntityType, &BlockMatData::rowType >, member< BlockMatData, EntityType, &BlockMatData::colType >, member< BlockMatData, int, &BlockMatData::rowSide >, member< BlockMatData, int, &BlockMatData::colSide > > >, ordered_non_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField >, member< BlockMatData, EntityType, &BlockMatData::rowType >, member< BlockMatData, EntityType, &BlockMatData::colType > > >, ordered_non_unique< composite_key< BlockMatData, member< BlockMatData, std::string, &BlockMatData::rowField >, member< BlockMatData, std::string, &BlockMatData::colField > > >, ordered_non_unique< member< BlockMatData, std::string, &BlockMatData::rowField > >, ordered_non_unique< member< BlockMatData, std::string, &BlockMatData::colField > > > > BlockMatContainor
boost::shared_ptr< TractionBcVec > bcData
OpSpatialConsistencyDivTerm(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > &data_ptr)
boost::shared_ptr< BcRotVec > bcSpatialRotationVecPtr
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
PhysicalEquations(const int size_active, const int size_dependent)
OpSpatialRotation(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > &data_ptr)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
OpDispBc_dx(const std::string &row_field_name, const std::string &col_field_name, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, boost::shared_ptr< BcDispVec > &bc_disp_ptr)
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
virtual MoFEMErrorCode assemble(int row_side, EntityType row_type, EntData &data)
std::vector< EntityHandle > & mapGaussPts
OpSpatialRotation_domega_dx(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const bool assemble_off=false)
OpSpatialPhysical_du_du(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const double alpha)
MoFEMErrorCode addDMs(const BitRefLevel bit=BitRefLevel().set(0))
OpSpatialPhysical_du_dP(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const bool assemble_off=false)
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:77
base class for all interface classes
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
VectorBoundedArray< int, 3 > VectorInt3
Definition: Types.hpp:83
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
OpSpatialPhysical_du_dx(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const bool assemble_off=false)
MoFEMErrorCode assemble(int side, EntityType type, EntData &data)
std::vector< BcRot > BcRotVec
MoFEMErrorCode integrate(EntData &data)
virtual MoFEMErrorCode integrate(int row_side, EntityType row_type, EntData &data)
MoFEMErrorCode addSpatialDispStressSchurMatrix(Mat Sw, AO aoSw)
virtual MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
BlockMatData(const std::string row_field, const std::string col_field, EntityType row_type, EntityType col_type, int row_side, int col_side, const MatrixDouble &m, const VectorInt row_ind, VectorInt col_ind)
EpElement(MoFEM::Interface &m_field)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
MoFEMErrorCode addBubbleSchurMatrix(Mat SBubble, AO aoSBubble)
FTensor::Tensor3< double, 3, 3, 3 > getDiffRotationFormVector(FTensor::Tensor1< T, 3 > &t_omega)
virtual MoFEMErrorCode assemble(EntData &data)
DM dmElasticSchurOmega
Sub problem of dmElastic Schur.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpSpatialRotation_domega_dP(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const bool assemble_off=false)
EshelbianCore(MoFEM::Interface &m_field)
std::bitset< BITINTERFACEUID_SIZE > BitIntefaceId
Definition: Types.hpp:56
MoFEMErrorCode integrate(EntData &data)
OpSpatialSchurBegin(const std::string &row_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr)
MoFEMErrorCode addFields(const EntityHandle meshset=0)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode assemble(int row_side, EntityType row_type, EntData &data)
static const MOFEMuuid IDD_MOFEMEshelbianCrackInterface
OpSpatialConsistency_dBubble_domega(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const bool assemble_off=false)
boost::shared_ptr< BcRotVec > bcRotPtr
std::vector< Range > TractionFreeBc
OpCalculateStrainEnergy(const std::string field_name, boost::shared_ptr< DataAtIntegrationPts > data_at_pts, boost::shared_ptr< double > &e)
MoFEMErrorCode integrate(EntData &row_data)
OpAssembleBasic(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const char type, const bool assemble_symmetry)
boost::shared_ptr< VectorDouble > VectorPtr
std::vector< TractionBc > TractionBcVec
OpTractionBc(std::string row_field, FeTractionBc &bc_fe)
OpSpatialConsistency_dBubble_dx(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const bool assemble_off=false)
static DTensor2Ptr get_VecTensor2(std::vector< double > &v)
OpSpatialRotation_domega_du(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const bool assemble_off=false)
boost::shared_ptr< MatrixDouble > MatrixPtr
OpSpatialConsistency_dP_dx(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const bool assemble_off=false)
OpSpatialRotation_domega_domega(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const bool assemble_off=false)
MoFEMErrorCode addStreachSchurMatrix(Mat Suu, AO aoSuu)
virtual MoFEMErrorCode evaluateLhs(EntData &data)=0
MatrixDouble K
local tangent matrix
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
FTensor::Tensor3< FTensor::PackPtr< double *, 1 >, 3, 3, 3 > getFTensor3FromMat(MatrixDouble &m)
OpSpatialConsistencyP(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > &data_ptr)
boost::shared_ptr< BcDispVec > bcSpatialDispVecPtr
double tol
virtual MoFEMErrorCode evaluateRhs(EntData &data)=0
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
OpAssembleBasic(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const char type)
MoFEMErrorCode postProcessResults(const int tag, const std::string file)
MoFEMErrorCode assemble(int row_side, EntityType row_type, EntData &data)
MoFEMErrorCode addBoundaryFiniteElement(const EntityHandle meshset=0)
MoFEMErrorCode addVolumeFiniteElement(const EntityHandle meshset=0)
OpCalculateRotationAndSpatialGradient(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > &data_ptr)
MoFEMErrorCode setGenericVolumeElementOps(const int tag, const bool add_elastic, const bool add_material, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore >> &fe_rhs, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore >> &fe_lhs)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
boost::shared_ptr< TractionFreeBc > bcSpatialFreeTraction
OpSpatialPhysical(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const double alpha)
TractionBc(std::string name, std::vector< double > &attr, Range &faces)
MoFEMErrorCode setUpTSElastic(TS ts, Mat m, Vec f)
void setInd(const VectorInt &row_ind, const VectorInt &col_ind) const
OpSpatialPrj_dx_dx(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
BcRot(std::string name, std::vector< double > &attr, Range &faces)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode addMaterial_HMHHStVenantKirchhoff(const int tape, const double lambda, const double mu, const double sigma_y)
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
MoFEMErrorCode getTractionFreeBc(const EntityHandle meshset, boost::shared_ptr< TractionFreeBc > &bc_ptr, const std::string disp_block_set_name, const std::string rot_block_set_name)
OpSpatialPrj_dx_dw(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr)
static DTensor0Ptr get_VecTensor0(std::vector< double > &v)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
FeTractionBc(MoFEM::Interface &m_field, const std::string field_name, boost::shared_ptr< TractionBcVec > &bc)
boost::shared_ptr< PhysicalEquations > physicalEquations
OpSpatialRotation_domega_dBubble(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const bool assemble_off=false)
MoFEMErrorCode solveElastic(TS ts, Vec x)
FTensor::Tensor2< adouble, 3, 3 > ATensor2
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > DTensor0Ptr
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
boost::shared_ptr< BcDispVec > bcDispPtr
DataForcesAndSourcesCore::EntData EntData
OpPostProcDataStructure(moab::Interface &post_proc_mesh, std::vector< EntityHandle > &map_gauss_pts, const std::string field_name, boost::shared_ptr< DataAtIntegrationPts > &data_ptr)
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
MoFEMErrorCode getBc(boost::shared_ptr< BC > &bc_vec_ptr, const std::string block_set_name, const int nb_attributes)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > elasticBcRhs
#define CHKERR
Inline error check.
Definition: definitions.h:596
MoFEMErrorCode getSpatialTractionFreeBc(const EntityHandle meshset=0)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
OpAssembleVolume(const std::string &field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const char type)
DM dM
Coupled problem all fields.
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
map< std::string, MatrixDouble > invBlockMat
FaceElementForcesAndSourcesCore::UserDataOperator FaceUserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
DM dmElasticSchurStreach
Sub problem of dmElastic Schur.
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:75
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:88
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
virtual MoFEMErrorCode assemble(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
EpElement(MoFEM::Interface &m_field, const std::string field_name, boost::shared_ptr< TractionBcVec > &bc)
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > elasticFeLhs
DM dmElasticSchurSpatialDisp
Sub problem of dmElastic Schur.
OpSpatialConsistency_dP_domega(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > elasticBcLhs
FTensor::Tensor3< FTensor::PackPtr< double *, 1 >, 3, 3, 3 > DTensor3Ptr
const VectorInt & getIndices() const
Get global indices of dofs on entity.
boost::shared_ptr< TractionBcVec > bcSpatialTraction
virtual OpJacobian * returnOpJacobian(const std::string &field_name, const int tag, const bool eval_rhs, const bool eval_lhs, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, boost::shared_ptr< PhysicalEquations > &physics_ptr)=0
MoFEMErrorCode integrate(EntData &row_data)
Data on single entity (This is passed as argument to DataOperator::doWork)
OpAssembleFace(const std::string &field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const char type)
OpSpatialConsistencyBubble(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > &data_ptr)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
OpSpatialEquilibrium_dw_dP(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const bool assemble_off=false)
BcDisp(std::string name, std::vector< double > &attr, Range &faces)
boost::shared_ptr< BcDispVec > bcDispPtr
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
OpRotationBc_dx(const std::string &row_field_name, const std::string &col_field_name, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, boost::shared_ptr< BcRotVec > &bc_rot_ptr)
MoFEMErrorCode integrate(EntData &data)
OpAssembleVolume(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > data_ptr, const char type, const bool assemble_symmetry)
OpSpatialPreconditionMass(const std::string &row_field, MatrixPtr m_ptr)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
static DTensor3Ptr get_vecTensor3(std::vector< double > &v)
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
VolumeElementForcesAndSourcesCore::UserDataOperator VolUserDataOperator
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:76
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
DM dmElasticSchurBubble
Sub problem of dmElastic Schur.
virtual MoFEMErrorCode recordTape(const int tag)=0
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
MoFEMErrorCode addOmegaSchurMatrix(Mat SOmega, AO aoSOmega)
moab::Interface & postProcMesh
MoFEMErrorCode integrate(EntData &row_data)
std::vector< BcDisp > BcDispVec
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
FTensor::Tensor2< double, 3, 3 > getRotationFormVector(FTensor::Tensor1< T, 3 > &t_omega)
FTensor::Tensor1< adouble, 3 > ATensor1
OpAssembleFace(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const char type, const bool assemble_symmetry)
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > elasticFeRhs
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, UnknownInterface **iface) const
Getting interface of core database.
OpSpatialEquilibrium(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > &data_ptr)
virtual MoFEMErrorCode integrate(EntData &data)
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, 3, 3 > DTensor2Ptr
MoFEMErrorCode addMaterial_HMHMooneyRivlin(const int tape, const double alpha, const double beta, const double lambda, const double sigma_y)
constexpr std::enable_if<(Dim0<=2 &&Dim1<=2), Tensor2_Expr< Levi_Civita< T >, T, Dim0, Dim1, i, j > >::type levi_civita(const Index< i, Dim0 > &, const Index< j, Dim1 > &)
levi_civita functions to make for easy adhoc use