v0.9.0
EshelbianPlasticity.cpp
Go to the documentation of this file.
1 /**
2  * \file EshelbianPlasticity.cpp
3  * \example EshelbianPlasticity.cpp
4  *
5  * \brief Eshelbian plasticity implementation
6  */
7 
8 #include <MoFEM.hpp>
9 using namespace MoFEM;
10 
11 #include <BasicFiniteElements.hpp>
13 
14 #include <EshelbianPlasticity.hpp>
15 #include <boost/math/constants/constants.hpp>
16 
17 namespace EshelbianPlasticity {
18 
19 MoFEMErrorCode EshelbianCore::query_interface(const MOFEMuuid &uuid,
20  UnknownInterface **iface) const {
22  *iface = NULL;
24  *iface = const_cast<EshelbianCore *>(this);
26  }
27  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
29 }
30 
31 MoFEMErrorCode OpJacobian::doWork(int side, EntityType type, EntData &data) {
33  if (type != MBVERTEX)
35 
36  if (evalRhs)
37  CHKERR evaluateRhs(data);
38 
39  if (evalLhs)
40  CHKERR evaluateLhs(data);
41 
43 }
44 
45 EshelbianCore::EshelbianCore(MoFEM::Interface &m_field)
46  : mField(m_field), dM(PETSC_NULL), dmElastic(PETSC_NULL),
47  dmElasticSchurStreach(PETSC_NULL), dmElasticSchurBubble(PETSC_NULL),
48  dmElasticSchurSpatialDisp(PETSC_NULL), dmElasticSchurOmega(PETSC_NULL),
49  dmMaterial(PETSC_NULL), piolaStress("P"), eshelbyStress("S"),
50  spatialDisp("w"), spatialPosition("x"), materialDisp("W"),
51  streachTensor("u"), rotAxis("omega"), materialGradient("G"),
52  tauField("TAU"), lambdaField("LAMBDA"), bubbleField("BUBBLE"),
53  elementVolumeName("EP"), naturalBcElement("NATURAL_BC"),
54  essentialBcElement("ESSENTIAL_BC") {
55 
56  ierr = getOptions();
57  CHKERRABORT(PETSC_COMM_WORLD, ierr);
58 }
59 
61  auto destroy_dm = [](DM &dm) {
62  if (dm != PETSC_NULL)
63  return DMDestroy(&dm);
64  else
65  return PetscErrorCode(0);
66  };
68  CHKERR destroy_dm(dmElasticSchurOmega);
69  CHKERR destroy_dm(dmElasticSchurBubble);
70  CHKERR destroy_dm(dmElasticSchurStreach);
71  CHKERR destroy_dm(dmMaterial);
72  CHKERR destroy_dm(dmElastic);
73  CHKERR destroy_dm(dM);
74 }
75 
78  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
79  "none");
80 
81  spaceOrder = 1;
82  CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
83  spaceOrder, &spaceOrder, PETSC_NULL);
84  spaceGhostOrder = 1;
85  CHKERR PetscOptionsInt("-space_ghost_order",
86  "approximation ghost frame oder for space", "",
87  spaceGhostOrder, &spaceGhostOrder, PETSC_NULL);
88  ghostFrameOn = PETSC_TRUE;
89  CHKERR PetscOptionsBool("-space_ghost_frame_on", "On ghost frame", "",
90  ghostFrameOn, &ghostFrameOn, PETSC_NULL);
91  materialOrder = 1;
92  CHKERR PetscOptionsInt("-material_order", "approximation oder for material",
93  "", materialOrder, &materialOrder, PETSC_NULL);
94 
95  alpha_u = 0;
96  CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alpha_u,
97  &alpha_u, PETSC_NULL);
98 
99  preconditioner_eps = 1e-3;
100  CHKERR PetscOptionsScalar("-preconditioner_eps", "preconditioner_eps", "",
102  PETSC_NULL);
103 
104  ierr = PetscOptionsEnd();
105  CHKERRG(ierr);
107 }
108 
111 
112  Range tets;
113  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
114  Range tets_skin_part;
115  Skinner skin(&mField.get_moab());
116  CHKERR skin.find_skin(0, tets, false, tets_skin_part);
117  ParallelComm *pcomm =
118  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
119  Range tets_skin;
120  CHKERR pcomm->filter_pstatus(tets_skin_part,
121  PSTATUS_SHARED | PSTATUS_MULTISHARED,
122  PSTATUS_NOT, -1, &tets_skin);
123 
124  auto subtract_faces_where_displacements_are_applied =
125  [&](const std::string disp_block_set_name) {
128  if (it->getName().compare(0, disp_block_set_name.length(),
129  disp_block_set_name) == 0) {
130  Range faces;
131  CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2,
132  faces, true);
133  tets_skin = subtract(tets_skin, faces);
134  }
135  }
137  };
138  CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_DISP_BC");
139  CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_ROTATION_BC");
140  CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_TRACTION_BC");
141 
142  Range faces;
143  CHKERR mField.get_moab().get_adjacencies(tets, 2, true, faces,
144  moab::Interface::UNION);
145  Range faces_not_on_the_skin = subtract(faces, tets_skin);
146 
147  auto add_hdiv_field = [&](const std::string field_name, const int order,
148  const int dim) {
151  MB_TAG_SPARSE, MF_ZERO);
152  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
153  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
154  CHKERR mField.set_field_order(faces_not_on_the_skin, field_name, order);
155  CHKERR mField.set_field_order(tets_skin, field_name, 0);
157  };
158 
159  auto add_hdiv_rt_field = [&](const std::string field_name, const int order,
160  const int dim) {
163  MB_TAG_DENSE, MF_ZERO);
164  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
165  CHKERR mField.set_field_order(meshset, MBTET, field_name, 0);
166  CHKERR mField.set_field_order(tets_skin, field_name, order);
168  };
169 
170  auto add_l2_field = [this, meshset](const std::string field_name,
171  const int order, const int dim) {
174  MB_TAG_DENSE, MF_ZERO);
175  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
176  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
178  };
179 
180  auto add_h1_field = [this, meshset](const std::string field_name,
181  const int order, const int dim) {
184  MB_TAG_DENSE, MF_ZERO);
185  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
186  CHKERR mField.set_field_order(meshset, MBVERTEX, field_name, 1);
187  CHKERR mField.set_field_order(meshset, MBEDGE, field_name, order);
188  CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
189  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
191  };
192 
193  auto add_bubble_field = [this, meshset](const std::string field_name,
194  const int order, const int dim) {
196  CHKERR mField.add_field(field_name, HDIV, USER_BASE, dim, MB_TAG_DENSE,
197  MF_ZERO);
198  // Modify field
199  auto field_ptr = mField.get_field_structure(field_name);
200  auto field_order_table =
201  const_cast<Field *>(field_ptr)->getFieldOrderTable();
202  auto get_cgg_bubble_order_zero = [](int p) { return 0; };
203  auto get_cgg_bubble_order_tet = [](int p) {
204  return NBVOLUMETET_CCG_BUBBLE(p);
205  };
206  field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
207  field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
208  field_order_table[MBTRI] = get_cgg_bubble_order_zero;
209  field_order_table[MBTET] = get_cgg_bubble_order_tet;
210  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
211  CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
212  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
214  };
215 
216  // spatial fields
217  CHKERR add_hdiv_field(piolaStress, spaceOrder, 3);
218  // CHKERR add_hdiv_rt_field(piolaStress + "_RT", spaceOrder, 3);
219  CHKERR add_bubble_field(bubbleField, spaceOrder, 1);
220  CHKERR add_l2_field(spatialDisp, spaceOrder - 1, 3);
221  CHKERR add_l2_field(rotAxis, spaceOrder - 1, 3);
222  CHKERR add_l2_field(streachTensor, spaceOrder, 6);
223  const bool check_if_field_exist = mField.check_field(spatialPosition);
224  CHKERR add_h1_field(spatialPosition, spaceGhostOrder, 3);
225 
226  // material fields
227  CHKERR add_hdiv_field(eshelbyStress, materialOrder, 3);
228  CHKERR add_l2_field(materialGradient, materialOrder - 1, 9);
229  // CHKERR add_l2_field(materialDisp, materialOrder - 1, 3);
230  // CHKERR add_l2_field(tauField, materialOrder - 1, 1);
231  // CHKERR add_l2_field(lambdaField, materialOrder - 1, 1);
232 
233  // Add history filedes
234  CHKERR add_l2_field(materialGradient + "0", materialOrder - 1, 9);
235 
237 
238  if (!check_if_field_exist) {
240  CHKERR mField.loop_dofs(spatialPosition, ent_method);
241  }
242 
244 }
245 
249 
250  // set finite element fields
251  auto add_field_to_fe = [this](const std::string fe,
252  const std::string field_name) {
258  };
259 
263 
264  CHKERR add_field_to_fe(elementVolumeName, piolaStress);
265  CHKERR add_field_to_fe(elementVolumeName, bubbleField);
266  CHKERR add_field_to_fe(elementVolumeName, eshelbyStress);
267  CHKERR add_field_to_fe(elementVolumeName, streachTensor);
268  CHKERR add_field_to_fe(elementVolumeName, rotAxis);
269  CHKERR add_field_to_fe(elementVolumeName, spatialDisp);
270  CHKERR add_field_to_fe(elementVolumeName, spatialPosition);
271  CHKERR add_field_to_fe(elementVolumeName, streachTensor);
272  CHKERR add_field_to_fe(elementVolumeName, materialGradient);
273  // CHKERR add_field_to_fe(elementVolumeName, materialDisp);
274  // CHKERR add_field_to_fe(elementVolumeName, tauField);
275  // CHKERR add_field_to_fe(elementVolumeName, lambdaField);
276 
277  // CHKERR mField.modify_finite_element_add_field_data(elementVolumeName,
278  // piolaStress + "_RT");
280  materialGradient + "0");
281 
282  // build finite elements data structures
284 
286 }
287 
291 
292  auto bc_elements_add_to_range = [&](const std::string disp_block_set_name,
293  Range &r) {
296  if (it->getName().compare(0, disp_block_set_name.length(),
297  disp_block_set_name) == 0) {
298  Range faces;
299  CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
300  true);
301  r.merge(faces);
302  }
303  }
305  };
306  Range natural_bc_elements;
307  CHKERR bc_elements_add_to_range("SPATIAL_DISP_BC", natural_bc_elements);
308  CHKERR bc_elements_add_to_range("SPATIAL_ROTATION_BC", natural_bc_elements);
309  Range essentail_bc_elements;
310  CHKERR bc_elements_add_to_range("SPATIAL_TRACTION_BC", essentail_bc_elements);
311 
312  // set finite element fields
313  auto add_field_to_fe = [this](const std::string fe,
314  const std::string field_name) {
320  };
321 
323  CHKERR mField.add_ents_to_finite_element_by_type(natural_bc_elements, MBTRI,
325  CHKERR add_field_to_fe(naturalBcElement, piolaStress);
326  CHKERR add_field_to_fe(naturalBcElement, eshelbyStress);
327  CHKERR add_field_to_fe(naturalBcElement, spatialPosition);
328 
330  CHKERR mField.add_ents_to_finite_element_by_type(essentail_bc_elements, MBTRI,
332  CHKERR add_field_to_fe(essentialBcElement, piolaStress);
333  CHKERR add_field_to_fe(essentialBcElement, eshelbyStress);
334  CHKERR add_field_to_fe(essentialBcElement, spatialPosition);
335  // CHKERR mField.modify_finite_element_add_field_data(essentialBcElement,
336  // piolaStress + "_RT");
337 
338  // build finite elements data structures
342 }
343 
346 
347  // find adjacencies between finite elements and dofs
349 
350  // Create coupled problem
351  CHKERR DMCreate(mField.get_comm(), &dM);
352  CHKERR DMSetType(dM, "DMMOFEM");
353  CHKERR DMMoFEMCreateMoFEM(dM, &mField, "ESHELBY_PLASTICITY", bit,
354  BitRefLevel().set());
355  CHKERR DMMoFEMSetDestroyProblem(dM, PETSC_TRUE);
356  CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
360  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
361  CHKERR DMSetUp(dM);
362  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
363 
364  auto remove_dofs_on_essential_spatial_stress_boundary =
365  [&](const std::string prb_name) {
367  for (int d : {0, 1, 2})
369  prb_name, piolaStress, (*bcSpatialFreeTraction)[d], d, d, NOISY,
370  true);
372  };
373  CHKERR remove_dofs_on_essential_spatial_stress_boundary("ESHELBY_PLASTICITY");
374 
375  // Create elastic sub-problem
376  CHKERR DMCreate(mField.get_comm(), &dmElastic);
377  CHKERR DMSetType(dmElastic, "DMMOFEM");
378  CHKERR DMMoFEMCreateSubDM(dmElastic, dM, "ELASTIC_PROBLEM");
385  if (ghostFrameOn)
392  CHKERR DMSetUp(dmElastic);
393 
394  // Create elastic streach-problem
396  CHKERR DMSetType(dmElasticSchurStreach, "DMMOFEM");
398  "ELASTIC_PROBLEM_STREACH_SCHUR");
404  if (ghostFrameOn)
406  spatialPosition.c_str());
412  CHKERR DMSetUp(dmElasticSchurStreach);
413 
414  // Create elastic bubble-problem
416  CHKERR DMSetType(dmElasticSchurBubble, "DMMOFEM");
418  "ELASTIC_PROBLEM_BUBBLE_SCHUR");
423  if (ghostFrameOn)
430  CHKERR DMSetUp(dmElasticSchurBubble);
431 
432  // Create elastic omega-problem
434  CHKERR DMSetType(dmElasticSchurOmega, "DMMOFEM");
436  "ELASTIC_PROBLEM_OMEGA_SCHUR");
440  if (ghostFrameOn)
447  CHKERR DMSetUp(dmElasticSchurOmega);
448 
449  // Create elastic tet_stress-problem
451  CHKERR DMSetType(dmElasticSchurSpatialDisp, "DMMOFEM");
453  "ELASTIC_PROBLEM_SPATIAL_DISP_SCHUR");
456  if (ghostFrameOn)
458  spatialPosition.c_str());
460  elementVolumeName.c_str());
463  essentialBcElement.c_str());
467 
468  {
469  PetscSection section;
470  CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
471  &section);
472  CHKERR DMSetDefaultSection(dmElastic, section);
473  CHKERR DMSetDefaultGlobalSection(dmElastic, section);
474  CHKERR PetscSectionDestroy(&section);
475  }
476 
478 }
479 
480 BcDisp::BcDisp(std::string name, std::vector<double> &attr, Range &faces)
481  : blockName(name), faces(faces) {
482  vals.resize(3, false);
483  flags.resize(3, false);
484  for (int ii = 0; ii != 3; ++ii) {
485  vals[ii] = attr[ii];
486  flags[ii] = static_cast<int>(attr[ii + 3]);
487  }
488 }
489 
490 BcRot::BcRot(std::string name, std::vector<double> &attr, Range &faces)
491  : blockName(name), faces(faces) {
492  vals.resize(3, false);
493  for (int ii = 0; ii != 3; ++ii) {
494  vals[ii] = attr[ii];
495  }
496  theta = attr[3];
497 }
498 
499 TractionBc::TractionBc(std::string name, std::vector<double> &attr,
500  Range &faces)
501  : blockName(name), faces(faces) {
502  vals.resize(3, false);
503  flags.resize(3, false);
504  for (int ii = 0; ii != 3; ++ii) {
505  vals[ii] = attr[ii];
506  flags[ii] = static_cast<int>(attr[ii + 3]);
507  }
508 }
509 
512  boost::shared_ptr<TractionFreeBc> &bc_ptr,
513  const std::string disp_block_set_name,
514  const std::string rot_block_set_name) {
516  Range tets;
517  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
518  Range tets_skin_part;
519  Skinner skin(&mField.get_moab());
520  CHKERR skin.find_skin(0, tets, false, tets_skin_part);
521  ParallelComm *pcomm =
522  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
523  Range tets_skin;
524  CHKERR pcomm->filter_pstatus(tets_skin_part,
525  PSTATUS_SHARED | PSTATUS_MULTISHARED,
526  PSTATUS_NOT, -1, &tets_skin);
527 
528  bc_ptr->resize(3);
529  for (int dd = 0; dd != 3; ++dd)
530  (*bc_ptr)[dd] = tets_skin;
531 
533  if (it->getName().compare(0, disp_block_set_name.length(),
534  disp_block_set_name) == 0) {
535  std::vector<double> block_attributes;
536  CHKERR it->getAttributes(block_attributes);
537  if (block_attributes.size() != 6) {
538  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
539  "In block %s six attributes are required for given BC "
540  "blockset (3 values + "
541  "3 flags) != %d",
542  it->getName().c_str(), block_attributes.size());
543  }
544  Range faces;
545  CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
546  true);
547  if (block_attributes[3] != 0)
548  (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
549  if (block_attributes[4] != 0)
550  (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
551  if (block_attributes[5] != 0)
552  (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
553  }
554  if (it->getName().compare(0, rot_block_set_name.length(),
555  rot_block_set_name) == 0) {
556  Range faces;
557  CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
558  true);
559  (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
560  (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
561  (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
562  }
563  }
564 
565  // for (int dd = 0; dd != 3; ++dd) {
566  // EntityHandle meshset;
567  // CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
568  // CHKERR mField.get_moab().add_entities(meshset, (*bc_ptr)[dd]);
569  // std::string file_name = disp_block_set_name +
570  // "_traction_free_bc_" + boost::lexical_cast<std::string>(dd) + ".vtk";
571  // CHKERR mField.get_moab().write_file(file_name.c_str(), " VTK ", "",
572  // &meshset, 1);
573  // CHKERR mField.get_moab().delete_entities(&meshset, 1);
574  // }
575 
577 }
578 
579 /**
580  * @brief Set integration rule on element
581  * \param order on row
582  * \param order on column
583  * \param order on data
584  *
585  * Use maximal oder on data in order to determine integration rule
586  *
587  */
588 struct VolRule {
589  int operator()(int p_row, int p_col, int p_data) const {
590  return 2 * (p_data + 1);
591  }
592 };
593 
594 struct FaceRule {
595  int operator()(int p_row, int p_col, int p_data) const {
596  return 2 * (p_data);
597  }
598 };
599 
601 
604 
606  BaseFunctionUnknownInterface **iface) const {
608  *iface = NULL;
609  if (uuid == IDD_TET_BASE_FUNCTION) {
610  *iface = const_cast<CGGUserPolynomialBase *>(this);
612  } else {
613  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "wrong interference");
614  }
615  CHKERR BaseFunction::query_interface(uuid, iface);
617  }
618 
620  boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
622 
624  CHKERR ctx_ptr->query_interface(IDD_TET_BASE_FUNCTION, &iface);
625  cTx = reinterpret_cast<EntPolynomialBaseCtx *>(iface);
626 
627  int nb_gauss_pts = pts.size2();
628  if (!nb_gauss_pts) {
630  }
631 
632  if (pts.size1() < 3) {
633  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
634  "Wrong dimension of pts, should be at least 3 rows with "
635  "coordinates");
636  }
637 
638  switch (cTx->sPace) {
639  case HDIV:
641  break;
642  default:
643  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
644  }
645 
647  }
648 
649 private:
651 
653 
656 
657  // This should be used only in case USER_BASE is selected
658  if (cTx->bAse != USER_BASE) {
659  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
660  "Wrong base, should be USER_BASE");
661  }
662  // get access to data structures on element
664  // Get approximation order on element. Note that bubble functions are only
665  // on tetrahedron.
666  const int order = data.dataOnEntities[MBTET][0].getDataOrder();
667  /// number of integration points
668  const int nb_gauss_pts = pts.size2();
669  // number of base functions
670 
671  // calculate shape functions, i.e. barycentric coordinates
672  shapeFun.resize(nb_gauss_pts, 4, false);
673  CHKERR ShapeMBTET(&*shapeFun.data().begin(), &pts(0, 0), &pts(1, 0),
674  &pts(2, 0), nb_gauss_pts);
675  // direvatives of shape functions
676  double diff_shape_fun[12];
677  CHKERR ShapeDiffMBTET(diff_shape_fun);
678 
679  const int nb_base_functions = NBVOLUMETET_CCG_BUBBLE(order);
680  // get base functions and set size
681  MatrixDouble &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
682  phi.resize(nb_gauss_pts, 9 * nb_base_functions, false);
683  // finally calculate base functions
685  &phi(0, 0), &phi(0, 1), &phi(0, 2),
686 
687  &phi(0, 3), &phi(0, 4), &phi(0, 5),
688 
689  &phi(0, 6), &phi(0, 7), &phi(0, 8));
690  CHKERR CGG_BubbleBase_MBTET(order, &shapeFun(0, 0), diff_shape_fun, t_phi,
691  nb_gauss_pts);
692 
694  }
695 };
696 
698  const int tag, const bool do_rhs, const bool do_lhs,
699  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe) {
701  fe = boost::make_shared<EpElement<VolumeElementForcesAndSourcesCore>>(mField);
702 
703  fe->getUserPolynomialBase() =
704  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
705 
706  // set integration rule
707  fe->getRuleHook = VolRule();
708 
709  if (!dataAtPts) {
710  dataAtPts =
711  boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
712  dataAtPts->approxPAtPts = boost::make_shared<MatrixDouble>();
713  dataAtPts->divPAtPts = boost::make_shared<MatrixDouble>();
714  dataAtPts->approxSigmaAtPts = boost::make_shared<MatrixDouble>();
715  dataAtPts->divSigmaAtPts = boost::make_shared<MatrixDouble>();
716  dataAtPts->wAtPts = boost::make_shared<MatrixDouble>();
717  dataAtPts->xAtPts = boost::make_shared<MatrixDouble>();
718  dataAtPts->xGradAtPts = boost::make_shared<MatrixDouble>();
719  dataAtPts->xDotGradAtPts = boost::make_shared<MatrixDouble>();
720  dataAtPts->xDetAtPts = boost::make_shared<VectorDouble>();
721  dataAtPts->xInvGradAtPts = boost::make_shared<MatrixDouble>();
722  dataAtPts->rotAxisAtPts = boost::make_shared<MatrixDouble>();
723  dataAtPts->rotDotAxisAtPts = boost::make_shared<MatrixDouble>();
724  dataAtPts->streachTensorAtPts = boost::make_shared<MatrixDouble>();
725  dataAtPts->streachDotTensorAtPts = boost::make_shared<MatrixDouble>();
726  dataAtPts->rotMatAtPts = boost::make_shared<MatrixDouble>();
727  dataAtPts->diffRotMatAtPts = boost::make_shared<MatrixDouble>();
728  dataAtPts->hAtPts = boost::make_shared<MatrixDouble>();
729  dataAtPts->hDeltaAtPts = boost::make_shared<MatrixDouble>();
730  dataAtPts->GAtPts = boost::make_shared<MatrixDouble>();
731  dataAtPts->G0AtPts = boost::make_shared<MatrixDouble>();
732  dataAtPts->wwMatPtr = boost::make_shared<MatrixDouble>();
733  dataAtPts->ooMatPtr = boost::make_shared<MatrixDouble>();
734  }
735 
736  // calculate fields values
737  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
738  piolaStress, dataAtPts->approxPAtPts));
739  // fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
740  // piolaStress + "_RT", dataAtPts->approxPAtPts, MBMAXTYPE));
741  fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
742  bubbleField, dataAtPts->approxPAtPts, MBMAXTYPE));
743  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
744  piolaStress, dataAtPts->divPAtPts));
745  // fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
746  // piolaStress + "_RT", dataAtPts->divPAtPts, MBMAXTYPE));
747  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
748  eshelbyStress, dataAtPts->approxSigmaAtPts));
749  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
750  eshelbyStress, dataAtPts->divSigmaAtPts));
751  fe->getOpPtrVector().push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
752  streachTensor, dataAtPts->streachTensorAtPts, MBTET));
753  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
754  rotAxis, dataAtPts->rotAxisAtPts, MBTET));
755  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
756  rotAxis, dataAtPts->rotDotAxisAtPts, MBTET));
757  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
758  spatialPosition, dataAtPts->xAtPts, MBVERTEX));
759  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
760  spatialPosition, dataAtPts->xGradAtPts));
761  fe->getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
762  materialGradient, dataAtPts->GAtPts, MBTET));
763  fe->getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
764  materialGradient + "0", dataAtPts->G0AtPts, MBTET));
765  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
766  spatialDisp, dataAtPts->wAtPts, MBTET));
767 
768  // velocities
769  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradientDot<3, 3>(
770  spatialPosition, dataAtPts->xDotGradAtPts));
771  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
772  rotAxis, dataAtPts->rotDotAxisAtPts, MBTET));
773  fe->getOpPtrVector().push_back(
775  streachTensor, dataAtPts->streachDotTensorAtPts, MBTET));
776 
777  // calculate other derived quantities
778  fe->getOpPtrVector().push_back(
780 
781  // evaluate integration points
782  fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
783  spatialDisp, tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
784 
786 }
787 
789  const int tag, const bool add_elastic, const bool add_material,
790  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe_rhs,
791  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe_lhs) {
793 
794  // Right hand side
795  CHKERR setBaseVolumeElementOps(tag, true, false, fe_rhs);
796 
797  // elastic
798  if (add_elastic) {
799  fe_rhs->getOpPtrVector().push_back(
801  fe_rhs->getOpPtrVector().push_back(
803  fe_rhs->getOpPtrVector().push_back(
805  fe_rhs->getOpPtrVector().push_back(
807  fe_rhs->getOpPtrVector().push_back(
809  fe_rhs->getOpPtrVector().push_back(
811  fe_rhs->getOpPtrVector().push_back(
813  }
814 
815  // Left hand side
816  CHKERR setBaseVolumeElementOps(tag, true, true, fe_lhs);
817 
818  // elastic
819  if (add_elastic) {
820 
821  // Schur
822  fe_lhs->getOpPtrVector().push_back(
824 
825  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_du(
827  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
829 
830  fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
832  fe_lhs->getOpPtrVector().push_back(
834  fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
835  bubbleField, rotAxis, dataAtPts, false));
836  fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dP_dx(
838  fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_dx(
840 
841  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
843  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
844  streachTensor, rotAxis, dataAtPts, false));
845  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dx(
847 
848  fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dP(
849  rotAxis, piolaStress, dataAtPts, false));
850  fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dBubble(
851  rotAxis, bubbleField, dataAtPts, false));
852  fe_lhs->getOpPtrVector().push_back(
854  fe_lhs->getOpPtrVector().push_back(
856  fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dx(
857  rotAxis, spatialPosition, dataAtPts, false));
858 
859  fe_lhs->getOpPtrVector().push_back(
861  fe_lhs->getOpPtrVector().push_back(
863 
864  // Schur
865  fe_lhs->getOpPtrVector().push_back(
867  fe_lhs->getOpPtrVector().push_back(
869  fe_lhs->getOpPtrVector().push_back(
871 
872  if (add_material) {
873  }
874  }
875 
877 }
878 
880  const bool add_elastic, const bool add_material,
881  boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> &fe_rhs,
882  boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> &fe_lhs) {
884 
885  fe_rhs =
886  boost::make_shared<EpElement<FaceElementForcesAndSourcesCore>>(mField);
887  fe_lhs =
888  boost::make_shared<EpElement<FaceElementForcesAndSourcesCore>>(mField);
889 
890  // set integration rule
891  fe_rhs->getRuleHook = FaceRule();
892  fe_lhs->getRuleHook = FaceRule();
893 
894  if (!dataAtPts->xAtPts)
895  dataAtPts->xAtPts = boost::make_shared<MatrixDouble>();
896 
897  if (add_elastic) {
898  fe_rhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
899  spatialPosition, dataAtPts->xAtPts, MBVERTEX));
900  fe_rhs->getOpPtrVector().push_back(
902  fe_rhs->getOpPtrVector().push_back(
904 
905  fe_lhs->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
906  spatialPosition, dataAtPts->xAtPts, MBVERTEX));
907  fe_lhs->getOpPtrVector().push_back(new OpDispBc_dx(
909  fe_lhs->getOpPtrVector().push_back(new OpRotationBc_dx(
911  }
912 
914 }
915 
919  elasticFeLhs);
922 }
923 
926  boost::shared_ptr<FEMethod> null;
927  boost::shared_ptr<EpElement<FeTractionBc>> spatial_traction_bc(
928  new EpElement<FeTractionBc>(mField, piolaStress/* + "_RT"*/,
930 
931  CHKERR DMMoFEMTSSetIFunction(dm, DM_NO_ELEMENT, null, spatial_traction_bc,
932  null);
936  spatial_traction_bc);
937 
938  schurAssembly = boost::make_shared<EpFEMethod>();
944 }
945 
948  boost::shared_ptr<TsCtx> ts_ctx;
950  CHKERR TSSetIFunction(ts, f, PETSC_NULL, PETSC_NULL);
951  CHKERR TSSetIJacobian(ts, m, m, PETSC_NULL, PETSC_NULL);
952  CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx.get(), PETSC_NULL);
953 
954  auto add_schur_streach_op = [this](auto &list, Mat S, AO aoS) {
956  for (auto &fe : list)
957  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
958  fe_cast->addStreachSchurMatrix(S, aoS);
959  else
960  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
962  };
963 
964  auto add_schur_streach_pre = [this](auto &list, Mat S, AO aoS) {
966  for (auto &fe : list)
967  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
968  fe_cast->addStreachSchurMatrix(S, aoS);
969  else
970  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
972  };
973 
974  auto add_schur_bubble_op = [this](auto &list, Mat S, AO aoS) {
976  for (auto &fe : list)
977  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
978  fe_cast->addBubbleSchurMatrix(S, aoS);
979  else
980  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
982  };
983 
984  auto add_schur_bubble_pre = [this](auto &list, Mat S, AO aoS) {
986  for (auto &fe : list)
987  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
988  fe_cast->addBubbleSchurMatrix(S, aoS);
989  else
990  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
992  };
993 
994  auto add_schur_spatial_disp_op = [this](auto &list, Mat S, AO aoS) {
996  for (auto &fe : list)
997  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
998  fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
999  else
1000  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
1002  };
1003 
1004  auto add_schur_spatial_disp_pre = [this](auto &list, Mat S, AO aoS) {
1006  for (auto &fe : list)
1007  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
1008  fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
1009  else
1010  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
1012  };
1013 
1014  auto add_schur_omega_op = [this](auto &list, Mat S, AO aoS) {
1016  for (auto &fe : list)
1017  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
1018  fe_cast->addOmegaSchurMatrix(S, aoS);
1019  else
1020  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
1022  };
1023 
1024  auto add_schur_omega_pre = [this](auto &list, Mat S, AO ao) {
1026  for (auto &fe : list)
1027  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
1028  fe_cast->addOmegaSchurMatrix(S, ao);
1029  else
1030  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
1032  };
1033 
1034  const MoFEM::Problem *schur_streach_prb_ptr;
1035  CHKERR DMMoFEMGetProblemPtr(dmElasticSchurStreach, &schur_streach_prb_ptr);
1036  if (auto sub_data = schur_streach_prb_ptr->subProblemData) {
1037  Mat Suu;
1038  CHKERR DMCreateMatrix(dmElasticSchurStreach, &Suu);
1039  AO aoSuu;
1040  CHKERR sub_data->getRowMap(&aoSuu);
1041 
1042  CHKERR add_schur_streach_op(ts_ctx->loops_to_do_IJacobian, Suu, aoSuu);
1043  CHKERR add_schur_streach_pre(ts_ctx->preProcess_IJacobian, Suu, aoSuu);
1044  CHKERR add_schur_streach_pre(ts_ctx->postProcess_IJacobian, Suu, aoSuu);
1045 
1046  CHKERR MatDestroy(&Suu);
1047  CHKERR AODestroy(&aoSuu);
1048 
1049  const MoFEM::Problem *schur_bubble_prb_ptr;
1050  CHKERR DMMoFEMGetProblemPtr(dmElasticSchurBubble, &schur_bubble_prb_ptr);
1051  if (auto bubble_data = schur_bubble_prb_ptr->subProblemData) {
1052  Mat SBubble;
1053  CHKERR DMCreateMatrix(dmElasticSchurBubble, &SBubble);
1054  AO aoSBubble;
1055  CHKERR bubble_data->getRowMap(&aoSBubble);
1056 
1057  CHKERR add_schur_bubble_op(ts_ctx->loops_to_do_IJacobian, SBubble,
1058  aoSBubble);
1059  CHKERR add_schur_bubble_pre(ts_ctx->preProcess_IJacobian, SBubble,
1060  aoSBubble);
1061  CHKERR add_schur_bubble_pre(ts_ctx->postProcess_IJacobian, SBubble,
1062  aoSBubble);
1063 
1064  CHKERR MatDestroy(&SBubble);
1065  CHKERR AODestroy(&aoSBubble);
1066 
1067  const MoFEM::Problem *schur_omega_prb_ptr;
1068  CHKERR DMMoFEMGetProblemPtr(dmElasticSchurOmega, &schur_omega_prb_ptr);
1069  if (auto tet_stress_data = schur_omega_prb_ptr->subProblemData) {
1070  Mat SOmega;
1071  CHKERR DMCreateMatrix(dmElasticSchurOmega, &SOmega);
1072  AO aoSOmega;
1073  CHKERR tet_stress_data->getRowMap(&aoSOmega);
1074 
1075  CHKERR add_schur_omega_op(ts_ctx->loops_to_do_IJacobian, SOmega,
1076  aoSOmega);
1077  CHKERR add_schur_omega_pre(ts_ctx->preProcess_IJacobian, SOmega,
1078  aoSOmega);
1079  CHKERR add_schur_omega_pre(ts_ctx->postProcess_IJacobian, SOmega,
1080  aoSOmega);
1081 
1082  const MoFEM::Problem *schur_spatial_disp_prb_ptr;
1084  &schur_spatial_disp_prb_ptr);
1085  if (auto spatial_disp_data =
1086  schur_spatial_disp_prb_ptr->subProblemData) {
1087 
1088  Mat Sw;
1089  CHKERR DMCreateMatrix(dmElasticSchurSpatialDisp, &Sw);
1090  AO aoSw;
1091  CHKERR spatial_disp_data->getRowMap(&aoSw);
1092 
1093  CHKERR add_schur_spatial_disp_op(ts_ctx->loops_to_do_IJacobian, Sw,
1094  aoSw);
1095  CHKERR add_schur_spatial_disp_pre(ts_ctx->preProcess_IJacobian, Sw,
1096  aoSw);
1097  CHKERR add_schur_spatial_disp_pre(ts_ctx->postProcess_IJacobian, Sw,
1098  aoSw);
1099 
1100  CHKERR MatDestroy(&Sw);
1101  CHKERR AODestroy(&aoSw);
1102  } else
1104  "Problem does not have sub-problem data");
1105 
1106  CHKERR MatDestroy(&SOmega);
1107  CHKERR AODestroy(&aoSOmega);
1108 
1109  } else
1111  "Problem does not have sub-problem data");
1112 
1113  } else
1115  "Problem does not have sub-problem data");
1116 
1117  } else
1119  "Problem does not have sub-problem data");
1120 
1121  struct Monitor : public FEMethod {
1122 
1123  using Ele = ForcesAndSourcesCore;
1124  using VolEle = VolumeElementForcesAndSourcesCore;
1126  using SetPtsData = FieldEvaluatorInterface::SetPtsData;
1127 
1128  EshelbianCore &eP;
1129  boost::shared_ptr<SetPtsData> dataFieldEval;
1130  boost::shared_ptr<VolEle> volPostProcEle;
1131  boost::shared_ptr<double> gEnergy;
1132 
1133  Monitor(EshelbianCore &ep)
1134  : eP(ep),
1135  dataFieldEval(ep.mField.getInterface<FieldEvaluatorInterface>()
1136  ->getData<VolEle>()),
1137  volPostProcEle(new VolEle(ep.mField)), gEnergy(new double) {
1138  ierr = ep.mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
1139  dataFieldEval, "EP");
1140  CHKERRABORT(PETSC_COMM_WORLD, ierr);
1141 
1142  auto no_rule = [](int, int, int) { return -1; };
1143 
1144  auto set_element_for_field_eval = [&]() {
1145  boost::shared_ptr<Ele> vol_ele(dataFieldEval->feMethodPtr.lock());
1146  vol_ele->getRuleHook = no_rule;
1147  vol_ele->getUserPolynomialBase() =
1148  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1149  vol_ele->getOpPtrVector().push_back(
1151  eP.dataAtPts->approxPAtPts));
1152  // vol_ele->getOpPtrVector().push_back(
1153  // new OpCalculateHVecTensorField<3, 3>(
1154  // eP.piolaStress + "_RT", eP.dataAtPts->approxPAtPts, MBMAXTYPE));
1155  vol_ele->getOpPtrVector().push_back(
1157  eP.bubbleField, eP.dataAtPts->approxPAtPts, MBMAXTYPE));
1158  vol_ele->getOpPtrVector().push_back(
1160  eP.streachTensor, eP.dataAtPts->streachTensorAtPts, MBTET));
1161  vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1162  eP.rotAxis, eP.dataAtPts->rotAxisAtPts, MBTET));
1163  vol_ele->getOpPtrVector().push_back(
1165  eP.materialGradient, eP.dataAtPts->GAtPts, MBTET));
1166  vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1167  eP.spatialDisp, eP.dataAtPts->wAtPts, MBTET));
1168  vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1169  eP.spatialPosition, eP.dataAtPts->xAtPts, MBVERTEX));
1170  vol_ele->getOpPtrVector().push_back(
1172  eP.dataAtPts->xGradAtPts));
1173  vol_ele->getOpPtrVector().push_back(
1175  eP.dataAtPts));
1176  };
1177 
1178 
1179  auto set_element_for_post_process = [&]() {
1180  volPostProcEle->getRuleHook = VolRule();
1181  volPostProcEle->getUserPolynomialBase() =
1182  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1183  volPostProcEle->getOpPtrVector().push_back(
1185  eP.dataAtPts->approxPAtPts));
1186  // volPostProcEle->getOpPtrVector().push_back(
1187  // new OpCalculateHVecTensorField<3, 3>(
1188  // eP.piolaStress + "_RT", eP.dataAtPts->approxPAtPts, MBMAXTYPE));
1189  volPostProcEle->getOpPtrVector().push_back(
1191  eP.bubbleField, eP.dataAtPts->approxPAtPts, MBMAXTYPE));
1192  volPostProcEle->getOpPtrVector().push_back(
1194  eP.streachTensor, eP.dataAtPts->streachTensorAtPts, MBTET));
1195  volPostProcEle->getOpPtrVector().push_back(
1197  eP.rotAxis, eP.dataAtPts->rotAxisAtPts, MBTET));
1198  volPostProcEle->getOpPtrVector().push_back(
1200  eP.materialGradient, eP.dataAtPts->GAtPts, MBTET));
1201  volPostProcEle->getOpPtrVector().push_back(
1203  eP.dataAtPts->wAtPts, MBTET));
1204  volPostProcEle->getOpPtrVector().push_back(
1206  eP.spatialPosition, eP.dataAtPts->xAtPts, MBVERTEX));
1207  volPostProcEle->getOpPtrVector().push_back(
1209  eP.dataAtPts->xGradAtPts));
1210  volPostProcEle->getOpPtrVector().push_back(
1212  eP.dataAtPts));
1213  volPostProcEle->getOpPtrVector().push_back(
1214  new OpCalculateStrainEnergy(eP.spatialDisp, eP.dataAtPts, gEnergy));
1215  };
1216 
1217  set_element_for_field_eval();
1218  set_element_for_post_process();
1219  }
1220 
1221  MoFEMErrorCode preProcess() { return 0; }
1222 
1223  MoFEMErrorCode operator()() { return 0; }
1224 
1225  MoFEMErrorCode postProcess() {
1227 
1229  1, "out_sol_elastic_" + boost::lexical_cast<std::string>(ts_step) +
1230  ".h5m");
1231 
1232  // Loop boundary elements with traction boundary conditions
1233  *gEnergy = 0;
1234  CHKERR eP.mField.loop_finite_elements(problemPtr->getName(), "EP",
1235  *volPostProcEle,nullptr);
1236 
1237  double body_energy;
1238  MPI_Allreduce(gEnergy.get(), &body_energy, 1, MPI_DOUBLE, MPI_SUM,
1239  eP.mField.get_comm());
1240  CHKERR PetscPrintf(eP.mField.get_comm(),
1241  "Step %d time %3.4g strain energy %3.6e\n", ts_step,
1242  ts_t, body_energy);
1243 
1244  auto post_proc_at_points = [&](std::array<double, 3> point,
1245  std::string str) {
1247 
1248  dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
1249 
1250  struct OpPrint : public VolOp {
1251 
1252  EshelbianCore &eP;
1253  std::array<double, 3> point;
1254  std::string str;
1255 
1256  OpPrint(EshelbianCore &ep, std::array<double, 3> &point,
1257  std::string &str)
1258  : VolOp(ep.spatialDisp, VolOp::OPROW), eP(ep), point(point),
1259  str(str) {}
1260 
1261  MoFEMErrorCode doWork(int side, EntityType type,
1264  if (type == MBTET) {
1265  if (getGaussPts().size2()) {
1266 
1267  auto t_h = getFTensor2FromMat<3, 3>(*(eP.dataAtPts->hAtPts));
1268  auto t_approx_P =
1269  getFTensor2FromMat<3, 3>(*(eP.dataAtPts->approxPAtPts));
1270 
1274  const double jac = dEterminant(t_h);
1276  t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
1277 
1278  auto add = [&]() {
1279  std::ostringstream s;
1280  s << str << " elem " << getFEEntityHandle() << " ";
1281  return s.str();
1282  };
1283 
1284  std::ostringstream print;
1285  print << add() << "comm rank " << eP.mField.get_comm_rank()
1286  << std::endl;
1287  print << add() << "point " << getVectorAdaptor(point.data(), 3)
1288  << std::endl;
1289  print << add() << "coords at gauss pts " << getCoordsAtGaussPts()
1290  << std::endl;
1291  print << add() << "w " << (*eP.dataAtPts->wAtPts) << std::endl;
1292  print << add() << "Piola " << (*eP.dataAtPts->approxPAtPts)
1293  << std::endl;
1294  print << add() << "Cauchy " << t_cauchy << std::endl;
1295  print << std::endl;
1296  CHKERR PetscSynchronizedPrintf(eP.mField.get_comm(), "%s",
1297  print.str().c_str());
1298 
1299  }
1300  }
1302  }
1303  };
1304 
1305  if (auto fe_ptr = dataFieldEval->feMethodPtr.lock()) {
1306 
1307  fe_ptr->getOpPtrVector().push_back(new OpPrint(eP, point, str));
1309  ->evalFEAtThePoint3D(point.data(), 1e-12, problemPtr->getName(),
1310  "EP", dataFieldEval,
1311  eP.mField.get_comm_rank(),
1313  fe_ptr->getOpPtrVector().pop_back();
1314  }
1315 
1317  };
1318 
1319  // Points for Cook beam
1320  std::array<double, 3> pointA = {48., 60., 5.};
1321  CHKERR post_proc_at_points(pointA, "Point A");
1322  PetscSynchronizedFlush(eP.mField.get_comm(), PETSC_STDOUT);
1323 
1324  std::array<double, 3> pointB = {48. / 2., 44. + (60. - 44.) / 2., 0.};
1325  CHKERR post_proc_at_points(pointB, "Point B");
1326  PetscSynchronizedFlush(eP.mField.get_comm(), PETSC_STDOUT);
1327 
1328  std::array<double, 3> pointC = {48. / 2., (44. - 0.) / 2., 0.};
1329  CHKERR post_proc_at_points(pointC, "Point C");
1330  PetscSynchronizedFlush(eP.mField.get_comm(), PETSC_STDOUT);
1331 
1333  }
1334  };
1335 
1336  boost::shared_ptr<FEMethod> monitor_ptr(new Monitor(*this));
1337  ts_ctx->get_loops_to_do_Monitor().push_back(
1339 
1340  CHKERR TSAppendOptionsPrefix(ts, "elastic_");
1341  CHKERR TSSetFromOptions(ts);
1343 }
1344 
1347 
1348  CHKERR TSSetDM(ts, dmElastic);
1349 
1350  SNES snes;
1351  CHKERR TSGetSNES(ts, &snes);
1352 
1353  PetscViewerAndFormat *vf;
1354  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
1355  PETSC_VIEWER_DEFAULT, &vf);
1356  CHKERR SNESMonitorSet(
1357  snes,
1358  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
1359  vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
1360 
1361  PetscSection section;
1362  CHKERR DMGetDefaultSection(dmElastic, &section);
1363  int num_fields;
1364  CHKERR PetscSectionGetNumFields(section, &num_fields);
1365  for (int ff = 0; ff != num_fields; ff++) {
1366  const char *field_name;
1367  CHKERR PetscSectionGetFieldName(section, ff, &field_name);
1368  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Field %d name %s\n", ff, field_name);
1369  }
1370 
1371  CHKERR DMoFEMMeshToLocalVector(dmElastic, x, INSERT_VALUES, SCATTER_FORWARD);
1372 
1373  // Adding field split solver
1374 
1375  KSP ksp;
1376  CHKERR SNESGetKSP(snes, &ksp);
1377  PC pc;
1378  CHKERR KSPGetPC(ksp, &pc);
1379  PetscBool is_uu_field_split;
1380  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_uu_field_split);
1381  if (is_uu_field_split) {
1382 
1383  const MoFEM::Problem *schur_uu_ptr;
1385  if (auto uu_data = schur_uu_ptr->subProblemData) {
1386 
1387  const MoFEM::Problem *prb_ptr;
1389  map<std::string, IS> is_map;
1390  for (int ff = 0; ff != num_fields; ff++) {
1391  const char *field_name;
1392  CHKERR PetscSectionGetFieldName(section, ff, &field_name);
1393  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1394  prb_ptr->getName(), ROW, field_name, 0, MAX_DOFS_ON_ENTITY,
1395  &is_map[field_name]);
1396  }
1397  // CHKERR mField.getInterface<ISManager>()
1398  // ->isCreateProblemFieldAndEntityType(
1399  // prb_ptr->getName(), ROW, piolaStress, MBTET, MBTET, 0,
1400  // MAX_DOFS_ON_ENTITY, &is_map["T_STRESS_ON_TETS"]);
1401  // CHKERR mField.getInterface<ISManager>()
1402  // ->isCreateProblemFieldAndEntityType(
1403  // prb_ptr->getName(), ROW, piolaStress, MBTRI, MBTRI, 0,
1404  // MAX_DOFS_ON_ENTITY, &is_map["T_STRESS_ON_TRIS"]);
1405 
1406  CHKERR uu_data->getRowIs(&is_map["E_IS_SUU"]);
1407 
1408  CHKERR PCFieldSplitSetIS(pc, NULL, is_map[streachTensor]);
1409  CHKERR PCFieldSplitSetIS(pc, NULL, is_map["E_IS_SUU"]);
1410 
1411  CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER,
1412  schurAssembly->Suu);
1413 
1414  CHKERR PCSetUp(pc);
1415  PetscInt n;
1416  KSP *uu_ksp;
1417  CHKERR PCFieldSplitGetSubKSP(pc, &n, &uu_ksp);
1418  PC bubble_pc;
1419  CHKERR KSPGetPC(uu_ksp[1], &bubble_pc);
1420  PetscBool is_bubble_field_split;
1421  PetscObjectTypeCompare((PetscObject)bubble_pc, PCFIELDSPLIT,
1422  &is_bubble_field_split);
1423  if (is_bubble_field_split) {
1424 
1425  const MoFEM::Problem *schur_bubble_ptr;
1426  CHKERR DMMoFEMGetProblemPtr(dmElasticSchurBubble, &schur_bubble_ptr);
1427  if (auto bubble_data = schur_bubble_ptr->subProblemData) {
1428 
1429  CHKERR bubble_data->getRowIs(&is_map["E_IS_BUBBLE"]);
1430 
1431  AO uu_ao;
1432  CHKERR uu_data->getRowMap(&uu_ao);
1433 
1434  CHKERR AOApplicationToPetscIS(uu_ao, is_map[bubbleField]);
1435  CHKERR PCFieldSplitSetIS(bubble_pc, NULL, is_map[bubbleField]);
1436  CHKERR PCFieldSplitSetIS(bubble_pc, NULL, is_map["E_IS_BUBBLE"]);
1437  CHKERR PCFieldSplitSetSchurPre(
1438  bubble_pc, PC_FIELDSPLIT_SCHUR_PRE_USER, schurAssembly->SBubble);
1439 
1440  CHKERR PCSetUp(bubble_pc);
1441  PetscInt bubble_n;
1442  KSP *bubble_ksp;
1443  CHKERR PCFieldSplitGetSubKSP(bubble_pc, &bubble_n, &bubble_ksp);
1444  PC omega_pc;
1445  CHKERR KSPGetPC(bubble_ksp[1], &omega_pc);
1446  PetscBool is_omega_field_split;
1447  PetscObjectTypeCompare((PetscObject)omega_pc, PCFIELDSPLIT,
1448  &is_omega_field_split);
1449 
1450  if (is_omega_field_split) {
1451 
1452  const MoFEM::Problem *schur_omega_ptr;
1453  CHKERR DMMoFEMGetProblemPtr(dmElasticSchurOmega, &schur_omega_ptr);
1454  if (auto omega_data = schur_omega_ptr->subProblemData) {
1455 
1456  AO bubble_ao;
1457  CHKERR bubble_data->getRowMap(&bubble_ao);
1458 
1459  CHKERR AOApplicationToPetscIS(uu_ao, is_map[rotAxis]);
1460  CHKERR AOApplicationToPetscIS(bubble_ao, is_map[rotAxis]);
1461  CHKERR omega_data->getRowIs(&is_map["E_IS_OMEGA"]);
1462 
1463  CHKERR PCFieldSplitSetIS(omega_pc, NULL, is_map[rotAxis]);
1464  CHKERR PCFieldSplitSetIS(omega_pc, NULL, is_map["E_IS_OMEGA"]);
1465  CHKERR PCFieldSplitSetSchurPre(omega_pc,
1466  PC_FIELDSPLIT_SCHUR_PRE_USER,
1467  schurAssembly->SOmega);
1468 
1469  CHKERR PCSetUp(omega_pc);
1470  PetscInt omega_n;
1471  KSP *omega_ksp;
1472  CHKERR PCFieldSplitGetSubKSP(omega_pc, &omega_n, &omega_ksp);
1473  PC w_pc;
1474  CHKERR KSPGetPC(omega_ksp[1], &w_pc);
1475  PetscBool is_w_field_split;
1476  PetscObjectTypeCompare((PetscObject)w_pc, PCFIELDSPLIT,
1477  &is_w_field_split);
1478  if (is_w_field_split) {
1479 
1480  const MoFEM::Problem *schur_w_ptr;
1482  &schur_w_ptr);
1483  if (auto w_data = schur_w_ptr->subProblemData) {
1484 
1485  AO omega_ao;
1486  CHKERR omega_data->getRowMap(&omega_ao);
1487 
1488  CHKERR AOApplicationToPetscIS(uu_ao, is_map[spatialDisp]);
1489  CHKERR AOApplicationToPetscIS(bubble_ao, is_map[spatialDisp]);
1490  CHKERR AOApplicationToPetscIS(omega_ao, is_map[spatialDisp]);
1491  CHKERR w_data->getRowIs(&is_map["E_IS_W"]);
1492 
1493  CHKERR PCFieldSplitSetIS(w_pc, NULL, is_map[spatialDisp]);
1494  CHKERR PCFieldSplitSetIS(w_pc, NULL, is_map["E_IS_W"]);
1495  CHKERR PCFieldSplitSetSchurPre(
1496  w_pc, PC_FIELDSPLIT_SCHUR_PRE_USER, schurAssembly->Sw);
1497 
1498  CHKERR AODestroy(&omega_ao);
1499  }
1500  }
1501 
1502  CHKERR PetscFree(omega_ksp);
1503  }
1504  }
1505  CHKERR PetscFree(bubble_ksp);
1506  CHKERR AODestroy(&uu_ao);
1507  }
1508  }
1509  CHKERR PetscFree(uu_ksp);
1510 
1511  for (auto &m : is_map)
1512  CHKERR ISDestroy(&m.second);
1513  }
1514  }
1515 
1516  CHKERR TSSolve(ts, x);
1518 }
1519 
1521  const std::string file) {
1523 
1524  if (!dataAtPts) {
1525  dataAtPts =
1526  boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
1527  dataAtPts->approxPAtPts = boost::make_shared<MatrixDouble>();
1528  dataAtPts->approxSigmaAtPts = boost::make_shared<MatrixDouble>();
1529  dataAtPts->divSigmaAtPts = boost::make_shared<MatrixDouble>();
1530  dataAtPts->wAtPts = boost::make_shared<MatrixDouble>();
1531  dataAtPts->xAtPts = boost::make_shared<MatrixDouble>();
1532  dataAtPts->xGradAtPts = boost::make_shared<MatrixDouble>();
1533  dataAtPts->xDetAtPts = boost::make_shared<VectorDouble>();
1534  dataAtPts->xInvGradAtPts = boost::make_shared<MatrixDouble>();
1535  dataAtPts->rotAxisAtPts = boost::make_shared<MatrixDouble>();
1536  dataAtPts->streachTensorAtPts = boost::make_shared<MatrixDouble>();
1537  dataAtPts->rotMatAtPts = boost::make_shared<MatrixDouble>();
1538  dataAtPts->diffRotMatAtPts = boost::make_shared<MatrixDouble>();
1539  dataAtPts->hAtPts = boost::make_shared<MatrixDouble>();
1540  dataAtPts->hDeltaAtPts = boost::make_shared<MatrixDouble>();
1541  dataAtPts->GAtPts = boost::make_shared<MatrixDouble>();
1542  dataAtPts->G0AtPts = boost::make_shared<MatrixDouble>();
1543  dataAtPts->hDeltaDetAtPts = boost::make_shared<VectorDouble>();
1544  }
1545 
1547  post_proc.getUserPolynomialBase() =
1548  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1549 
1551  post_proc.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
1552  piolaStress, dataAtPts->approxPAtPts));
1553  // post_proc.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
1554  // piolaStress + "_RT", dataAtPts->approxPAtPts, MBMAXTYPE));
1555  post_proc.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
1556  bubbleField, dataAtPts->approxPAtPts, MBMAXTYPE));
1557  post_proc.getOpPtrVector().push_back(
1559  streachTensor, dataAtPts->streachTensorAtPts, MBTET));
1560  post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1561  rotAxis, dataAtPts->rotAxisAtPts, MBTET));
1562  post_proc.getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
1563  materialGradient, dataAtPts->GAtPts, MBTET));
1564  post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1565  spatialDisp, dataAtPts->wAtPts, MBTET));
1566  post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1567  spatialPosition, dataAtPts->xAtPts, MBVERTEX));
1568  post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
1569  spatialPosition, dataAtPts->xGradAtPts));
1570 
1571  // evaluate derived quantities
1572  post_proc.getOpPtrVector().push_back(
1574 
1575  // evaluate integration points
1576  post_proc.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
1577  spatialDisp, tag, true, false, dataAtPts, physicalEquations));
1578  post_proc.getOpPtrVector().push_back(new OpPostProcDataStructure(
1579  post_proc.postProcMesh, post_proc.mapGaussPts, spatialDisp, dataAtPts));
1580 
1581  CHKERR DMoFEMLoopFiniteElements(dM, elementVolumeName.c_str(), &post_proc);
1582  CHKERR post_proc.writeFile(file.c_str());
1584 }
1585 
1586 } // namespace EshelbianPlasticity
Get field gradients time derivative at integration pts for scalar filed rank 0, i....
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 query_interface(const MOFEMuuid &uuid, BaseFunctionUnknownInterface **iface) const
structure for User Loop Methods on finite elementsIt can be used to calculate stiffness matrices,...
Get time direvatives of values at integration pts for tensor filed rank 1, i.e. vector field mofem_fo...
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:22
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
MoFEM interface unique ID.
Problem manager is used to build and partition problems.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEMErrorCode setBaseVolumeElementOps(const int tag, const bool do_rhs, const bool do_lhs, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore >> &fe)
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
field with continuous normal traction
Definition: definitions.h:173
user implemented approximation base
Definition: definitions.h:154
virtual moab::Interface & get_moab()=0
MoFEMErrorCode setElasticElementOps(const int tag)
boost::shared_ptr< EpFEMethod > schurAssembly
data structure for finite element entityIt keeps that about indices of degrees of freedom,...
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:318
boost::shared_ptr< BcRotVec > bcSpatialRotationVecPtr
Section manager is used to create indexes and sectionsFIXME: ISManager is not properly testsed by ato...
Definition: ISManager.hpp:36
moab::Interface & postProcMesh
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:434
boost::shared_ptr< SubProblemData > subProblemData
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problemIt if true is assumed that matrix has the same indexing on rows and columns....
Definition: DMMMoFEM.cpp:367
MoFEMErrorCode getValueHdivForCGGBubble(MatrixDouble &pts)
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
MoFEMErrorCode addDMs(const BitRefLevel bit=BitRefLevel().set(0))
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:414
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
auto & getUserPolynomialBase()
Get the User Polynomial Base object.
Class used to pass element data to calculate base functions on tet,triangle,edge.
const FieldApproximationBase bAse
base class for all interface classes
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:544
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:190
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
keeps basic data about problemThis is low level structure with information about problem,...
MoFEMErrorCode removeDofsOnEntities(const std::string problem_name, const std::string field_name, const Range ents, const int lo_coeff=0, const int hi_coeff=MAX_DOFS_ON_ENTITY, int verb=VERBOSE, const bool debug=false)
Remove DOFs from problemRemove DOFs from problem which are on entities on the given range and given f...
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition: DMMMoFEM.cpp:761
Get values at integration pts for tensor filed rank 1, i.e. vector field.
static const MOFEMuuid IDD_TET_BASE_FUNCTION
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:496
DM dmElasticSchurOmega
Sub problem of dmElastic Schur.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Projection of edge entities with one mid-node on hierarchical basis.
Field evaluator interface.
MoFEMErrorCode addFields(const EntityHandle meshset=0)
static const MOFEMuuid IDD_MOFEMEshelbianCrackInterface
Implementation of tonsorial bubble base div(v) = 0.
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMMoFEM.cpp:166
MoFEMErrorCode writeFile(const std::string file_name)
wrote results in (MOAB) format, use "file_name.h5m"
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
Definition: DMMMoFEM.cpp:187
int operator()(int p_row, int p_col, int p_data) const
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
virtual int get_comm_rank() const =0
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:331
boost::shared_ptr< BcDispVec > bcSpatialDispVecPtr
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
Base class if inherited used to calculate base functions.
MoFEMErrorCode postProcessResults(const int tag, const std::string file)
MoFEMErrorCode addBoundaryFiniteElement(const EntityHandle meshset=0)
MoFEMErrorCode addVolumeFiniteElement(const EntityHandle meshset=0)
static double dEterminant(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:385
DataForcesAndSourcesCore & dAta
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)
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
boost::shared_ptr< TractionFreeBc > bcSpatialFreeTraction
TractionBc(std::string name, std::vector< double > &attr, Range &faces)
MoFEMErrorCode setUpTSElastic(TS ts, Mat m, Vec f)
MoFEMErrorCode CGG_BubbleBase_MBTET(const int p, const double *N, const double *diffN, FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > &phi, const int gdim)
Calculate CGGT tonsorial bubble base.
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:146
BcRot(std::string name, std::vector< double > &attr, Range &faces)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:43
Post processing.
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)
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
boost::shared_ptr< PhysicalEquations > physicalEquations
MoFEMErrorCode solveElastic(TS ts, Vec x)
Calculate divergence of tonsorial field using vectorial base.
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMMoFEM.cpp:109
FieldOrderTable & getFieldOrderTable()
Get the Field Order Table.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
#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::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > elasticBcRhs
#define CHKERR
Inline error check.
Definition: definitions.h:596
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
DM dM
Coupled problem all fields.
std::vector< EntityHandle > mapGaussPts
ForcesAndSourcesCore::UserDataOperator UserDataOperator
DM dmElasticSchurStreach
Sub problem of dmElastic Schur.
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:285
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Set integration rule on element.
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMMoFEM.cpp:990
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > elasticFeLhs
DM dmElasticSchurSpatialDisp
Sub problem of dmElastic Schur.
Eshelbian plasticity interface.
boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > elasticBcLhs
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:298
virtual bool check_field(const std::string &name) const =0
check if field is in database
boost::shared_ptr< TractionBcVec > bcSpatialTraction
Data on single entity (This is passed as argument to DataOperator::doWork)
structure to get information form mofem into DataForcesAndSourcesCore
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMMoFEM.cpp:349
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
BcDisp(std::string name, std::vector< double > &attr, Range &faces)
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMMoFEM.cpp:708
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:337
virtual MPI_Comm & get_comm() const =0
continuous field
Definition: definitions.h:171
std::string getName() const
virtual const Field * get_field_structure(const std::string &name)=0
get field structure
DM dmElasticSchurBubble
Sub problem of dmElastic Schur.
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:971
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
boost::shared_ptr< SPD > getData(const double *ptr=nullptr, const int nb_eval_points=0, const double eps=1e-12, VERBOSITY_LEVELS verb=QUIET)
Get the Data object.
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > elasticFeRhs
DEPRECATED MoFEMErrorCode loop_finite_elements(const Problem *problem_ptr, const std::string &fe_name, FEMethod &method, int lower_rank, int upper_rank, MoFEMTypes bh, int verb=DEFAULT_VERBOSITY)
int operator()(int p_row, int p_col, int p_data) const
field with C-1 continuity
Definition: definitions.h:174