v0.9.0
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 
38 
39 struct EpElementBase {
40  Mat Suu;
41  AO aoSuu;
42  Mat SBubble;
44  Mat SOmega;
46  Mat Sw;
47  AO aoSw;
49  : Suu(PETSC_NULL), aoSuu(PETSC_NULL), SBubble(PETSC_NULL),
50  aoSBubble(PETSC_NULL), Sw(PETSC_NULL), aoSw(PETSC_NULL),
51  SOmega(PETSC_NULL), aoSOmega(PETSC_NULL) {}
53  if (Suu)
54  MatDestroy(&Suu);
55  if (aoSuu)
56  AODestroy(&aoSuu);
57  if (SBubble)
58  MatDestroy(&SBubble);
59  if (aoSBubble)
60  AODestroy(&aoSBubble);
61  if (Sw)
62  MatDestroy(&Sw);
63  if (aoSw)
64  AODestroy(&aoSw);
65  if (SOmega)
66  MatDestroy(&SOmega);
67  if (aoSOmega)
68  AODestroy(&aoSOmega);
69  }
70 
73  if (this->Suu)
74  MatDestroy(&this->Suu);
75  if (this->aoSuu)
76  AODestroy(&this->aoSuu);
77  this->Suu = Suu;
78  this->aoSuu = aoSuu;
79  PetscObjectReference((PetscObject)this->Suu);
80  PetscObjectReference((PetscObject)this->aoSuu);
82  }
83 
86  if (this->SBubble)
87  MatDestroy(&this->SBubble);
88  if (this->aoSBubble)
89  AODestroy(&this->aoSBubble);
90  this->SBubble = SBubble;
91  this->aoSBubble = aoSBubble;
92  PetscObjectReference((PetscObject)this->SBubble);
93  PetscObjectReference((PetscObject)this->aoSBubble);
95  }
96 
99  if (this->Sw)
100  MatDestroy(&this->Sw);
101  if (this->aoSw)
102  AODestroy(&this->aoSw);
103  this->Sw = Sw;
104  this->aoSw = aoSw;
105  PetscObjectReference((PetscObject)this->Sw);
106  PetscObjectReference((PetscObject)this->aoSw);
108  }
109 
112  if (this->SOmega)
113  MatDestroy(&this->SOmega);
114  if (this->aoSOmega)
115  AODestroy(&this->aoSOmega);
116  this->SOmega = SOmega;
117  this->aoSOmega = aoSOmega;
118  PetscObjectReference((PetscObject)this->SOmega);
119  PetscObjectReference((PetscObject)this->aoSOmega);
121  }
122 };
123 
124 template <typename E> struct EpElement : public E, EpElementBase {
125  EpElement(MoFEM::Interface &m_field) : E(m_field), EpElementBase() {}
126 };
127 
128 template <> struct EpElement<BasicMethod> : public FEMethod, EpElementBase {
130 };
131 
132 template <> struct EpElement<FEMethod> : public FEMethod, EpElementBase {
134 };
135 
136 struct EpFEMethod : EpElement<FEMethod> {
140  if (Suu)
141  CHKERR MatZeroEntries(Suu);
142  if (SBubble)
143  CHKERR MatZeroEntries(SBubble);
144  if (Sw)
145  CHKERR MatZeroEntries(Sw);
146  if (SOmega)
147  CHKERR MatZeroEntries(SOmega);
149  }
150 
153  auto assemble = [](Mat &a) {
155  if (a) {
156  CHKERR MatAssemblyBegin(a, MAT_FINAL_ASSEMBLY);
157  CHKERR MatAssemblyEnd(a, MAT_FINAL_ASSEMBLY);
158  }
160  };
161  CHKERR assemble(Suu);
162  CHKERR assemble(SBubble);
163  CHKERR assemble(SOmega);
164  CHKERR assemble(Sw);
165  // std::string wait;
166  // CHKERR MatView(SOmega, PETSC_VIEWER_DRAW_WORLD);
167  // std::cin >> wait;
168  // CHKERR MatView(Sw, PETSC_VIEWER_DRAW_WORLD);
169  // std::cin >> wait;
170  // CHKERR MatView(SBubble, PETSC_VIEWER_DRAW_WORLD);
171  // std::cin >> wait;
172  // CHKERR MatView(Suu, PETSC_VIEWER_DRAW_WORLD);
173  // std::cin >> wait;
174  // CHKERR MatView(ts_B, PETSC_VIEWER_DRAW_WORLD);
175  // std::cin >> wait;
177  }
178 };
179 
180 struct PhysicalEquations;
182 
200 
208 
229 
230  std::array<MatrixDouble, 6> expLogUHessian;
231  std::array<MatrixDouble, 3> rotationHessian;
232 
233  // Not really data at integration points, used to calculate Schur complement
234 
235  struct BlockMatData {
236 
237  std::string rowField;
238  std::string colField;
239  EntityType rowType;
240  EntityType colType;
241  int rowSide;
242  int colSide;
246 
248 
249  BlockMatData(const std::string row_field, const std::string col_field,
250  EntityType row_type, EntityType col_type, int row_side,
251  int col_side, const MatrixDouble &m, const VectorInt row_ind,
252  VectorInt col_ind)
253  : rowField(row_field), colField(col_field), rowType(row_type),
254  colType(col_type), rowSide(row_side), colSide(col_side),
255  setAtElement(true) {
256 
257  M.resize(m.size1(), m.size2(), false);
258  noalias(M) = m;
259  rowInd.resize(row_ind.size(), false);
260  noalias(rowInd) = row_ind;
261  colInd.resize(col_ind.size(), false);
262  noalias(colInd) = col_ind;
263  }
264 
265  void setInd(const VectorInt &row_ind, const VectorInt &col_ind) const {
266  auto &const_row_ind = const_cast<VectorInt &>(rowInd);
267  auto &const_col_ind = const_cast<VectorInt &>(colInd);
268  const_row_ind.resize(row_ind.size(), false);
269  noalias(const_row_ind) = row_ind;
270  const_col_ind.resize(col_ind.size(), false);
271  noalias(const_col_ind) = col_ind;
272  }
273 
274  void setMat(const MatrixDouble &m) const {
275  auto &const_m = const_cast<MatrixDouble &>(M);
276  const_m.resize(m.size1(), m.size2(), false);
277  noalias(const_m) = m;
278  }
279 
280  void addMat(const MatrixDouble &m) const {
281  auto &const_m = const_cast<MatrixDouble &>(M);
282  const_m += m;
283  }
284 
285  void clearMat() const {
286  auto &const_m = const_cast<MatrixDouble &>(M);
287  const_m.clear();
288  }
289 
290  void setSetAtElement() const {
291  bool &set = const_cast<bool &>(setAtElement);
292  set = true;
293  }
294 
295  void unSetAtElement() const {
296  bool &set = const_cast<bool &>(setAtElement);
297  set = false;
298  }
299  };
300 
301  typedef multi_index_container<
302  BlockMatData,
303  indexed_by<
304 
305  ordered_unique<
306 
307  composite_key<
308  BlockMatData,
309 
310  member<BlockMatData, std::string, &BlockMatData::rowField>,
311  member<BlockMatData, std::string, &BlockMatData::colField>,
312  member<BlockMatData, EntityType, &BlockMatData::rowType>,
313  member<BlockMatData, EntityType, &BlockMatData::colType>,
314  member<BlockMatData, int, &BlockMatData::rowSide>,
315  member<BlockMatData, int, &BlockMatData::colSide>
316 
317  >>,
318  ordered_non_unique<
319 
320  composite_key<
321  BlockMatData,
322 
323  member<BlockMatData, std::string, &BlockMatData::rowField>,
324  member<BlockMatData, std::string, &BlockMatData::colField>,
325  member<BlockMatData, EntityType, &BlockMatData::rowType>,
326  member<BlockMatData, EntityType, &BlockMatData::colType>
327 
328  >>,
329  ordered_non_unique<
330 
331  composite_key<
332  BlockMatData,
333 
334  member<BlockMatData, std::string, &BlockMatData::rowField>,
335  member<BlockMatData, std::string, &BlockMatData::colField>
336 
337  >>,
338  ordered_non_unique<
339  member<BlockMatData, std::string, &BlockMatData::rowField>>,
340  ordered_non_unique<
341  member<BlockMatData, std::string, &BlockMatData::colField>>
342 
343  >>
345 
349 
350  boost::shared_ptr<PhysicalEquations> physicsPtr;
351 };
352 
353 struct OpJacobian;
354 
356 
361 
364 
365  PhysicalEquations() = delete;
366  PhysicalEquations(const int size_active, const int size_dependent)
367  : activeVariables(size_active, 0), dependentVariables(size_dependent, 0),
368  dependentVariablesDirevatives(size_dependent * size_active, 0) {}
369  virtual ~PhysicalEquations() {}
370 
371  virtual MoFEMErrorCode recordTape(const int tag, DTensor2Ptr *t_h) = 0;
372 
373  virtual OpJacobian *
374  returnOpJacobian(const std::string &field_name, const int tag,
375  const bool eval_rhs, const bool eval_lhs,
376  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
377  boost::shared_ptr<PhysicalEquations> &physics_ptr) = 0;
378 
379  std::vector<double> activeVariables;
380  std::vector<double> dependentVariables;
381  std::vector<double> dependentVariablesDirevatives;
382 
383  /** \name Active variables */
384 
385  /**@{*/
386 
387  template <int S>
388  inline static DTensor2Ptr get_VecTensor2(std::vector<double> &v) {
389  return DTensor2Ptr(&v[S + 0], &v[S + 1], &v[S + 2], &v[S + 3], &v[S + 4],
390  &v[S + 5], &v[S + 6], &v[S + 7], &v[S + 8]);
391  }
392 
393  template <int S>
394  inline static DTensor0Ptr get_VecTensor0(std::vector<double> &v) {
395  return DTensor0Ptr(&v[S + 0]);
396  }
397 
398  template <int S0, int S1, int S2>
399  inline static DTensor3Ptr get_vecTensor3(std::vector<double> &v) {
400  constexpr int A00 = 18 * S0 + 18 * 0 + 9 * S1 + 3 * S2;
401  constexpr int A01 = 18 * S0 + 18 * 1 + 9 * S1 + 3 * S2;
402  constexpr int A02 = 18 * S0 + 18 * 2 + 9 * S1 + 3 * S2;
403  constexpr int A10 = 18 * S0 + 18 * 3 + 9 * S1 + 3 * S2;
404  constexpr int A11 = 18 * S0 + 18 * 4 + 9 * S1 + 3 * S2;
405  constexpr int A12 = 18 * S0 + 18 * 5 + 9 * S1 + 3 * S2;
406  constexpr int A20 = 18 * S0 + 18 * 6 + 9 * S1 + 3 * S2;
407  constexpr int A21 = 18 * S0 + 18 * 7 + 9 * S1 + 3 * S2;
408  constexpr int A22 = 18 * S0 + 18 * 8 + 9 * S1 + 3 * S2;
409  return DTensor3Ptr(
410 
411  &v[A00 + 0], &v[A00 + 1], &v[A00 + 2], &v[A01 + 0], &v[A01 + 1],
412  &v[A01 + 2], &v[A02 + 0], &v[A02 + 1], &v[A02 + 2],
413 
414  &v[A10 + 0], &v[A10 + 1], &v[A10 + 2], &v[A11 + 0], &v[A11 + 1],
415  &v[A11 + 2], &v[A12 + 0], &v[A12 + 1], &v[A12 + 2],
416 
417  &v[A20 + 0], &v[A20 + 1], &v[A20 + 2], &v[A21 + 0], &v[A21 + 1],
418  &v[A21 + 2], &v[A22 + 0], &v[A22 + 1], &v[A22 + 2]
419 
420  );
421  }
422 
423  /**@}*/
424 
425  /** \name Active variables */
426 
427  /**@{*/
428 
429  inline DTensor2Ptr get_h() { return get_VecTensor2<0>(activeVariables); }
430  inline DTensor2Ptr get_H() { return get_VecTensor2<9>(activeVariables); }
431 
432  /**@}*/
433 
434  /** \name Dependent variables */
435 
436  /**@{*/
437 
438  inline DTensor2Ptr get_P() { return get_VecTensor2<0>(dependentVariables); }
440  return get_VecTensor2<9>(dependentVariables);
441  }
442  inline DTensor0Ptr get_Phi() {
443  return get_VecTensor0<18>(dependentVariables);
444  }
445  inline double &get_PhiRef() { return dependentVariables[18]; }
447  return get_VecTensor2<9 + 9 + 1>(dependentVariables);
448  }
449 
450  /**@}*/
451 
452  /** \name Derivatives of dependent variables */
453 
454  /**@{*/
455 
457  return get_vecTensor3<0, 0, 0>(dependentVariablesDirevatives);
458  }
460  return get_vecTensor3<0, 1, 0>(dependentVariablesDirevatives);
461  }
463  return get_vecTensor3<0, 0, 1>(dependentVariablesDirevatives);
464  }
466  return get_vecTensor3<0, 1, 1>(dependentVariablesDirevatives);
467  }
469  return get_vecTensor3<0, 0, 2>(dependentVariablesDirevatives);
470  }
472  return get_vecTensor3<0, 1, 2>(dependentVariablesDirevatives);
473  }
474 
476  return get_vecTensor3<9, 0, 0>(dependentVariablesDirevatives);
477  }
479  return get_vecTensor3<9, 1, 0>(dependentVariablesDirevatives);
480  }
482  return get_vecTensor3<9, 0, 1>(dependentVariablesDirevatives);
483  }
485  return get_vecTensor3<9, 1, 1>(dependentVariablesDirevatives);
486  }
488  return get_vecTensor3<9, 0, 2>(dependentVariablesDirevatives);
489  }
491  return get_vecTensor3<9, 1, 2>(dependentVariablesDirevatives);
492  }
493 
495  return get_VecTensor2<18 * (9 + 9) + 0>(dependentVariablesDirevatives);
496  }
498  return get_VecTensor2<18 * (9 + 9) + 9>(dependentVariablesDirevatives);
499  }
500 
502  return get_vecTensor3<9 + 9 + 1, 0, 0>(dependentVariablesDirevatives);
503  }
505  return get_vecTensor3<9 + 9 + 1, 1, 0>(dependentVariablesDirevatives);
506  }
508  return get_vecTensor3<9 + 9 + 1, 0, 1>(dependentVariablesDirevatives);
509  }
511  return get_vecTensor3<9 + 9 + 1, 1, 1>(dependentVariablesDirevatives);
512  }
514  return get_vecTensor3<9 + 9 + 1, 0, 2>(dependentVariablesDirevatives);
515  }
517  return get_vecTensor3<9 + 9 + 1, 1, 2>(dependentVariablesDirevatives);
518  }
519 
520  /**@}*/
521 };
522 
523 struct BcDisp {
524  BcDisp(std::string name, std::vector<double> &attr, Range &faces);
525  std::string blockName;
526  Range faces;
529 };
530 typedef std::vector<BcDisp> BcDispVec;
531 
532 struct BcRot {
533  BcRot(std::string name, std::vector<double> &attr, Range &faces);
534  std::string blockName;
535  Range faces;
537  double theta;
538 };
539 typedef std::vector<BcRot> BcRotVec;
540 
541 typedef std::vector<Range> TractionFreeBc;
542 
543 struct TractionBc {
544  TractionBc(std::string name, std::vector<double> &attr, Range &faces);
545  std::string blockName;
546  Range faces;
549 };
550 typedef std::vector<TractionBc> TractionBcVec;
551 
552 struct OpJacobian : public UserDataOperator {
553  const int tAg; ///< adol-c tape
554  const bool evalRhs;
555  const bool evalLhs;
556  boost::shared_ptr<DataAtIntegrationPts>
557  dataAtPts; ///< data at integration pts
558  boost::shared_ptr<PhysicalEquations>
559  physicsPtr; ///< material physical equations
560 
561  OpJacobian(const std::string &field_name, const int tag, const bool eval_rhs,
562  const bool eval_lhs,
563  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
564  boost::shared_ptr<PhysicalEquations> &physics_ptr)
565  : UserDataOperator(field_name, OPROW), tAg(tag), evalRhs(eval_rhs),
566  evalLhs(eval_lhs), dataAtPts(data_ptr), physicsPtr(physics_ptr) {}
567 
568  virtual MoFEMErrorCode evaluateRhs(EntData &data) = 0;
569  virtual MoFEMErrorCode evaluateLhs(EntData &data) = 0;
570 
571  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
572 };
573 
574 template <typename T> struct OpAssembleBasic : public T {
575 
576  const bool assembleSymmetry;
577 
578  boost::shared_ptr<DataAtIntegrationPts>
579  dataAtPts; ///< data at integration pts
580 
581  OpAssembleBasic(const std::string &field_name,
582  boost::shared_ptr<DataAtIntegrationPts> data_ptr,
583  const char type)
584  : T(field_name, type), dataAtPts(data_ptr), assembleSymmetry(false) {}
585 
586  OpAssembleBasic(const std::string &row_field, const std::string &col_field,
587  boost::shared_ptr<DataAtIntegrationPts> data_ptr,
588  const char type, const bool assemble_symmetry)
589  : T(row_field, col_field, type, false), dataAtPts(data_ptr),
590  assembleSymmetry(assemble_symmetry) {}
591 
592  VectorDouble nF; ///< local right hand side vector
593  MatrixDouble K; ///< local tangent matrix
595 
598  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
600  }
601 
602  virtual MoFEMErrorCode integrate(int row_side, EntityType row_type,
603  EntData &data) {
605  CHKERR integrate(data);
607  }
608 
609  virtual MoFEMErrorCode assemble(EntData &data) {
611  double *vec_ptr = &*nF.begin();
612  int nb_dofs = data.getIndices().size();
613  int *ind_ptr = &*data.getIndices().begin();
614  CHKERR VecSetValues(T::getFEMethod()->ts_F, nb_dofs, ind_ptr, vec_ptr,
615  ADD_VALUES);
617  }
618 
619  virtual MoFEMErrorCode assemble(int row_side, EntityType row_type,
620  EntData &data) {
622  CHKERR assemble(data);
624  }
625 
626  virtual MoFEMErrorCode integrate(EntData &row_data, EntData &col_data) {
628  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
630  }
631 
632  virtual MoFEMErrorCode assemble(int row_side, int col_side,
633  EntityType row_type, EntityType col_type,
634  EntData &row_data, EntData &col_data) {
636 
637  auto &bmc = dataAtPts->blockMatContainor;
638  int *row_ind_ptr = &*row_data.getIndices().begin();
639  int *col_ind_ptr = &*col_data.getIndices().begin();
640  int row_nb_dofs = row_data.getIndices().size();
641  int col_nb_dofs = col_data.getIndices().size();
642 
643  auto add_block = [&](auto &row_name, auto &col_name, auto row_side,
644  auto col_side, auto row_type, auto col_type,
645  const auto &m, const auto &row_ind,
646  const auto &col_ind) {
647  auto it = bmc.get<0>().find(boost::make_tuple(
648  row_name, col_name, row_type, col_type, row_side, col_side));
649  if (it != bmc.get<0>().end()) {
650  it->setMat(m);
651  it->setInd(row_ind, col_ind);
652  it->setSetAtElement();
653  } else
655  row_name, col_name, row_type, col_type, row_side, col_side, m,
656  row_ind, col_ind));
657  };
658 
659  add_block(this->rowFieldName, this->colFieldName, row_side, col_side,
660  row_type, col_type, K, row_data.getIndices(),
661  col_data.getIndices());
662  if (assembleSymmetry) {
663  transposeK.resize(col_nb_dofs, row_nb_dofs, false);
664  noalias(transposeK) = trans(K);
665  add_block(this->colFieldName, this->rowFieldName, col_side, row_side,
666  col_type, row_type, transposeK, col_data.getIndices(),
667  row_data.getIndices());
668  }
669 
670  double *mat_ptr = &*K.data().begin();
671  CHKERR MatSetValues(T::getFEMethod()->ts_B, row_nb_dofs, row_ind_ptr,
672  col_nb_dofs, col_ind_ptr, mat_ptr, ADD_VALUES);
673  if (assembleSymmetry) {
674  double *mat_ptr = &*transposeK.data().begin();
675  CHKERR MatSetValues(T::getFEMethod()->ts_B, col_nb_dofs, col_ind_ptr,
676  row_nb_dofs, row_ind_ptr, mat_ptr, ADD_VALUES);
677  }
678 
680  }
681 
682  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
684  if (data.getIndices().empty())
686  nF.resize(data.getIndices().size(), false);
687  nF.clear();
688  CHKERR integrate(side, type, data);
689  CHKERR assemble(side, type, data);
691  }
692 
693  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
694  EntityType col_type, EntData &row_data,
695  EntData &col_data) {
697  if (row_data.getIndices().empty())
699  if (col_data.getIndices().empty())
701  K.resize(row_data.getIndices().size(), col_data.getIndices().size(), false);
702  K.clear();
703  CHKERR integrate(row_data, col_data);
704  CHKERR assemble(row_side, col_side, row_type, col_type, row_data, col_data);
706  }
707 };
708 
709 struct OpAssembleVolume : public OpAssembleBasic<VolUserDataOperator> {
710  OpAssembleVolume(const std::string &field,
711  boost::shared_ptr<DataAtIntegrationPts> data_ptr,
712  const char type)
713  : OpAssembleBasic<VolUserDataOperator>(field, data_ptr, type) {}
714 
715  OpAssembleVolume(const std::string &row_field, const std::string &col_field,
716  boost::shared_ptr<DataAtIntegrationPts> data_ptr,
717  const char type, const bool assemble_symmetry)
718  : OpAssembleBasic<VolUserDataOperator>(row_field, col_field, data_ptr,
719  type, assemble_symmetry) {}
720 };
721 
722 struct OpAssembleFace : public OpAssembleBasic<FaceUserDataOperator> {
723  OpAssembleFace(const std::string &field,
724  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
725  const char type)
726  : OpAssembleBasic<FaceUserDataOperator>(field, data_ptr, type) {}
727 
728  OpAssembleFace(const std::string &row_field, const std::string &col_field,
729  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
730  const char type, const bool assemble_symmetry)
731  : OpAssembleBasic<FaceUserDataOperator>(row_field, col_field, data_ptr,
732  type, assemble_symmetry) {}
733 };
734 
736  boost::shared_ptr<DataAtIntegrationPts>
737  dataAtPts; ///< data at integration pts
739  const std::string &field_name,
740  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
741  : VolUserDataOperator(field_name, OPROW), dataAtPts(data_ptr) {}
742  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
743 };
744 
746  const double alpha_w;
747  OpSpatialEquilibrium(const std::string &field_name,
748  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
749  const double alpha)
750  : OpAssembleVolume(field_name, data_ptr, OPROW), alpha_w(alpha) {}
752 };
753 
755  OpSpatialRotation(const std::string &field_name,
756  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
757  : OpAssembleVolume(field_name, data_ptr, OPROW) {}
759 };
760 
762  const double alpha_u;
763  OpSpatialPhysical(const std::string &field_name,
764  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
765  const double alpha)
766  : OpAssembleVolume(field_name, data_ptr, OPROW), alpha_u(alpha) {}
768 };
769 
771  OpSpatialConsistencyP(const std::string &field_name,
772  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
773  : OpAssembleVolume(field_name, data_ptr, OPROW) {}
775 };
776 
778  OpSpatialConsistencyBubble(const std::string &field_name,
779  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
780  : OpAssembleVolume(field_name, data_ptr, OPROW) {}
782 };
783 
785  OpSpatialConsistencyDivTerm(const std::string &field_name,
786  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
787  : OpAssembleVolume(field_name, data_ptr, OPROW) {}
789 };
790 
791 struct OpDispBc : public OpAssembleFace {
792  boost::shared_ptr<BcDispVec> bcDispPtr;
793  OpDispBc(const std::string &field_name,
794  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
795  boost::shared_ptr<BcDispVec> &bc_disp_ptr)
796  : OpAssembleFace(field_name, data_ptr, OPROW), bcDispPtr(bc_disp_ptr) {}
798 };
799 
800 struct OpDispBc_dx : public OpAssembleFace {
801  boost::shared_ptr<BcDispVec> bcDispPtr;
802  OpDispBc_dx(const std::string &row_field_name,
803  const std::string &col_field_name,
804  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
805  boost::shared_ptr<BcDispVec> &bc_disp_ptr)
806  : OpAssembleFace(row_field_name, col_field_name, data_ptr, OPROWCOL,
807  false),
808  bcDispPtr(bc_disp_ptr) {}
809  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
810 };
811 
812 struct OpRotationBc : public OpAssembleFace {
813  boost::shared_ptr<BcRotVec> bcRotPtr;
814  OpRotationBc(const std::string &field_name,
815  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
816  boost::shared_ptr<BcRotVec> &bc_rot_ptr)
817  : OpAssembleFace(field_name, data_ptr, OPROW), bcRotPtr(bc_rot_ptr) {}
819 };
820 
822  boost::shared_ptr<BcRotVec> bcRotPtr;
823  OpRotationBc_dx(const std::string &row_field_name,
824  const std::string &col_field_name,
825  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
826  boost::shared_ptr<BcRotVec> &bc_rot_ptr)
827  : OpAssembleFace(row_field_name, col_field_name, data_ptr, OPROWCOL,
828  false),
829  bcRotPtr(bc_rot_ptr) {}
830  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
831 };
832 
833 struct FeTractionBc;
834 
836  OpTractionBc(std::string row_field, FeTractionBc &bc_fe)
837  : FaceUserDataOperator(row_field, FaceUserDataOperator::OPROW),
838  bcFe(bc_fe) {}
839  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
840 
841 protected:
845 };
846 
847 struct FeTractionBc : public FEMethod {
849  FeTractionBc(MoFEM::Interface &m_field, const std::string field_name,
850  boost::shared_ptr<TractionBcVec> &bc)
851  : FEMethod(), mField(m_field), bcData(bc), fieldName(field_name) {}
853  MoFEMErrorCode postProcess() { return 0; }
854 
855  friend MoFEMErrorCode OpTractionBc::doWork(int side, EntityType type,
856  EntData &data);
857 
858 protected:
860  boost::shared_ptr<TractionBcVec> bcData;
861  std::string fieldName;
862 };
863 
864 template <>
866  EpElement(MoFEM::Interface &m_field, const std::string field_name,
867  boost::shared_ptr<TractionBcVec> &bc)
868  : FeTractionBc(m_field, field_name, bc), EpElementBase() {}
873  }
874 };
875 
877  OpSpatialEquilibrium_dw_dP(const std::string &row_field,
878  const std::string &col_field,
879  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
880  const bool assemble_off = false)
881  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
882  assemble_off) {
883  sYmm = false;
884  }
885  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
886 };
887 
889  const double alpha_w;
890  OpSpatialEquilibrium_dw_dw(const std::string &row_field,
891  const std::string &col_field,
892  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
893  const double alpha)
894  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
895  false), alpha_w(alpha) {
896  sYmm = true;
897  }
898  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
899 };
900 
902  const double alpha_u;
903  OpSpatialPhysical_du_du(const std::string &row_field,
904  const std::string &col_field,
905  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
906  const double alpha)
907  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false),
908  alpha_u(alpha) {
909  sYmm = false;
910  }
911  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
912 };
913 
915  OpSpatialPhysical_du_dP(const std::string &row_field,
916  const std::string &col_field,
917  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
918  const bool assemble_off = false)
919  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
920  assemble_off) {
921  sYmm = false;
922  }
923 
924  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
925 };
926 
929  const std::string &row_field, const std::string &col_field,
930  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
931  const bool assemble_off = false)
932  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
933  assemble_off) {
934  sYmm = false;
935  }
936 
937  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
938 };
939 
941  OpSpatialPhysical_du_domega(const std::string &row_field,
942  const std::string &col_field,
943  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
944  const bool assemble_off = false)
945  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
946  assemble_off) {
947  sYmm = false;
948  }
949  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
950 };
951 
953  OpSpatialPhysical_du_dx(const std::string &row_field,
954  const std::string &col_field,
955  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
956  const bool assemble_off = false)
957  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
958  assemble_off) {
959  sYmm = false;
960  }
961  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
962 };
963 
965  OpSpatialRotation_domega_dP(const std::string &row_field,
966  const std::string &col_field,
967  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
968  const bool assemble_off = false)
969  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
970  assemble_off) {
971  sYmm = false;
972  }
973  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
974 };
975 
978  const std::string &row_field, const std::string &col_field,
979  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
980  const bool assemble_off = false)
981  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
982  assemble_off) {
983  sYmm = false;
984  }
985  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
986 };
987 
990  const std::string &row_field, const std::string &col_field,
991  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
992  const bool assemble_off = false)
993  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
994  assemble_off) {
995  sYmm = false;
996  }
997  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
998 };
999 
1001  OpSpatialRotation_domega_du(const std::string &row_field,
1002  const std::string &col_field,
1003  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1004  const bool assemble_off = false)
1005  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
1006  assemble_off) {
1007  sYmm = false;
1008  }
1009  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1010 };
1011 
1013  OpSpatialRotation_domega_dx(const std::string &row_field,
1014  const std::string &col_field,
1015  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1016  const bool assemble_off = false)
1017  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
1018  assemble_off) {
1019  sYmm = false;
1020  }
1021  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1022 };
1023 
1026  const std::string &row_field, const std::string &col_field,
1027  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
1028  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false) {
1029  sYmm = false;
1030  }
1031  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1032 };
1033 
1036  const std::string &row_field, const std::string &col_field,
1037  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1038  const bool assemble_off = false)
1039  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
1040  assemble_off) {
1041  sYmm = false;
1042  }
1043  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1044 };
1045 
1047  OpSpatialConsistency_dP_dx(const std::string &row_field,
1048  const std::string &col_field,
1049  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1050  const bool assemble_off = false)
1051  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
1052  assemble_off) {
1053  sYmm = false;
1054  }
1055  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1056 };
1057 
1060  const std::string &row_field, const std::string &col_field,
1061  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1062  const bool assemble_off = false)
1063  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
1064  assemble_off) {
1065  sYmm = false;
1066  }
1067  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1068 };
1069 
1071 
1072  moab::Interface &postProcMesh;
1073  std::vector<EntityHandle> &mapGaussPts;
1074  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
1075 
1076  OpPostProcDataStructure(moab::Interface &post_proc_mesh,
1077  std::vector<EntityHandle> &map_gauss_pts,
1078  const std::string field_name,
1079  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
1080  : VolUserDataOperator(field_name, UserDataOperator::OPROW),
1081  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
1082  dataAtPts(data_ptr) {}
1083 
1084  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
1085 };
1086 
1088  OpSpatialPrj(const std::string &row_field,
1089  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
1090  : OpAssembleVolume(row_field, data_ptr, OPROW) {}
1091  MoFEMErrorCode integrate(EntData &row_data);
1092 };
1093 
1095  OpSpatialPrj_dx_dx(const std::string &row_field, const std::string &col_field,
1096  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
1097  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false) {
1098  // FIXME: That is symmetric
1099  sYmm = false;
1100  }
1101  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1102 };
1103 
1105  OpSpatialPrj_dx_dw(const std::string &row_field, const std::string &col_field,
1106  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
1107  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false) {
1108  sYmm = false;
1109  }
1110  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1111 };
1112 
1114  OpSpatialSchurBegin(const std::string &row_field,
1115  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
1116  : OpAssembleVolume(row_field, data_ptr, OPROW) {
1117  sYmm = false;
1118  }
1119  MoFEMErrorCode integrate(EntData &row_data) { return 0; }
1120  MoFEMErrorCode assemble(int row_side, EntityType row_type, EntData &data);
1121 };
1122 
1124  OpSpatialPreconditionMass(const std::string &row_field, MatrixPtr m_ptr)
1125  : OpAssembleVolume(row_field, nullptr, OPROW), mPtr(m_ptr) {
1126  sYmm = false;
1127  }
1128  MoFEMErrorCode integrate(EntData &row_data);
1129  MoFEMErrorCode assemble(int row_side, EntityType row_type, EntData &data) {
1130  return 0;
1131  }
1132 
1133 private:
1135 };
1136 
1138  OpSpatialSchurEnd(const std::string &row_field,
1139  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1140  const double eps)
1141  : OpAssembleVolume(row_field, data_ptr, OPROW), eps(eps) {
1142  sYmm = false;
1143  }
1144 
1145  MoFEMErrorCode integrate(EntData &row_data) { return 0; }
1146  MoFEMErrorCode assemble(int side, EntityType type, EntData &data);
1147  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1148  return assemble(side, type, data);
1149  }
1150 
1151 private:
1152  const double eps;
1153 
1157 
1158  map<std::string, MatrixDouble> invBlockMat;
1159  map<std::string, DataAtIntegrationPts::BlockMatContainor> blockMat;
1160 };
1161 
1163 
1164  OpCalculateStrainEnergy(const std::string field_name,
1165  boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1166  boost::shared_ptr<double> &e)
1167  : VolUserDataOperator(field_name, VolUserDataOperator::OPROW),
1168  dataAtPts(data_at_pts), energy(e) {}
1169 
1170  MoFEMErrorCode doWork(int side, EntityType type,
1172 
1173 private:
1174  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
1175  boost::shared_ptr<double> energy;
1176 };
1177 
1179 
1180  /**
1181  * \brief Getting interface of core database
1182  * @param uuid unique ID of interface
1183  * @param iface returned pointer to interface
1184  * @return error code
1185  */
1186  MoFEMErrorCode query_interface(const MOFEMuuid &uuid,
1187  UnknownInterface **iface) const;
1188 
1190 
1191  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
1192  boost::shared_ptr<PhysicalEquations> physicalEquations;
1193 
1194  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> elasticFeRhs;
1195  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> elasticFeLhs;
1196  boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> elasticBcLhs;
1197  boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> elasticBcRhs;
1198  boost::shared_ptr<EpFEMethod> schurAssembly;
1199 
1200  SmartPetscObj<DM> dM; ///< Coupled problem all fields
1201  SmartPetscObj<DM> dmElastic; ///< Elastic problem
1202  SmartPetscObj<DM> dmElasticSchurStreach; ///< Sub problem of dmElastic Schur
1203  SmartPetscObj<DM> dmElasticSchurBubble; ///< Sub problem of dmElastic Schur
1204  SmartPetscObj<DM>
1205  dmElasticSchurSpatialDisp; ///< Sub problem of dmElastic Schur
1206  SmartPetscObj<DM> dmElasticSchurOmega; ///< Sub problem of dmElastic Schur
1207  SmartPetscObj<DM> dmMaterial; ///< Material problem
1208 
1209  const std::string piolaStress;
1210  const std::string eshelbyStress;
1211  const std::string spatialDisp;
1212  const std::string materialDisp;
1213  const std::string streachTensor;
1214  const std::string rotAxis;
1215  const std::string materialGradient;
1216  const std::string tauField;
1217  const std::string lambdaField;
1218  const std::string bubbleField;
1219 
1220  const std::string elementVolumeName;
1221  const std::string naturalBcElement;
1222  const std::string essentialBcElement;
1223 
1224  EshelbianCore(MoFEM::Interface &m_field);
1225  virtual ~EshelbianCore();
1226 
1229  double alpha_u;
1230  double alpha_w;
1232 
1234 
1235  boost::shared_ptr<BcDispVec> bcSpatialDispVecPtr;
1236  boost::shared_ptr<BcRotVec> bcSpatialRotationVecPtr;
1237  boost::shared_ptr<TractionBcVec> bcSpatialTraction;
1238  boost::shared_ptr<TractionFreeBc> bcSpatialFreeTraction;
1239 
1240  template <typename BC>
1241  MoFEMErrorCode getBc(boost::shared_ptr<BC> &bc_vec_ptr,
1242  const std::string block_set_name,
1243  const int nb_attributes) {
1246  auto block_name = it->getName();
1247  if (block_name.compare(0, block_set_name.length(), block_set_name) == 0) {
1248  std::vector<double> block_attributes;
1249  CHKERR it->getAttributes(block_attributes);
1250  if (block_attributes.size() != nb_attributes) {
1251  SETERRQ3(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1252  "In block %s expected %d attributes, but given %d",
1253  it->getName().c_str(), nb_attributes,
1254  block_attributes.size());
1255  }
1256  Range faces;
1257  CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
1258  true);
1259  bc_vec_ptr->emplace_back(block_name, block_attributes, faces);
1260  }
1261  }
1263  }
1264 
1266  bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
1267  return getBc(bcSpatialDispVecPtr, "SPATIAL_DISP_BC", 6);
1268  }
1269 
1271  bcSpatialRotationVecPtr = boost::make_shared<BcRotVec>();
1272  return getBc(bcSpatialRotationVecPtr, "SPATIAL_ROTATION_BC", 4);
1273  }
1274 
1276  bcSpatialTraction = boost::make_shared<TractionBcVec>();
1277  return getBc(bcSpatialTraction, "SPATIAL_TRACTION_BC", 6);
1278  }
1279 
1281  boost::shared_ptr<TractionFreeBc> &bc_ptr,
1282  const std::string disp_block_set_name,
1283  const std::string rot_block_set_name);
1284  inline MoFEMErrorCode
1287  boost::shared_ptr<TractionFreeBc>(new TractionFreeBc());
1288  return getTractionFreeBc(meshset, bcSpatialFreeTraction, "SPATIAL_DISP_BC",
1289  "SPATIAL_ROTATION_BC");
1290  }
1291 
1292  MoFEMErrorCode addFields(const EntityHandle meshset = 0);
1295  MoFEMErrorCode addDMs(const BitRefLevel bit = BitRefLevel().set(0));
1296 
1298  const double lambda,
1299  const double mu,
1300  const double sigma_y);
1301 
1302  MoFEMErrorCode addMaterial_HMHMooneyRivlin(const int tape, const double alpha,
1303  const double beta,
1304  const double lambda,
1305  const double sigma_y);
1306 
1308  const int tag, const bool do_rhs, const bool do_lhs,
1309  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe);
1310 
1312  const int tag, const bool add_elastic, const bool add_material,
1313  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe_rhs,
1314  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe_lhs);
1315 
1317  const bool add_elastic, const bool add_material,
1318  boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> &fe_rhs,
1319  boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> &fe_lhs);
1320 
1321  MoFEMErrorCode setElasticElementOps(const int tag);
1323 
1324  MoFEMErrorCode setUpTSElastic(TS ts, Mat m, Vec f);
1325  MoFEMErrorCode solveElastic(TS ts, Vec x);
1326 
1327  MoFEMErrorCode postProcessResults(const int tag, const std::string file);
1328 };
1329 
1330 } // namespace EshelbianPlasticity
1331 
1332 #endif //__ESHELBIAN_PLASTICITY_HPP__
OpSpatialEquilibrium_dw_dw(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const double alpha)
boost::shared_ptr< PhysicalEquations > physicsPtr
material physical equations
SmartPetscObj< DM > dmElastic
Elastic problem.
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)
OpSpatialEquilibrium(const std::string &field_name, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const double alpha)
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
SmartPetscObj< DM > dmElasticSchurOmega
Sub problem of dmElastic Schur.
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:505
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:74
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:481
VectorBoundedArray< int, 3 > VectorInt3
Definition: Types.hpp:82
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
data at integration pts
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
SmartPetscObj< DM > dmMaterial
Material problem.
OpSpatialPhysical_du_dx(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const bool assemble_off=false)
std::array< MatrixDouble, 6 > expLogUHessian
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:512
MoFEMErrorCode addBubbleSchurMatrix(Mat SBubble, AO aoSBubble)
virtual MoFEMErrorCode assemble(EntData &data)
std::array< MatrixDouble, 3 > rotationHessian
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
const VectorInt & getIndices() const
Get global indices of dofs on entity.
std::vector< TractionBc > TractionBcVec
OpTractionBc(std::string row_field, FeTractionBc &bc_fe)
virtual MoFEMErrorCode recordTape(const int tag, DTensor2Ptr *t_h)=0
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
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)
SmartPetscObj< DM > dmElasticSchurStreach
Sub problem of dmElastic Schur.
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)
SmartPetscObj< DM > dmElasticSchurSpatialDisp
Sub problem of dmElastic Schur.
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:600
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)
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)
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:71
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:87
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)
SmartPetscObj< DM > dM
Coupled problem all fields.
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > elasticFeLhs
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
SmartPetscObj< DM > dmElasticSchurBubble
Sub problem of dmElastic Schur.
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:411
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:72
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
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)
boost::shared_ptr< PhysicalEquations > physicsPtr
moab::Interface & postProcMesh
MoFEMErrorCode integrate(EntData &row_data)
std::vector< BcDisp > BcDispVec
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
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.
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)