v0.10.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(getTSB(), PETSC_VIEWER_DRAW_WORLD);
175  // std::cin >> wait;
177  }
178 };
179 
180 struct PhysicalEquations;
182  : public boost::enable_shared_from_this<DataAtIntegrationPts> {
183 
192 
203 
211 
232 
234  return boost::shared_ptr<MatrixDouble>(shared_from_this(),
236  }
238  return boost::shared_ptr<MatrixDouble>(shared_from_this(), &approxPAtPts);
239  }
240 
242  return boost::shared_ptr<MatrixDouble>(shared_from_this(), &divPAtPts);
243  }
244 
246  return boost::shared_ptr<MatrixDouble>(shared_from_this(), &divSigmaAtPts);
247  }
248 
250  return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wAtPts);
251  }
252 
254  return boost::shared_ptr<MatrixDouble>(shared_from_this(), &wDotAtPts);
255  }
256 
258  return boost::shared_ptr<MatrixDouble>(shared_from_this(),
260  }
261 
263  return boost::shared_ptr<MatrixDouble>(shared_from_this(),
265  }
266 
268  return boost::shared_ptr<MatrixDouble>(shared_from_this(),
270  }
271 
273  return boost::shared_ptr<MatrixDouble>(shared_from_this(), &rotAxisAtPts);
274  }
275 
277  return boost::shared_ptr<MatrixDouble>(shared_from_this(),
278  &rotAxisDotAtPts);
279  }
280 
282  return boost::shared_ptr<MatrixDouble>(shared_from_this(), &GAtPts);
283  }
284 
286  return boost::shared_ptr<MatrixDouble>(shared_from_this(), &GAtPts);
287  }
288 
289  std::array<MatrixDouble, 6> expLogUHessian;
290  std::array<MatrixDouble, 3> rotationHessian;
291 
292  // Not really data at integration points, used to calculate Schur complement
293 
294  struct BlockMatData {
295 
296  std::string rowField;
297  std::string colField;
298  EntityType rowType;
299  EntityType colType;
300  int rowSide;
301  int colSide;
305 
307 
308  BlockMatData(const std::string row_field, const std::string col_field,
309  EntityType row_type, EntityType col_type, int row_side,
310  int col_side, const MatrixDouble &m, const VectorInt row_ind,
311  VectorInt col_ind)
312  : rowField(row_field), colField(col_field), rowType(row_type),
313  colType(col_type), rowSide(row_side), colSide(col_side),
314  setAtElement(true) {
315 
316  M.resize(m.size1(), m.size2(), false);
317  noalias(M) = m;
318  rowInd.resize(row_ind.size(), false);
319  noalias(rowInd) = row_ind;
320  colInd.resize(col_ind.size(), false);
321  noalias(colInd) = col_ind;
322  }
323 
324  void setInd(const VectorInt &row_ind, const VectorInt &col_ind) const {
325  auto &const_row_ind = const_cast<VectorInt &>(rowInd);
326  auto &const_col_ind = const_cast<VectorInt &>(colInd);
327  const_row_ind.resize(row_ind.size(), false);
328  noalias(const_row_ind) = row_ind;
329  const_col_ind.resize(col_ind.size(), false);
330  noalias(const_col_ind) = col_ind;
331  }
332 
333  void setMat(const MatrixDouble &m) const {
334  auto &const_m = const_cast<MatrixDouble &>(M);
335  const_m.resize(m.size1(), m.size2(), false);
336  noalias(const_m) = m;
337  }
338 
339  void addMat(const MatrixDouble &m) const {
340  auto &const_m = const_cast<MatrixDouble &>(M);
341  const_m += m;
342  }
343 
344  void clearMat() const {
345  auto &const_m = const_cast<MatrixDouble &>(M);
346  const_m.clear();
347  }
348 
349  void setSetAtElement() const {
350  bool &set = const_cast<bool &>(setAtElement);
351  set = true;
352  }
353 
354  void unSetAtElement() const {
355  bool &set = const_cast<bool &>(setAtElement);
356  set = false;
357  }
358  };
359 
360  typedef multi_index_container<
361  BlockMatData,
362  indexed_by<
363 
364  ordered_unique<
365 
366  composite_key<
367  BlockMatData,
368 
369  member<BlockMatData, std::string, &BlockMatData::rowField>,
370  member<BlockMatData, std::string, &BlockMatData::colField>,
371  member<BlockMatData, EntityType, &BlockMatData::rowType>,
372  member<BlockMatData, EntityType, &BlockMatData::colType>,
373  member<BlockMatData, int, &BlockMatData::rowSide>,
374  member<BlockMatData, int, &BlockMatData::colSide>
375 
376  >>,
377  ordered_non_unique<
378 
379  composite_key<
380  BlockMatData,
381 
382  member<BlockMatData, std::string, &BlockMatData::rowField>,
383  member<BlockMatData, std::string, &BlockMatData::colField>,
384  member<BlockMatData, EntityType, &BlockMatData::rowType>,
385  member<BlockMatData, EntityType, &BlockMatData::colType>
386 
387  >>,
388  ordered_non_unique<
389 
390  composite_key<
391  BlockMatData,
392 
393  member<BlockMatData, std::string, &BlockMatData::rowField>,
394  member<BlockMatData, std::string, &BlockMatData::colField>
395 
396  >>,
397  ordered_non_unique<
398  member<BlockMatData, std::string, &BlockMatData::rowField>>,
399  ordered_non_unique<
400  member<BlockMatData, std::string, &BlockMatData::colField>>
401 
402  >>
404 
408 
409  boost::shared_ptr<PhysicalEquations> physicsPtr;
410 };
411 
412 struct OpJacobian;
413 
415 
420 
423 
424  PhysicalEquations() = delete;
425  PhysicalEquations(const int size_active, const int size_dependent)
426  : activeVariables(size_active, 0), dependentVariables(size_dependent, 0),
427  dependentVariablesDirevatives(size_dependent * size_active, 0) {}
428  virtual ~PhysicalEquations() {}
429 
430  virtual MoFEMErrorCode recordTape(const int tag, DTensor2Ptr *t_h) = 0;
431 
432  virtual OpJacobian *
433  returnOpJacobian(const std::string &field_name, const int tag,
434  const bool eval_rhs, const bool eval_lhs,
435  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
436  boost::shared_ptr<PhysicalEquations> &physics_ptr) = 0;
437 
438  std::vector<double> activeVariables;
439  std::vector<double> dependentVariables;
440  std::vector<double> dependentVariablesDirevatives;
441 
442  /** \name Active variables */
443 
444  /**@{*/
445 
446  template <int S>
447  inline static DTensor2Ptr get_VecTensor2(std::vector<double> &v) {
448  return DTensor2Ptr(&v[S + 0], &v[S + 1], &v[S + 2], &v[S + 3], &v[S + 4],
449  &v[S + 5], &v[S + 6], &v[S + 7], &v[S + 8]);
450  }
451 
452  template <int S>
453  inline static DTensor0Ptr get_VecTensor0(std::vector<double> &v) {
454  return DTensor0Ptr(&v[S + 0]);
455  }
456 
457  template <int S0, int S1, int S2>
458  inline static DTensor3Ptr get_vecTensor3(std::vector<double> &v) {
459  constexpr int A00 = 18 * S0 + 18 * 0 + 9 * S1 + 3 * S2;
460  constexpr int A01 = 18 * S0 + 18 * 1 + 9 * S1 + 3 * S2;
461  constexpr int A02 = 18 * S0 + 18 * 2 + 9 * S1 + 3 * S2;
462  constexpr int A10 = 18 * S0 + 18 * 3 + 9 * S1 + 3 * S2;
463  constexpr int A11 = 18 * S0 + 18 * 4 + 9 * S1 + 3 * S2;
464  constexpr int A12 = 18 * S0 + 18 * 5 + 9 * S1 + 3 * S2;
465  constexpr int A20 = 18 * S0 + 18 * 6 + 9 * S1 + 3 * S2;
466  constexpr int A21 = 18 * S0 + 18 * 7 + 9 * S1 + 3 * S2;
467  constexpr int A22 = 18 * S0 + 18 * 8 + 9 * S1 + 3 * S2;
468  return DTensor3Ptr(
469 
470  &v[A00 + 0], &v[A00 + 1], &v[A00 + 2], &v[A01 + 0], &v[A01 + 1],
471  &v[A01 + 2], &v[A02 + 0], &v[A02 + 1], &v[A02 + 2],
472 
473  &v[A10 + 0], &v[A10 + 1], &v[A10 + 2], &v[A11 + 0], &v[A11 + 1],
474  &v[A11 + 2], &v[A12 + 0], &v[A12 + 1], &v[A12 + 2],
475 
476  &v[A20 + 0], &v[A20 + 1], &v[A20 + 2], &v[A21 + 0], &v[A21 + 1],
477  &v[A21 + 2], &v[A22 + 0], &v[A22 + 1], &v[A22 + 2]
478 
479  );
480  }
481 
482  /**@}*/
483 
484  /** \name Active variables */
485 
486  /**@{*/
487 
488  inline DTensor2Ptr get_h() { return get_VecTensor2<0>(activeVariables); }
489  inline DTensor2Ptr get_H() { return get_VecTensor2<9>(activeVariables); }
490 
491  /**@}*/
492 
493  /** \name Dependent variables */
494 
495  /**@{*/
496 
497  inline DTensor2Ptr get_P() { return get_VecTensor2<0>(dependentVariables); }
499  return get_VecTensor2<9>(dependentVariables);
500  }
501  inline DTensor0Ptr get_Phi() {
502  return get_VecTensor0<18>(dependentVariables);
503  }
504  inline double &get_PhiRef() { return dependentVariables[18]; }
506  return get_VecTensor2<9 + 9 + 1>(dependentVariables);
507  }
508 
509  /**@}*/
510 
511  /** \name Derivatives of dependent variables */
512 
513  /**@{*/
514 
516  return get_vecTensor3<0, 0, 0>(dependentVariablesDirevatives);
517  }
519  return get_vecTensor3<0, 1, 0>(dependentVariablesDirevatives);
520  }
522  return get_vecTensor3<0, 0, 1>(dependentVariablesDirevatives);
523  }
525  return get_vecTensor3<0, 1, 1>(dependentVariablesDirevatives);
526  }
528  return get_vecTensor3<0, 0, 2>(dependentVariablesDirevatives);
529  }
531  return get_vecTensor3<0, 1, 2>(dependentVariablesDirevatives);
532  }
533 
535  return get_vecTensor3<9, 0, 0>(dependentVariablesDirevatives);
536  }
538  return get_vecTensor3<9, 1, 0>(dependentVariablesDirevatives);
539  }
541  return get_vecTensor3<9, 0, 1>(dependentVariablesDirevatives);
542  }
544  return get_vecTensor3<9, 1, 1>(dependentVariablesDirevatives);
545  }
547  return get_vecTensor3<9, 0, 2>(dependentVariablesDirevatives);
548  }
550  return get_vecTensor3<9, 1, 2>(dependentVariablesDirevatives);
551  }
552 
554  return get_VecTensor2<18 * (9 + 9) + 0>(dependentVariablesDirevatives);
555  }
557  return get_VecTensor2<18 * (9 + 9) + 9>(dependentVariablesDirevatives);
558  }
559 
561  return get_vecTensor3<9 + 9 + 1, 0, 0>(dependentVariablesDirevatives);
562  }
564  return get_vecTensor3<9 + 9 + 1, 1, 0>(dependentVariablesDirevatives);
565  }
567  return get_vecTensor3<9 + 9 + 1, 0, 1>(dependentVariablesDirevatives);
568  }
570  return get_vecTensor3<9 + 9 + 1, 1, 1>(dependentVariablesDirevatives);
571  }
573  return get_vecTensor3<9 + 9 + 1, 0, 2>(dependentVariablesDirevatives);
574  }
576  return get_vecTensor3<9 + 9 + 1, 1, 2>(dependentVariablesDirevatives);
577  }
578 
579  /**@}*/
580 };
581 
582 struct BcDisp {
583  BcDisp(std::string name, std::vector<double> &attr, Range &faces);
584  std::string blockName;
585  Range faces;
588 };
589 typedef std::vector<BcDisp> BcDispVec;
590 
591 struct BcRot {
592  BcRot(std::string name, std::vector<double> &attr, Range &faces);
593  std::string blockName;
594  Range faces;
596  double theta;
597 };
598 typedef std::vector<BcRot> BcRotVec;
599 
600 typedef std::vector<Range> TractionFreeBc;
601 
602 struct TractionBc {
603  TractionBc(std::string name, std::vector<double> &attr, Range &faces);
604  std::string blockName;
605  Range faces;
608 };
609 typedef std::vector<TractionBc> TractionBcVec;
610 
611 struct OpJacobian : public UserDataOperator {
612  const int tAg; ///< adol-c tape
613  const bool evalRhs;
614  const bool evalLhs;
615  boost::shared_ptr<DataAtIntegrationPts>
616  dataAtPts; ///< data at integration pts
617  boost::shared_ptr<PhysicalEquations>
618  physicsPtr; ///< material physical equations
619 
620  OpJacobian(const std::string &field_name, const int tag, const bool eval_rhs,
621  const bool eval_lhs,
622  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
623  boost::shared_ptr<PhysicalEquations> &physics_ptr)
624  : UserDataOperator(field_name, OPROW), tAg(tag), evalRhs(eval_rhs),
625  evalLhs(eval_lhs), dataAtPts(data_ptr), physicsPtr(physics_ptr) {}
626 
627  virtual MoFEMErrorCode evaluateRhs(EntData &data) = 0;
628  virtual MoFEMErrorCode evaluateLhs(EntData &data) = 0;
629 
630  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
631 };
632 
633 template <typename T> struct OpAssembleBasic : public T {
634 
635  const bool assembleSymmetry;
636 
637  boost::shared_ptr<DataAtIntegrationPts>
638  dataAtPts; ///< data at integration pts
639 
640  OpAssembleBasic(const std::string &field_name,
641  boost::shared_ptr<DataAtIntegrationPts> data_ptr,
642  const char type)
643  : T(field_name, type), dataAtPts(data_ptr), assembleSymmetry(false) {}
644 
645  OpAssembleBasic(const std::string &row_field, const std::string &col_field,
646  boost::shared_ptr<DataAtIntegrationPts> data_ptr,
647  const char type, const bool assemble_symmetry)
648  : T(row_field, col_field, type, false), dataAtPts(data_ptr),
649  assembleSymmetry(assemble_symmetry) {}
650 
651  VectorDouble nF; ///< local right hand side vector
652  MatrixDouble K; ///< local tangent matrix
654 
657  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
659  }
660 
661  virtual MoFEMErrorCode integrate(int row_side, EntityType row_type,
662  EntData &data) {
664  CHKERR integrate(data);
666  }
667 
668  virtual MoFEMErrorCode assemble(EntData &data) {
670  double *vec_ptr = &*nF.begin();
671  int nb_dofs = data.getIndices().size();
672  int *ind_ptr = &*data.getIndices().begin();
673  CHKERR VecSetValues(this->getTSf(), nb_dofs, ind_ptr, vec_ptr, ADD_VALUES);
675  }
676 
677  virtual MoFEMErrorCode assemble(int row_side, EntityType row_type,
678  EntData &data) {
680  CHKERR assemble(data);
682  }
683 
684  virtual MoFEMErrorCode integrate(EntData &row_data, EntData &col_data) {
686  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
688  }
689 
690  virtual MoFEMErrorCode assemble(int row_side, int col_side,
691  EntityType row_type, EntityType col_type,
692  EntData &row_data, EntData &col_data) {
694 
695  auto &bmc = dataAtPts->blockMatContainor;
696  int *row_ind_ptr = &*row_data.getIndices().begin();
697  int *col_ind_ptr = &*col_data.getIndices().begin();
698  int row_nb_dofs = row_data.getIndices().size();
699  int col_nb_dofs = col_data.getIndices().size();
700 
701  auto add_block = [&](auto &row_name, auto &col_name, auto row_side,
702  auto col_side, auto row_type, auto col_type,
703  const auto &m, const auto &row_ind,
704  const auto &col_ind) {
705  auto it = bmc.get<0>().find(boost::make_tuple(
706  row_name, col_name, row_type, col_type, row_side, col_side));
707  if (it != bmc.get<0>().end()) {
708  it->setMat(m);
709  it->setInd(row_ind, col_ind);
710  it->setSetAtElement();
711  } else
713  row_name, col_name, row_type, col_type, row_side, col_side, m,
714  row_ind, col_ind));
715  };
716 
717  add_block(this->rowFieldName, this->colFieldName, row_side, col_side,
718  row_type, col_type, K, row_data.getIndices(),
719  col_data.getIndices());
720  if (assembleSymmetry) {
721  transposeK.resize(col_nb_dofs, row_nb_dofs, false);
722  noalias(transposeK) = trans(K);
723  add_block(this->colFieldName, this->rowFieldName, col_side, row_side,
724  col_type, row_type, transposeK, col_data.getIndices(),
725  row_data.getIndices());
726  }
727 
728  double *mat_ptr = &*K.data().begin();
729  CHKERR MatSetValues(this->getTSB(), row_nb_dofs, row_ind_ptr, col_nb_dofs,
730  col_ind_ptr, mat_ptr, ADD_VALUES);
731  if (assembleSymmetry) {
732  double *mat_ptr = &*transposeK.data().begin();
733  CHKERR MatSetValues(this->getTSB(), col_nb_dofs, col_ind_ptr,
734  row_nb_dofs, row_ind_ptr, mat_ptr, ADD_VALUES);
735  }
736 
738  }
739 
740  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
742  if (data.getIndices().empty())
744  nF.resize(data.getIndices().size(), false);
745  nF.clear();
746  CHKERR integrate(side, type, data);
747  CHKERR assemble(side, type, data);
749  }
750 
751  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
752  EntityType col_type, EntData &row_data,
753  EntData &col_data) {
755  if (row_data.getIndices().empty())
757  if (col_data.getIndices().empty())
759  K.resize(row_data.getIndices().size(), col_data.getIndices().size(), false);
760  K.clear();
761  CHKERR integrate(row_data, col_data);
762  CHKERR assemble(row_side, col_side, row_type, col_type, row_data, col_data);
764  }
765 };
766 
767 struct OpAssembleVolume : public OpAssembleBasic<VolUserDataOperator> {
768  OpAssembleVolume(const std::string &field,
769  boost::shared_ptr<DataAtIntegrationPts> data_ptr,
770  const char type)
771  : OpAssembleBasic<VolUserDataOperator>(field, data_ptr, type) {}
772 
773  OpAssembleVolume(const std::string &row_field, const std::string &col_field,
774  boost::shared_ptr<DataAtIntegrationPts> data_ptr,
775  const char type, const bool assemble_symmetry)
776  : OpAssembleBasic<VolUserDataOperator>(row_field, col_field, data_ptr,
777  type, assemble_symmetry) {}
778 };
779 
780 struct OpAssembleFace : public OpAssembleBasic<FaceUserDataOperator> {
781  OpAssembleFace(const std::string &field,
782  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
783  const char type)
784  : OpAssembleBasic<FaceUserDataOperator>(field, data_ptr, type) {}
785 
786  OpAssembleFace(const std::string &row_field, const std::string &col_field,
787  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
788  const char type, const bool assemble_symmetry)
789  : OpAssembleBasic<FaceUserDataOperator>(row_field, col_field, data_ptr,
790  type, assemble_symmetry) {}
791 };
792 
794  boost::shared_ptr<DataAtIntegrationPts>
795  dataAtPts; ///< data at integration pts
797  const std::string &field_name,
798  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
799  : VolUserDataOperator(field_name, OPROW), dataAtPts(data_ptr) {}
800  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
801 };
802 
804  const double alpha_w;
805  OpSpatialEquilibrium(const std::string &field_name,
806  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
807  const double alpha)
808  : OpAssembleVolume(field_name, data_ptr, OPROW), alpha_w(alpha) {}
810 };
811 
813  OpSpatialRotation(const std::string &field_name,
814  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
815  : OpAssembleVolume(field_name, data_ptr, OPROW) {}
817 };
818 
820  const double alpha_u;
821  OpSpatialPhysical(const std::string &field_name,
822  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
823  const double alpha)
824  : OpAssembleVolume(field_name, data_ptr, OPROW), alpha_u(alpha) {}
826 };
827 
829  OpSpatialConsistencyP(const std::string &field_name,
830  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
831  : OpAssembleVolume(field_name, data_ptr, OPROW) {}
833 };
834 
836  OpSpatialConsistencyBubble(const std::string &field_name,
837  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
838  : OpAssembleVolume(field_name, data_ptr, OPROW) {}
840 };
841 
843  OpSpatialConsistencyDivTerm(const std::string &field_name,
844  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
845  : OpAssembleVolume(field_name, data_ptr, OPROW) {}
847 };
848 
849 struct OpDispBc : public OpAssembleFace {
850  boost::shared_ptr<BcDispVec> bcDispPtr;
851  OpDispBc(const std::string &field_name,
852  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
853  boost::shared_ptr<BcDispVec> &bc_disp_ptr)
854  : OpAssembleFace(field_name, data_ptr, OPROW), bcDispPtr(bc_disp_ptr) {}
856 };
857 
858 struct OpDispBc_dx : public OpAssembleFace {
859  boost::shared_ptr<BcDispVec> bcDispPtr;
860  OpDispBc_dx(const std::string &row_field_name,
861  const std::string &col_field_name,
862  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
863  boost::shared_ptr<BcDispVec> &bc_disp_ptr)
864  : OpAssembleFace(row_field_name, col_field_name, data_ptr, OPROWCOL,
865  false),
866  bcDispPtr(bc_disp_ptr) {}
867  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
868 };
869 
870 struct OpRotationBc : public OpAssembleFace {
871  boost::shared_ptr<BcRotVec> bcRotPtr;
872  OpRotationBc(const std::string &field_name,
873  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
874  boost::shared_ptr<BcRotVec> &bc_rot_ptr)
875  : OpAssembleFace(field_name, data_ptr, OPROW), bcRotPtr(bc_rot_ptr) {}
877 };
878 
880  boost::shared_ptr<BcRotVec> bcRotPtr;
881  OpRotationBc_dx(const std::string &row_field_name,
882  const std::string &col_field_name,
883  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
884  boost::shared_ptr<BcRotVec> &bc_rot_ptr)
885  : OpAssembleFace(row_field_name, col_field_name, data_ptr, OPROWCOL,
886  false),
887  bcRotPtr(bc_rot_ptr) {}
888  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
889 };
890 
891 struct FeTractionBc;
892 
894  OpTractionBc(std::string row_field, FeTractionBc &bc_fe)
895  : FaceUserDataOperator(row_field, FaceUserDataOperator::OPROW),
896  bcFe(bc_fe) {}
897  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
898 
899 protected:
903 };
904 
905 struct FeTractionBc : public FEMethod {
907  FeTractionBc(MoFEM::Interface &m_field, const std::string field_name,
908  boost::shared_ptr<TractionBcVec> &bc)
909  : FEMethod(), mField(m_field), bcData(bc), fieldName(field_name) {}
911  MoFEMErrorCode postProcess() { return 0; }
912 
913  friend MoFEMErrorCode OpTractionBc::doWork(int side, EntityType type,
914  EntData &data);
915 
916 protected:
918  boost::shared_ptr<TractionBcVec> bcData;
919  std::string fieldName;
920 };
921 
922 template <>
924  EpElement(MoFEM::Interface &m_field, const std::string field_name,
925  boost::shared_ptr<TractionBcVec> &bc)
926  : FeTractionBc(m_field, field_name, bc), EpElementBase() {}
931  }
932 };
933 
935  OpSpatialEquilibrium_dw_dP(const std::string &row_field,
936  const std::string &col_field,
937  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
938  const bool assemble_off = false)
939  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
940  assemble_off) {
941  sYmm = false;
942  }
943  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
944 };
945 
947  const double alpha_w;
948  OpSpatialEquilibrium_dw_dw(const std::string &row_field,
949  const std::string &col_field,
950  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
951  const double alpha)
952  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
953  false), alpha_w(alpha) {
954  sYmm = true;
955  }
956  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
957 };
958 
960  const double alpha_u;
961  OpSpatialPhysical_du_du(const std::string &row_field,
962  const std::string &col_field,
963  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
964  const double alpha)
965  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false),
966  alpha_u(alpha) {
967  sYmm = false;
968  }
969  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
970 };
971 
973  OpSpatialPhysical_du_dP(const std::string &row_field,
974  const std::string &col_field,
975  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
976  const bool assemble_off = false)
977  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
978  assemble_off) {
979  sYmm = false;
980  }
981 
982  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
983 };
984 
987  const std::string &row_field, const std::string &col_field,
988  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
989  const bool assemble_off = false)
990  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
991  assemble_off) {
992  sYmm = false;
993  }
994 
995  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
996 };
997 
999  OpSpatialPhysical_du_domega(const std::string &row_field,
1000  const std::string &col_field,
1001  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1002  const bool assemble_off = false)
1003  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
1004  assemble_off) {
1005  sYmm = false;
1006  }
1007  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1008 };
1009 
1011  OpSpatialPhysical_du_dx(const std::string &row_field,
1012  const std::string &col_field,
1013  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1014  const bool assemble_off = false)
1015  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
1016  assemble_off) {
1017  sYmm = false;
1018  }
1019  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1020 };
1021 
1023  OpSpatialRotation_domega_dP(const std::string &row_field,
1024  const std::string &col_field,
1025  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1026  const bool assemble_off)
1027  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL,
1028  assemble_off) {
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)
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 
1048  const std::string &row_field, 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  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false) {
1063  sYmm = false;
1064  }
1065  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1066 };
1067 
1070  const std::string &row_field, const std::string &col_field,
1071  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
1072  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false) {
1073  sYmm = false;
1074  }
1075  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1076 };
1077 
1079 
1081  std::vector<EntityHandle> &mapGaussPts;
1082  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
1083 
1085  std::vector<EntityHandle> &map_gauss_pts,
1086  const std::string field_name,
1087  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
1088  : VolUserDataOperator(field_name, UserDataOperator::OPROW),
1089  postProcMesh(post_proc_mesh), mapGaussPts(map_gauss_pts),
1090  dataAtPts(data_ptr) {}
1091 
1092  MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
1093 };
1094 
1096  OpSpatialPrj(const std::string &row_field,
1097  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
1098  : OpAssembleVolume(row_field, data_ptr, OPROW) {}
1099  MoFEMErrorCode integrate(EntData &row_data);
1100 };
1101 
1103  OpSpatialPrj_dx_dx(const std::string &row_field, const std::string &col_field,
1104  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
1105  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false) {
1106  // FIXME: That is symmetric
1107  sYmm = false;
1108  }
1109  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1110 };
1111 
1113  OpSpatialPrj_dx_dw(const std::string &row_field, const std::string &col_field,
1114  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
1115  : OpAssembleVolume(row_field, col_field, data_ptr, OPROWCOL, false) {
1116  sYmm = false;
1117  }
1118  MoFEMErrorCode integrate(EntData &row_data, EntData &col_data);
1119 };
1120 
1122  OpSpatialSchurBegin(const std::string &row_field,
1123  boost::shared_ptr<DataAtIntegrationPts> &data_ptr)
1124  : OpAssembleVolume(row_field, data_ptr, OPROW) {
1125  sYmm = false;
1126  }
1127  MoFEMErrorCode integrate(EntData &row_data) { return 0; }
1128  MoFEMErrorCode assemble(int row_side, EntityType row_type, EntData &data);
1129 };
1130 
1132  OpSpatialPreconditionMass(const std::string &row_field, MatrixPtr m_ptr)
1133  : OpAssembleVolume(row_field, nullptr, OPROW), mPtr(m_ptr) {
1134  sYmm = false;
1135  }
1136  MoFEMErrorCode integrate(EntData &row_data);
1137  MoFEMErrorCode assemble(int row_side, EntityType row_type, EntData &data) {
1138  return 0;
1139  }
1140 
1141 private:
1143 };
1144 
1146  OpSpatialSchurEnd(const std::string &row_field,
1147  boost::shared_ptr<DataAtIntegrationPts> &data_ptr,
1148  const double eps)
1149  : OpAssembleVolume(row_field, data_ptr, OPROW), eps(eps) {
1150  sYmm = false;
1151  }
1152 
1153  MoFEMErrorCode integrate(EntData &row_data) { return 0; }
1154  MoFEMErrorCode assemble(int side, EntityType type, EntData &data);
1155  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
1156  return assemble(side, type, data);
1157  }
1158 
1159 private:
1160  const double eps;
1161 
1165 
1166  map<std::string, MatrixDouble> invBlockMat;
1167  map<std::string, DataAtIntegrationPts::BlockMatContainor> blockMat;
1168 };
1169 
1171 
1172  OpCalculateStrainEnergy(const std::string field_name,
1173  boost::shared_ptr<DataAtIntegrationPts> data_at_pts,
1174  boost::shared_ptr<double> &e)
1175  : VolUserDataOperator(field_name, VolUserDataOperator::OPROW),
1176  dataAtPts(data_at_pts), energy(e) {}
1177 
1178  MoFEMErrorCode doWork(int side, EntityType type,
1180 
1181 private:
1182  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
1183  boost::shared_ptr<double> energy;
1184 };
1185 
1187 
1189 
1190  MoFEMErrorCode doWork(int side, EntityType type,
1192 };
1193 
1195 
1196  /**
1197  * \brief Getting interface of core database
1198  * @param uuid unique ID of interface
1199  * @param iface returned pointer to interface
1200  * @return error code
1201  */
1202  MoFEMErrorCode query_interface(const MOFEMuuid &uuid,
1203  UnknownInterface **iface) const;
1204 
1206 
1207  boost::shared_ptr<DataAtIntegrationPts> dataAtPts;
1208  boost::shared_ptr<PhysicalEquations> physicalEquations;
1209 
1210  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> elasticFeRhs;
1211  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> elasticFeLhs;
1212  boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> elasticBcLhs;
1213  boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> elasticBcRhs;
1214  boost::shared_ptr<EpFEMethod> schurAssembly;
1215 
1216  SmartPetscObj<DM> dM; ///< Coupled problem all fields
1217  SmartPetscObj<DM> dmElastic; ///< Elastic problem
1218  SmartPetscObj<DM> dmElasticSchurStreach; ///< Sub problem of dmElastic Schur
1219  SmartPetscObj<DM> dmElasticSchurBubble; ///< Sub problem of dmElastic Schur
1220  SmartPetscObj<DM>
1221  dmElasticSchurSpatialDisp; ///< Sub problem of dmElastic Schur
1222  SmartPetscObj<DM> dmElasticSchurOmega; ///< Sub problem of dmElastic Schur
1223  SmartPetscObj<DM> dmMaterial; ///< Material problem
1224 
1225  const std::string piolaStress;
1226  const std::string eshelbyStress;
1227  const std::string spatialDisp;
1228  const std::string materialDisp;
1229  const std::string streachTensor;
1230  const std::string rotAxis;
1231  const std::string materialGradient;
1232  const std::string tauField;
1233  const std::string lambdaField;
1234  const std::string bubbleField;
1235 
1236  const std::string elementVolumeName;
1237  const std::string naturalBcElement;
1238  const std::string essentialBcElement;
1239 
1240  EshelbianCore(MoFEM::Interface &m_field);
1241  virtual ~EshelbianCore();
1242 
1245  double alpha_u;
1246  double alpha_w;
1248 
1250 
1251  boost::shared_ptr<BcDispVec> bcSpatialDispVecPtr;
1252  boost::shared_ptr<BcRotVec> bcSpatialRotationVecPtr;
1253  boost::shared_ptr<TractionBcVec> bcSpatialTraction;
1254  boost::shared_ptr<TractionFreeBc> bcSpatialFreeTraction;
1255 
1256  template <typename BC>
1257  MoFEMErrorCode getBc(boost::shared_ptr<BC> &bc_vec_ptr,
1258  const std::string block_set_name,
1259  const int nb_attributes) {
1262  auto block_name = it->getName();
1263  if (block_name.compare(0, block_set_name.length(), block_set_name) == 0) {
1264  std::vector<double> block_attributes;
1265  CHKERR it->getAttributes(block_attributes);
1266  if (block_attributes.size() != nb_attributes) {
1267  SETERRQ3(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1268  "In block %s expected %d attributes, but given %d",
1269  it->getName().c_str(), nb_attributes,
1270  block_attributes.size());
1271  }
1272  Range faces;
1273  CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
1274  true);
1275  bc_vec_ptr->emplace_back(block_name, block_attributes, faces);
1276  }
1277  }
1279  }
1280 
1282  bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
1283  return getBc(bcSpatialDispVecPtr, "SPATIAL_DISP_BC", 6);
1284  }
1285 
1287  bcSpatialRotationVecPtr = boost::make_shared<BcRotVec>();
1288  return getBc(bcSpatialRotationVecPtr, "SPATIAL_ROTATION_BC", 4);
1289  }
1290 
1292  bcSpatialTraction = boost::make_shared<TractionBcVec>();
1293  return getBc(bcSpatialTraction, "SPATIAL_TRACTION_BC", 6);
1294  }
1295 
1297  boost::shared_ptr<TractionFreeBc> &bc_ptr,
1298  const std::string disp_block_set_name,
1299  const std::string rot_block_set_name);
1300  inline MoFEMErrorCode
1303  boost::shared_ptr<TractionFreeBc>(new TractionFreeBc());
1304  return getTractionFreeBc(meshset, bcSpatialFreeTraction, "SPATIAL_DISP_BC",
1305  "SPATIAL_ROTATION_BC");
1306  }
1307 
1308  MoFEMErrorCode addFields(const EntityHandle meshset = 0);
1311  MoFEMErrorCode addDMs(const BitRefLevel bit = BitRefLevel().set(0));
1312 
1314  const double lambda,
1315  const double mu,
1316  const double sigma_y);
1317 
1318  MoFEMErrorCode addMaterial_HMHMooneyRivlin(const int tape, const double alpha,
1319  const double beta,
1320  const double lambda,
1321  const double sigma_y);
1322 
1324  const int tag, const bool do_rhs, const bool do_lhs,
1325  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe);
1326 
1328  const int tag, const bool add_elastic, const bool add_material,
1329  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe_rhs,
1330  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe_lhs);
1331 
1333  const bool add_elastic, const bool add_material,
1334  boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> &fe_rhs,
1335  boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> &fe_lhs);
1336 
1337  MoFEMErrorCode setElasticElementOps(const int tag);
1339 
1340  MoFEMErrorCode setUpTSElastic(TS ts, Mat m, Vec f);
1341  MoFEMErrorCode solveElastic(TS ts, Vec x);
1342 
1343  MoFEMErrorCode postProcessResults(const int tag, const std::string file);
1344 };
1345 
1346 } // namespace EshelbianPlasticity
1347 
1348 #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)
const double T
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)
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)
Deprecated interface functions.
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)
OpSpatialRotation_domega_dBubble(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const bool assemble_off)
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:509
virtual MoFEMErrorCode assemble(int row_side, EntityType row_type, EntData &data)
std::vector< EntityHandle > & mapGaussPts
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:76
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:485
VectorBoundedArray< int, 3 > VectorInt3
Definition: Types.hpp:87
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
OpSpatialConsistency_dBubble_domega(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr)
MoFEMErrorCode integrate(EntData &data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::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:516
MoFEMErrorCode addBubbleSchurMatrix(Mat SBubble, AO aoSBubble)
virtual MoFEMErrorCode assemble(EntData &data)
std::array< MatrixDouble, 3 > rotationHessian
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
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
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.
static Index< 'm', 3 > m
std::vector< TractionBc > TractionBcVec
OpTractionBc(std::string row_field, FeTractionBc &bc_fe)
virtual MoFEMErrorCode recordTape(const int tag, DTensor2Ptr *t_h)=0
static DTensor2Ptr get_VecTensor2(std::vector< double > &v)
boost::shared_ptr< MatrixDouble > MatrixPtr
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
ForcesAndSourcesCore::UserDataOperator UserDataOperator
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
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:67
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 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
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 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:604
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
MoFEMErrorCode integrate(EntData &row_data, EntData &col_data)
ublas::vector< int, IntAllocator > VectorInt
Definition: Types.hpp:73
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:92
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
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
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)
OpSpatialRotation_domega_dP(const std::string &row_field, const std::string &col_field, boost::shared_ptr< DataAtIntegrationPts > &data_ptr, const bool assemble_off)
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:415
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:74
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
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)
field with C-1 continuity
Definition: definitions.h:180