v0.14.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 #define SINGULARITY
9 #include <MoFEM.hpp>
10 using namespace MoFEM;
11 
12 #include <CGGTonsorialBubbleBase.hpp>
13 
14 #include <EshelbianPlasticity.hpp>
15 #include <boost/math/constants/constants.hpp>
16 
17 #include <cholesky.hpp>
18 #ifdef PYTHON_SDF
19  #include <boost/python.hpp>
20  #include <boost/python/def.hpp>
21  #include <boost/python/numpy.hpp>
22 namespace bp = boost::python;
23 namespace np = boost::python::numpy;
24 
25  #pragma message "With PYTHON_SDF"
26 
27 #else
28 
29  #pragma message "Without PYTHON_SDF"
30 
31 #endif
32 
33 #include <EshelbianContact.hpp>
34 
35 extern "C" {
36 #include <phg-quadrule/quad.h>
37 }
38 
39 
40 namespace EshelbianPlasticity {
41 struct VolUserDataOperatorStabAssembly
42  : public VolumeElementForcesAndSourcesCore::UserDataOperator {
43  using VolumeElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
44 };
45 struct FaceUserDataOperatorStabAssembly
46  : public FaceElementForcesAndSourcesCore::UserDataOperator {
47  using FaceElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
48 };
49 
50 } // namespace EshelbianPlasticity
51 
52 static auto get_range_from_block(MoFEM::Interface &m_field,
53  const std::string block_name, int dim) {
54  Range r;
55 
56  auto mesh_mng = m_field.getInterface<MeshsetsManager>();
57  auto bcs = mesh_mng->getCubitMeshsetPtr(
58 
59  std::regex((boost::format("%s(.*)") % block_name).str())
60 
61  );
62 
63  for (auto bc : bcs) {
64  Range faces;
65  CHK_MOAB_THROW(bc->getMeshsetIdEntitiesByDimension(m_field.get_moab(), dim,
66  faces, true),
67  "get meshset ents");
68  r.merge(faces);
69  }
70 
71  return r;
72 };
73 
74 static auto save_range(moab::Interface &moab, const std::string name,
75  const Range r) {
76  MoFEMFunctionBegin;
77  auto out_meshset = get_temp_meshset_ptr(moab);
78  CHKERR moab.add_entities(*out_meshset, r);
79  if (r.size()) {
80  CHKERR moab.write_file(name.c_str(), "VTK", "", out_meshset->get_ptr(), 1);
81  } else {
82  MOFEM_LOG("SELF", Sev::warning) << "Empty range for " << name;
83  }
84  MoFEMFunctionReturn(0);
85 };
86 
87 template <>
88 typename MoFEM::OpBaseImpl<
89  SCHUR,
90  EshelbianPlasticity::VolUserDataOperatorStabAssembly>::MatSetValuesHook
91  MoFEM::OpBaseImpl<SCHUR,
92  EshelbianPlasticity::VolUserDataOperatorStabAssembly>::
93  matSetValuesHook = [](ForcesAndSourcesCore::UserDataOperator *op_ptr,
94  const EntitiesFieldData::EntData &row_data,
95  const EntitiesFieldData::EntData &col_data,
96  MatrixDouble &m) {
97  return MatSetValues<AssemblyTypeSelector<SCHUR>>(
98  op_ptr->getKSPB(), row_data, col_data, m, ADD_VALUES);
99  };
100 
101 namespace EshelbianPlasticity {
102 
103 struct SetIntegrationAtFrontVolume {
104 
105  SetIntegrationAtFrontVolume() = default;
106 
107  MoFEMErrorCode operator()(ForcesAndSourcesCore *fe_raw_ptr, int order_row,
108  int order_col, int order_data) {
109  MoFEMFunctionBegin;
110 
111  auto &m_field = fe_raw_ptr->mField;
112  auto vol_fe_ptr = dynamic_cast<VolFe *>(fe_raw_ptr);
113  auto fe_handle = vol_fe_ptr->getFEEntityHandle();
114 
115  auto set_base_quadrature = [&]() {
116  MoFEMFunctionBegin;
117  int rule = 2 * order_data + 1;
118  if (rule < QUAD_3D_TABLE_SIZE) {
119  if (QUAD_3D_TABLE[rule]->dim != 3) {
120  SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
121  "wrong dimension");
122  }
123  if (QUAD_3D_TABLE[rule]->order < rule) {
124  SETERRQ2(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
125  "wrong order %d != %d", QUAD_3D_TABLE[rule]->order, rule);
126  }
127  const size_t nb_gauss_pts = QUAD_3D_TABLE[rule]->npoints;
128  auto &gauss_pts = vol_fe_ptr->gaussPts;
129  gauss_pts.resize(4, nb_gauss_pts, false);
130  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[1], 4,
131  &gauss_pts(0, 0), 1);
132  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[2], 4,
133  &gauss_pts(1, 0), 1);
134  cblas_dcopy(nb_gauss_pts, &QUAD_3D_TABLE[rule]->points[3], 4,
135  &gauss_pts(2, 0), 1);
136  cblas_dcopy(nb_gauss_pts, QUAD_3D_TABLE[rule]->weights, 1,
137  &gauss_pts(3, 0), 1);
138  auto &data = vol_fe_ptr->dataOnElement[H1];
139  data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4,
140  false);
141  double *shape_ptr =
142  &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
143  cblas_dcopy(4 * nb_gauss_pts, QUAD_3D_TABLE[rule]->points, 1,
144  shape_ptr, 1);
145  } else {
146  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
147  "rule > quadrature order %d < %d", rule,
148  QUAD_3D_TABLE_SIZE);
149  }
150  MoFEMFunctionReturn(0);
151  };
152 
153  CHKERR set_base_quadrature();
154 
155  if (frontNodes->size() && EshelbianCore::setSingularity) {
156 
157  auto get_singular_nodes = [&]() {
158  int num_nodes;
159  const EntityHandle *conn;
160  CHKERR m_field.get_moab().get_connectivity(fe_handle, conn, num_nodes,
161  true);
162  std::bitset<numNodes> singular_nodes;
163  for (auto nn = 0; nn != numNodes; ++nn) {
164  if (frontNodes->find(conn[nn]) != frontNodes->end()) {
165  singular_nodes.set(nn);
166  } else {
167  singular_nodes.reset(nn);
168  }
169  }
170  return singular_nodes;
171  };
172 
173  auto get_singular_edges = [&]() {
174  std::bitset<numEdges> singular_edges;
175  for (int ee = 0; ee != numEdges; ee++) {
176  EntityHandle edge;
177  CHKERR m_field.get_moab().side_element(fe_handle, 1, ee, edge);
178  if (frontEdges->find(edge) != frontEdges->end()) {
179  singular_edges.set(ee);
180  } else {
181  singular_edges.reset(ee);
182  }
183  }
184  return singular_edges;
185  };
186 
187  auto set_gauss_pts = [&](auto &ref_gauss_pts) {
188  MoFEMFunctionBegin;
189  vol_fe_ptr->gaussPts.swap(ref_gauss_pts);
190  const size_t nb_gauss_pts = vol_fe_ptr->gaussPts.size2();
191  // auto &dataH1 = vol_fe_ptr->dataOnElement[H1];
192  // data->dataOnEntities[MBVERTEX][0].getN(NOBASE).resize(nb_gauss_pts, 4);
193  // double *shape_ptr =
194  // &*data->dataOnEntities[MBVERTEX][0].getN(NOBASE).data().begin();
195  // CHKERR ShapeMBTET(shape_ptr, &vol_fe_ptr->gaussPts(0, 0),
196  // &vol_fe_ptr->gaussPts(1, 0), &gaussPts(2, 0),
197  // nb_gauss_pts);
198  MoFEMFunctionReturn(0);
199  };
200 
201  auto singular_nodes = get_singular_nodes();
202  if (singular_nodes.count()) {
203  auto it_map_ref_coords = mapRefCoords.find(singular_nodes.to_ulong());
204  if (it_map_ref_coords != mapRefCoords.end()) {
205  CHKERR set_gauss_pts(it_map_ref_coords->second);
206  MoFEMFunctionReturnHot(0);
207  } else {
208 
209  auto refine_quadrature = [&]() {
210  MoFEMFunctionBegin;
211 
212  const int max_level = refinementLevels;
213  EntityHandle tet;
214 
215  moab::Core moab_ref;
216  double base_coords[] = {0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1};
217  EntityHandle nodes[4];
218  for (int nn = 0; nn != 4; nn++)
219  CHKERR moab_ref.create_vertex(&base_coords[3 * nn], nodes[nn]);
220  CHKERR moab_ref.create_element(MBTET, nodes, 4, tet);
221  MoFEM::CoreTmp<-1> mofem_ref_core(moab_ref, PETSC_COMM_SELF, -2);
222  MoFEM::Interface &m_field_ref = mofem_ref_core;
223  CHKERR m_field_ref.getInterface<BitRefManager>()
224  ->setBitRefLevelByDim(0, 3, BitRefLevel().set(0));
225 
226  Range nodes_at_front;
227  for (int nn = 0; nn != numNodes; nn++) {
228  if (singular_nodes[nn]) {
229  EntityHandle ent;
230  CHKERR moab_ref.side_element(tet, 0, nn, ent);
231  nodes_at_front.insert(ent);
232  }
233  }
234 
235  auto singular_edges = get_singular_edges();
236 
237  EntityHandle meshset;
238  CHKERR moab_ref.create_meshset(MESHSET_SET, meshset);
239  for (int ee = 0; ee != 6; ee++) {
240  if (singular_edges[ee]) {
241  EntityHandle ent;
242  CHKERR moab_ref.side_element(tet, 1, ee, ent);
243  CHKERR moab_ref.add_entities(meshset, &ent, 1);
244  }
245  }
246 
247  // refine mesh
248  auto *m_ref = m_field_ref.getInterface<MeshRefinement>();
249  for (int ll = 0; ll != max_level; ll++) {
250  Range edges;
251  CHKERR m_field_ref.getInterface<BitRefManager>()
252  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(ll),
253  BitRefLevel().set(), MBEDGE,
254  edges);
255  Range ref_edges;
256  CHKERR moab_ref.get_adjacencies(
257  nodes_at_front, 1, true, ref_edges, moab::Interface::UNION);
258  ref_edges = intersect(ref_edges, edges);
259  Range ents;
260  CHKERR moab_ref.get_entities_by_type(meshset, MBEDGE, ents, true);
261  ref_edges = intersect(ref_edges, ents);
262  Range tets;
263  CHKERR m_field_ref.getInterface<BitRefManager>()
264  ->getEntitiesByTypeAndRefLevel(
265  BitRefLevel().set(ll), BitRefLevel().set(), MBTET, tets);
266  CHKERR m_ref->addVerticesInTheMiddleOfEdges(
267  ref_edges, BitRefLevel().set(ll + 1));
268  CHKERR m_ref->refineTets(tets, BitRefLevel().set(ll + 1));
269  CHKERR m_field_ref.getInterface<BitRefManager>()
270  ->updateMeshsetByEntitiesChildren(meshset,
271  BitRefLevel().set(ll + 1),
272  meshset, MBEDGE, true);
273  }
274 
275  // get ref coords
276  Range tets;
277  CHKERR m_field_ref.getInterface<BitRefManager>()
278  ->getEntitiesByTypeAndRefLevel(BitRefLevel().set(max_level),
279  BitRefLevel().set(), MBTET,
280  tets);
281  MatrixDouble ref_coords(tets.size(), 12, false);
282  int tt = 0;
283  for (Range::iterator tit = tets.begin(); tit != tets.end();
284  tit++, tt++) {
285  int num_nodes;
286  const EntityHandle *conn;
287  CHKERR moab_ref.get_connectivity(*tit, conn, num_nodes, false);
288  CHKERR moab_ref.get_coords(conn, num_nodes, &ref_coords(tt, 0));
289  }
290 
291  auto &data = vol_fe_ptr->dataOnElement[H1];
292  const size_t nb_gauss_pts = vol_fe_ptr->gaussPts.size2();
293  MatrixDouble ref_gauss_pts(4, nb_gauss_pts * ref_coords.size1());
294  MatrixDouble &shape_n =
295  data->dataOnEntities[MBVERTEX][0].getN(NOBASE);
296  int gg = 0;
297  for (size_t tt = 0; tt != ref_coords.size1(); tt++) {
298  double *tet_coords = &ref_coords(tt, 0);
299  double det = Tools::tetVolume(tet_coords);
300  det *= 6;
301  for (size_t ggg = 0; ggg != nb_gauss_pts; ++ggg, ++gg) {
302  for (int dd = 0; dd != 3; dd++) {
303  ref_gauss_pts(dd, gg) =
304  shape_n(ggg, 0) * tet_coords[3 * 0 + dd] +
305  shape_n(ggg, 1) * tet_coords[3 * 1 + dd] +
306  shape_n(ggg, 2) * tet_coords[3 * 2 + dd] +
307  shape_n(ggg, 3) * tet_coords[3 * 3 + dd];
308  }
309  ref_gauss_pts(3, gg) = vol_fe_ptr->gaussPts(3, ggg) * det;
310  }
311  }
312 
313  mapRefCoords[singular_nodes.to_ulong()].swap(ref_gauss_pts);
314  CHKERR set_gauss_pts(mapRefCoords[singular_nodes.to_ulong()]);
315 
316  MoFEMFunctionReturn(0);
317  };
318  }
319  }
320  }
321 
322  MoFEMFunctionReturn(0);
323  }
324 
325  struct VolFe : public VolumeElementForcesAndSourcesCore {
326  using VolumeElementForcesAndSourcesCore::dataOnElement;
327  using VolumeElementForcesAndSourcesCore::VolumeElementForcesAndSourcesCore;
328  };
329 
330  static inline constexpr int numNodes = 4;
331  static inline constexpr int numEdges = 6;
332  static inline constexpr int refinementLevels = 3;
333 
334  static inline boost::shared_ptr<Range> frontNodes;
335  static inline boost::shared_ptr<Range> frontEdges;
336  static inline std::map<long int, MatrixDouble> mapRefCoords;
337 };
338 
339 double EshelbianCore::exponentBase = exp(1);
340 boost::function<double(const double)> EshelbianCore::f = EshelbianCore::f_log_e;
341 boost::function<double(const double)> EshelbianCore::d_f =
342  EshelbianCore::d_f_log_e;
343 boost::function<double(const double)> EshelbianCore::dd_f =
344  EshelbianCore::dd_f_log_e;
345 boost::function<double(const double)> EshelbianCore::inv_f =
346  EshelbianCore::inv_f_log;
347 boost::function<double(const double)> EshelbianCore::inv_d_f =
348  EshelbianCore::inv_d_f_log;
349 boost::function<double(const double)> EshelbianCore::inv_dd_f =
350  EshelbianCore::inv_dd_f_log;
351 
352 MoFEMErrorCode
353 EshelbianCore::query_interface(boost::typeindex::type_index type_index,
354  UnknownInterface **iface) const {
355  *iface = const_cast<EshelbianCore *>(this);
356  return 0;
357 }
358 
359 MoFEMErrorCode OpJacobian::doWork(int side, EntityType type, EntData &data) {
360  MoFEMFunctionBegin;
361 
362  if (evalRhs)
363  CHKERR evaluateRhs(data);
364 
365  if (evalLhs)
366  CHKERR evaluateLhs(data);
367 
368  MoFEMFunctionReturn(0);
369 }
370 
371 EshelbianCore::EshelbianCore(MoFEM::Interface &m_field) : mField(m_field) {
372  CHK_THROW_MESSAGE(getOptions(), "getOptions failed");
373 }
374 
375 EshelbianCore::~EshelbianCore() = default;
376 
377 MoFEMErrorCode EshelbianCore::getOptions() {
378  MoFEMFunctionBegin;
379  const char *list_rots[] = {"small", "moderate", "large", "no_h1"};
380  PetscInt choice_rot = EshelbianCore::rotSelector;
381  PetscInt choice_grad = EshelbianCore::gradApproximator;
382 
383  const char *list_stretches[] = {"linear", "log"};
384  PetscInt choice_stretch = StretchSelector::LOG;
385 
386  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
387  "none");
388  CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
389  spaceOrder, &spaceOrder, PETSC_NULL);
390  CHKERR PetscOptionsInt("-space_h1_order", "approximation oder for space", "",
391  spaceH1Order, &spaceH1Order, PETSC_NULL);
392  CHKERR PetscOptionsInt("-material_order", "approximation oder for material",
393  "", materialOrder, &materialOrder, PETSC_NULL);
394  CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alphaU,
395  &alphaU, PETSC_NULL);
396  CHKERR PetscOptionsScalar("-viscosity_alpha_w", "viscosity", "", alphaW,
397  &alphaW, PETSC_NULL);
398  CHKERR PetscOptionsScalar("-density_alpha_rho", "density", "", alphaRho,
399  &alphaRho, PETSC_NULL);
400  CHKERR PetscOptionsEList("-rotations", "rotations", "", list_rots,
401  LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
402  PETSC_NULL);
403  CHKERR PetscOptionsEList("-grad", "gradient of defamation approximate", "",
404  list_rots, NO_H1_CONFIGURATION + 1,
405  list_rots[choice_grad], &choice_grad, PETSC_NULL);
406  CHKERR PetscOptionsScalar("-exponent_base", "exponent_base", "", exponentBase,
407  &EshelbianCore::exponentBase, PETSC_NULL);
408  CHKERR PetscOptionsEList(
409  "-stretches", "stretches", "", list_stretches, StretchSelector::LOG + 1,
410  list_stretches[choice_stretch], &choice_stretch, PETSC_NULL);
411 
412  CHKERR PetscOptionsBool("-no_stretch", "do not solve for stretch", "",
413  noStretch, &noStretch, PETSC_NULL);
414  CHKERR PetscOptionsBool("-set_singularity", "set singularity", "",
415  setSingularity, &setSingularity, PETSC_NULL);
416 
417  // contact parameters
418  CHKERR PetscOptionsInt("-contact_max_post_proc_ref_level", "refinement level",
419  "", contactRefinementLevels, &contactRefinementLevels,
420  PETSC_NULL);
421 
422  ierr = PetscOptionsEnd();
423  CHKERRG(ierr);
424 
425  EshelbianCore::rotSelector = static_cast<RotSelector>(choice_rot);
426  EshelbianCore::gradApproximator = static_cast<RotSelector>(choice_grad);
427  EshelbianCore::stretchSelector = static_cast<StretchSelector>(choice_stretch);
428 
429  switch (EshelbianCore::stretchSelector) {
430  case StretchSelector::LINEAR:
431  EshelbianCore::f = f_linear;
432  EshelbianCore::d_f = d_f_linear;
433  EshelbianCore::dd_f = dd_f_linear;
434  EshelbianCore::inv_f = inv_f_linear;
435  EshelbianCore::inv_d_f = inv_d_f_linear;
436  EshelbianCore::inv_dd_f = inv_dd_f_linear;
437  break;
438  case StretchSelector::LOG:
439  if (EshelbianCore::exponentBase != exp(1)) {
440  EshelbianCore::f = f_log;
441  EshelbianCore::d_f = d_f_log;
442  EshelbianCore::dd_f = dd_f_log;
443  EshelbianCore::inv_f = inv_f_log;
444  EshelbianCore::inv_d_f = inv_d_f_log;
445  EshelbianCore::inv_dd_f = inv_dd_f_log;
446  } else {
447  EshelbianCore::f = f_log_e;
448  EshelbianCore::d_f = d_f_log_e;
449  EshelbianCore::dd_f = dd_f_log_e;
450  EshelbianCore::inv_f = inv_f_log;
451  EshelbianCore::inv_d_f = inv_d_f_log;
452  EshelbianCore::inv_dd_f = inv_dd_f_log;
453  }
454  break;
455  default:
456  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown stretch");
457  break;
458  };
459 
460  MOFEM_LOG("EP", Sev::inform) << "spaceOrder " << spaceOrder;
461  MOFEM_LOG("EP", Sev::inform) << "spaceH1Order " << spaceH1Order;
462  MOFEM_LOG("EP", Sev::inform) << "materialOrder " << materialOrder;
463  MOFEM_LOG("EP", Sev::inform) << "alphaU " << alphaU;
464  MOFEM_LOG("EP", Sev::inform) << "alphaW " << alphaW;
465  MOFEM_LOG("EP", Sev::inform) << "alphaRho " << alphaRho;
466  MOFEM_LOG("EP", Sev::inform)
467  << "Rotations " << list_rots[EshelbianCore::rotSelector];
468  MOFEM_LOG("EP", Sev::inform) << "Gradient of deformation "
469  << list_rots[EshelbianCore::gradApproximator];
470  if (exponentBase != exp(1))
471  MOFEM_LOG("EP", Sev::inform)
472  << "Base exponent " << EshelbianCore::exponentBase;
473  else
474  MOFEM_LOG("EP", Sev::inform) << "Base exponent e";
475  MOFEM_LOG("EP", Sev::inform) << "Stretch " << list_stretches[choice_stretch];
476 
477  if (spaceH1Order == -1)
478  spaceH1Order = spaceOrder;
479 
480  MoFEMFunctionReturn(0);
481 }
482 
483 MoFEMErrorCode EshelbianCore::addFields(const EntityHandle meshset) {
484  MoFEMFunctionBegin;
485 
486  auto get_tets = [&]() {
487  Range tets;
488  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
489  return tets;
490  };
491 
492  auto get_tets_skin = [&]() {
493  Range tets_skin_part;
494  Skinner skin(&mField.get_moab());
495  CHKERR skin.find_skin(0, get_tets(), false, tets_skin_part);
496  ParallelComm *pcomm =
497  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
498  Range tets_skin;
499  CHKERR pcomm->filter_pstatus(tets_skin_part,
500  PSTATUS_SHARED | PSTATUS_MULTISHARED,
501  PSTATUS_NOT, -1, &tets_skin);
502  return tets_skin;
503  };
504 
505  auto subtract_boundary_conditions = [&](auto &&tets_skin) {
506  // That mean, that hybrid field on all faces on which traction is applied,
507  // on other faces, or enforcing displacements as
508  // natural boundary condition.
509  if (bcSpatialTraction)
510  for (auto &v : *bcSpatialTraction) {
511  tets_skin = subtract(tets_skin, v.faces);
512  }
513  return tets_skin;
514  };
515 
516  auto add_blockset = [&](auto block_name, auto &&tets_skin) {
517  auto crack_faces =
518  get_range_from_block(mField, "block_name", SPACE_DIM - 1);
519  tets_skin.merge(crack_faces);
520  return tets_skin;
521  };
522 
523  auto subtract_blockset = [&](auto block_name, auto &&tets_skin) {
524  auto contact_range =
525  get_range_from_block(mField, block_name, SPACE_DIM - 1);
526  tets_skin = subtract(tets_skin, contact_range);
527  return tets_skin;
528  };
529 
530  auto get_stress_trace_faces = [&](auto &&tets_skin) {
531  Range faces;
532  CHKERR mField.get_moab().get_adjacencies(get_tets(), SPACE_DIM - 1, true,
533  faces, moab::Interface::UNION);
534  Range trace_faces = subtract(faces, tets_skin);
535  return trace_faces;
536  };
537 
538  auto tets = get_tets();
539 
540  // remove also contact faces, i.e. that is also kind of hybrid field but
541  // named but used to enforce contact conditions
542  auto trace_faces = get_stress_trace_faces(
543 
544  subtract_blockset("CONTACT",
545  subtract_boundary_conditions(get_tets_skin()))
546 
547  );
548 
549  contactFaces = boost::make_shared<Range>(intersect(
550  trace_faces, get_range_from_block(mField, "CONTACT", SPACE_DIM - 1)));
551  skeletonFaces =
552  boost::make_shared<Range>(subtract(trace_faces, *contactFaces));
553  materialSkeletonFaces =
554  boost::make_shared<Range>(get_stress_trace_faces(Range()));
555 
556 #ifndef NDEBUG
557  if (contactFaces->size())
558  CHKERR save_range(mField.get_moab(), "contact_faces.vtk", *contactFaces);
559  if (skeletonFaces->size())
560  CHKERR save_range(mField.get_moab(), "skeleton_faces.vtk", *skeletonFaces);
561  if (skeletonFaces->size())
562  CHKERR save_range(mField.get_moab(), "material_skeleton_faces.vtk",
563  *materialSkeletonFaces);
564 #endif
565 
566  auto add_broken_hdiv_field = [this, meshset](const std::string field_name,
567  const int order) {
568  MoFEMFunctionBegin;
569 
570  constexpr FieldApproximationBase base = DEMKOWICZ_JACOBI_BASE;
571 
572  auto get_side_map_hdiv = [&]() {
573  return std::vector<
574 
575  std::pair<EntityType,
576  std::function<MoFEMErrorCode(BaseFunction::DofsSideMap &)>
577 
578  >>{
579 
580  {MBTET,
581  [&](BaseFunction::DofsSideMap &dofs_side_map) -> MoFEMErrorCode {
582  return TetPolynomialBase::setDofsSideMap(HDIV, DISCONTINUOUS, base,
583  dofs_side_map);
584  }}
585 
586  };
587  };
588 
589  CHKERR mField.add_broken_field(field_name, HDIV, base, SPACE_DIM,
590  get_side_map_hdiv(), MB_TAG_DENSE, MF_ZERO);
591  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
592  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
593  MoFEMFunctionReturn(0);
594  };
595 
596  auto add_l2_field = [this, meshset](const std::string field_name,
597  const int order, const int dim) {
598  MoFEMFunctionBegin;
599  CHKERR mField.add_field(field_name, L2, AINSWORTH_LEGENDRE_BASE, dim,
600  MB_TAG_DENSE, MF_ZERO);
601  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
602  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
603  MoFEMFunctionReturn(0);
604  };
605 
606  auto add_h1_field = [this, meshset](const std::string field_name,
607  const int order, const int dim) {
608  MoFEMFunctionBegin;
609  CHKERR mField.add_field(field_name, H1, AINSWORTH_LEGENDRE_BASE, dim,
610  MB_TAG_DENSE, MF_ZERO);
611  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
612  CHKERR mField.set_field_order(meshset, MBVERTEX, field_name, 1);
613  CHKERR mField.set_field_order(meshset, MBEDGE, field_name, order);
614  CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
615  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
616  MoFEMFunctionReturn(0);
617  };
618 
619  auto add_l2_field_by_range = [this](const std::string field_name,
620  const int order, const int dim,
621  const int field_dim, Range &&r) {
622  MoFEMFunctionBegin;
623  CHKERR mField.add_field(field_name, L2, AINSWORTH_LEGENDRE_BASE, field_dim,
624  MB_TAG_DENSE, MF_ZERO);
625  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(r);
626  CHKERR mField.add_ents_to_field_by_dim(r, dim, field_name);
627  CHKERR mField.set_field_order(r, field_name, order);
628  MoFEMFunctionReturn(0);
629  };
630 
631  auto add_bubble_field = [this, meshset](const std::string field_name,
632  const int order, const int dim) {
633  MoFEMFunctionBegin;
634  CHKERR mField.add_field(field_name, HDIV, USER_BASE, dim, MB_TAG_DENSE,
635  MF_ZERO);
636  // Modify field
637  auto field_ptr = mField.get_field_structure(field_name);
638  auto field_order_table =
639  const_cast<Field *>(field_ptr)->getFieldOrderTable();
640  auto get_cgg_bubble_order_zero = [](int p) { return 0; };
641  auto get_cgg_bubble_order_tet = [](int p) {
642  return NBVOLUMETET_CCG_BUBBLE(p);
643  };
644  field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
645  field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
646  field_order_table[MBTRI] = get_cgg_bubble_order_zero;
647  field_order_table[MBTET] = get_cgg_bubble_order_tet;
648  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
649  CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
650  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
651  MoFEMFunctionReturn(0);
652  };
653 
654  auto add_user_l2_field = [this, meshset](const std::string field_name,
655  const int order, const int dim) {
656  MoFEMFunctionBegin;
657  CHKERR mField.add_field(field_name, L2, USER_BASE, dim, MB_TAG_DENSE,
658  MF_ZERO);
659  // Modify field
660  auto field_ptr = mField.get_field_structure(field_name);
661  auto field_order_table =
662  const_cast<Field *>(field_ptr)->getFieldOrderTable();
663  auto zero_dofs = [](int p) { return 0; };
664  auto dof_l2_tet = [](int p) { return NBVOLUMETET_L2(p); };
665  field_order_table[MBVERTEX] = zero_dofs;
666  field_order_table[MBEDGE] = zero_dofs;
667  field_order_table[MBTRI] = zero_dofs;
668  field_order_table[MBTET] = dof_l2_tet;
669  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
670  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
671  MoFEMFunctionReturn(0);
672  };
673 
674  // spatial fields
675  CHKERR add_broken_hdiv_field(piolaStress, spaceOrder);
676  CHKERR add_bubble_field(bubbleField, spaceOrder, 1);
677  CHKERR add_l2_field(spatialL2Disp, spaceOrder - 1, 3);
678  CHKERR add_user_l2_field(rotAxis, spaceOrder - 1, 3);
679  CHKERR add_user_l2_field(stretchTensor, noStretch ? -1 : spaceOrder, 6);
680 
681  if (!skeletonFaces)
682  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No skeleton faces");
683  if (!contactFaces)
684  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No contact faces");
685 
686  CHKERR add_l2_field_by_range(hybridSpatialDisp, spaceOrder - 1, 2, 3,
687  Range(*skeletonFaces));
688  CHKERR add_l2_field_by_range(contactDisp, spaceOrder - 1, 2, 3,
689  Range(*contactFaces));
690 
691  // spatial displacement
692  CHKERR add_h1_field(spatialH1Disp, spaceH1Order, 3);
693  // material positions
694  CHKERR add_h1_field(materialH1Positions, 2, 3);
695 
696  // Eshelby stress
697  CHKERR add_broken_hdiv_field(eshelbyStress, spaceOrder);
698  CHKERR add_l2_field(materialL2Disp, spaceOrder - 1, 3);
699  CHKERR add_l2_field_by_range(hybridMaterialDisp, spaceOrder - 1, 2, 3,
700  Range(*materialSkeletonFaces));
701 
702  CHKERR mField.build_fields();
703 
704  MoFEMFunctionReturn(0);
705 }
706 
707 MoFEMErrorCode EshelbianCore::projectGeometry(const EntityHandle meshset) {
708  MoFEMFunctionBegin;
709 
710  auto project_ho_geometry = [&]() {
711  Projection10NodeCoordsOnField ent_method(mField, materialH1Positions);
712  return mField.loop_dofs(materialH1Positions, ent_method);
713  };
714  CHKERR project_ho_geometry();
715 
716  auto get_skin = [&](auto &body_ents) {
717  Skinner skin(&mField.get_moab());
718  Range skin_ents;
719  CHKERR skin.find_skin(0, body_ents, false, skin_ents);
720  return skin_ents;
721  };
722 
723  auto filter_true_skin = [&](auto &&skin) {
724  Range boundary_ents;
725  ParallelComm *pcomm =
726  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
727  CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
728  PSTATUS_NOT, -1, &boundary_ents);
729  return boundary_ents;
730  };
731 
732  auto get_crack_front_edges = [&]() {
733  auto &moab = mField.get_moab();
734  auto crack_skin = get_skin(*crackFaces);
735  Range body_ents;
736  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
737  auto body_skin = filter_true_skin(get_skin(body_ents));
738  Range body_skin_edges;
739  CHKERR moab.get_adjacencies(body_skin, SPACE_DIM - 2, true, body_skin_edges,
740  moab::Interface::UNION);
741 
742  crack_skin = subtract(crack_skin, body_skin_edges);
743  std::map<int, Range> received_ents;
744  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
745  crack_skin, &received_ents);
746 
747  Range to_remove;
748  for (auto &m : received_ents) {
749  if (m.first != mField.get_comm_rank()) {
750  // if the same edges appear on two or more processors, they are not at
751  // the front
752  to_remove.merge(intersect(crack_skin, m.second));
753  }
754  }
755  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(to_remove);
756  crack_skin = subtract(crack_skin, to_remove);
757  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
758  crack_skin);
759 
760  return boost::make_shared<Range>(crack_skin);
761  };
762 
763  auto get_adj_front_edges = [&](auto &front_edges) {
764  auto &moab = mField.get_moab();
765  Range front_crack_nodes;
766  CHKERR moab.get_connectivity(front_edges, front_crack_nodes, true);
767  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
768  front_crack_nodes);
769  Range crack_front_edges;
770  CHKERR moab.get_adjacencies(front_crack_nodes, SPACE_DIM - 2, false,
771  crack_front_edges, moab::Interface::UNION);
772  crack_front_edges = subtract(crack_front_edges, front_edges);
773  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
774  crack_front_edges);
775 
776  Range crack_front_edges_nodes;
777  CHKERR moab.get_connectivity(crack_front_edges, crack_front_edges_nodes,
778  true);
779  crack_front_edges_nodes =
780  subtract(crack_front_edges_nodes, front_crack_nodes);
781  Range crack_front_edges_with_both_nodes_not_at_front;
782  CHKERR moab.get_adjacencies(crack_front_edges_nodes, SPACE_DIM - 2, false,
783  crack_front_edges_with_both_nodes_not_at_front,
784  moab::Interface::UNION);
785  crack_front_edges_with_both_nodes_not_at_front = intersect(
786  crack_front_edges_with_both_nodes_not_at_front, crack_front_edges);
787 
788  return std::make_pair(boost::make_shared<Range>(front_crack_nodes),
789  boost::make_shared<Range>(
790  crack_front_edges_with_both_nodes_not_at_front));
791  };
792 
793  crackFaces = boost::make_shared<Range>(
794  get_range_from_block(mField, "CRACK", SPACE_DIM - 1));
795  frontEdges = get_crack_front_edges();
796  auto [front_vertices, front_adj_edges] = get_adj_front_edges(*frontEdges);
797  frontVertices = front_vertices;
798  frontAdjEdges = front_adj_edges;
799 
800  auto set_singular_dofs = [&](auto &front_adj_edges, auto &front_vertices) {
801  MoFEMFunctionBegin;
802  auto &moab = mField.get_moab();
803 
804  double eps = 1;
805  double beta = 0;
806  CHKERR PetscOptionsGetScalar(PETSC_NULL, "-singularity_eps", &beta,
807  PETSC_NULL);
808  MOFEM_LOG("EP", Sev::inform) << "Singularity eps " << beta;
809  eps -= beta;
810 
811  for (auto edge : front_adj_edges) {
812  int num_nodes;
813  const EntityHandle *conn;
814  CHKERR moab.get_connectivity(edge, conn, num_nodes, false);
815  double coords[6];
816  CHKERR moab.get_coords(conn, num_nodes, coords);
817  const double dir[3] = {coords[3] - coords[0], coords[4] - coords[1],
818  coords[5] - coords[2]};
819  double dof[3] = {0, 0, 0};
820  if (front_vertices.find(conn[0]) != front_vertices.end()) {
821  for (int dd = 0; dd != 3; dd++) {
822  dof[dd] = -dir[dd] * eps;
823  }
824  } else if (front_vertices.find(conn[1]) != front_vertices.end()) {
825  for (int dd = 0; dd != 3; dd++) {
826  dof[dd] = +dir[dd] * eps;
827  }
828  }
829  for (_IT_GET_DOFS_FIELD_BY_NAME_AND_ENT_FOR_LOOP_(
830  mField, materialH1Positions, edge, dit)) {
831  const int idx = dit->get()->getEntDofIdx();
832  if (idx > 2) {
833  dit->get()->getFieldData() = 0;
834  } else {
835  dit->get()->getFieldData() = dof[idx];
836  }
837  }
838  }
839 
840  MoFEMFunctionReturn(0);
841  };
842 
843  if (setSingularity)
844  CHKERR set_singular_dofs(*frontAdjEdges, *frontVertices);
845 
846  MoFEMFunctionReturn(0);
847 }
848 
849 MoFEMErrorCode
850 EshelbianCore::addVolumeFiniteElement(const EntityHandle meshset) {
851  MoFEMFunctionBegin;
852 
853  // set finite element fields
854  auto add_field_to_fe = [this](const std::string fe,
855  const std::string field_name) {
856  MoFEMFunctionBegin;
857  CHKERR mField.modify_finite_element_add_field_row(fe, field_name);
858  CHKERR mField.modify_finite_element_add_field_col(fe, field_name);
859  CHKERR mField.modify_finite_element_add_field_data(fe, field_name);
860  MoFEMFunctionReturn(0);
861  };
862 
863  if (!mField.check_finite_element(elementVolumeName)) {
864  CHKERR mField.add_finite_element(elementVolumeName, MF_ZERO);
865  CHKERR mField.add_ents_to_finite_element_by_type(meshset, MBTET,
866  elementVolumeName);
867 
868  CHKERR add_field_to_fe(elementVolumeName, piolaStress);
869  CHKERR add_field_to_fe(elementVolumeName, bubbleField);
870  if (!noStretch)
871  CHKERR add_field_to_fe(elementVolumeName, stretchTensor);
872  CHKERR add_field_to_fe(elementVolumeName, rotAxis);
873  CHKERR add_field_to_fe(elementVolumeName, spatialL2Disp);
874  CHKERR add_field_to_fe(elementVolumeName, spatialH1Disp);
875  CHKERR add_field_to_fe(elementVolumeName, contactDisp);
876  CHKERR mField.modify_finite_element_add_field_data(elementVolumeName,
877  materialH1Positions);
878 
879  // build finite elements data structures
880  CHKERR mField.build_finite_elements(elementVolumeName);
881  }
882 
883  if (!mField.check_finite_element(materialVolumeElement)) {
884 
885  Range front_edges = get_range_from_block(mField, "FRONT", SPACE_DIM - 2);
886 
887  Range front_elements;
888  for (auto l = 0; l != frontLayers; ++l) {
889  Range front_elements_layer;
890  CHKERR mField.get_moab().get_adjacencies(front_edges, SPACE_DIM, true,
891  front_elements_layer,
892  moab::Interface::UNION);
893  front_elements.merge(front_elements_layer);
894  front_edges.clear();
895  CHKERR mField.get_moab().get_adjacencies(front_elements_layer,
896  SPACE_DIM - 2, true, front_edges,
897  moab::Interface::UNION);
898  }
899 
900  CHKERR mField.add_finite_element(materialVolumeElement, MF_ZERO);
901  CHKERR mField.add_ents_to_finite_element_by_type(front_elements, MBTET,
902  materialVolumeElement);
903  CHKERR add_field_to_fe(materialVolumeElement, eshelbyStress);
904  CHKERR add_field_to_fe(materialVolumeElement, materialL2Disp);
905  CHKERR mField.build_finite_elements(materialVolumeElement);
906  }
907 
908  MoFEMFunctionReturn(0);
909 }
910 
911 MoFEMErrorCode
912 EshelbianCore::addBoundaryFiniteElement(const EntityHandle meshset) {
913  MoFEMFunctionBegin;
914 
915  auto set_fe_adjacency = [&](auto fe_name) {
916  MoFEMFunctionBegin;
917  parentAdjSkeletonFunctionDim2 =
918  boost::make_shared<ParentFiniteElementAdjacencyFunctionSkeleton<2>>(
919  bitAdjParent, bitAdjParentMask, bitAdjEnt, bitAdjEntMask);
920  CHKERR mField.modify_finite_element_adjacency_table(
921  fe_name, MBTRI, *parentAdjSkeletonFunctionDim2);
922  MoFEMFunctionReturn(0);
923  };
924 
925  // set finite element fields
926  auto add_field_to_fe = [this](const std::string fe,
927  const std::string field_name) {
928  MoFEMFunctionBegin;
929  CHKERR mField.modify_finite_element_add_field_row(fe, field_name);
930  CHKERR mField.modify_finite_element_add_field_col(fe, field_name);
931  CHKERR mField.modify_finite_element_add_field_data(fe, field_name);
932  CHKERR mField.modify_finite_element_add_field_data(fe, materialH1Positions);
933  MoFEMFunctionReturn(0);
934  };
935 
936  if (!mField.check_finite_element(naturalBcElement)) {
937 
938  Range natural_bc_elements;
939  if (bcSpatialDispVecPtr) {
940  for (auto &v : *bcSpatialDispVecPtr) {
941  natural_bc_elements.merge(v.faces);
942  }
943  }
944  if (bcSpatialRotationVecPtr) {
945  for (auto &v : *bcSpatialRotationVecPtr) {
946  natural_bc_elements.merge(v.faces);
947  }
948  }
949  if (bcSpatialTraction) {
950  for (auto &v : *bcSpatialTraction) {
951  natural_bc_elements.merge(v.faces);
952  }
953  }
954 
955  CHKERR mField.add_finite_element(naturalBcElement, MF_ZERO);
956  CHKERR mField.add_ents_to_finite_element_by_type(natural_bc_elements, MBTRI,
957  naturalBcElement);
958  CHKERR add_field_to_fe(naturalBcElement, hybridSpatialDisp);
959  CHKERR mField.build_finite_elements(naturalBcElement);
960  }
961 
962  auto get_skin = [&](auto &body_ents) {
963  Skinner skin(&mField.get_moab());
964  Range skin_ents;
965  CHKERR skin.find_skin(0, body_ents, false, skin_ents);
966  return skin_ents;
967  };
968 
969  auto filter_true_skin = [&](auto &&skin) {
970  Range boundary_ents;
971  ParallelComm *pcomm =
972  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
973  CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
974  PSTATUS_NOT, -1, &boundary_ents);
975  return boundary_ents;
976  };
977 
978  if (!mField.check_finite_element(skinElement)) {
979 
980  Range body_ents;
981  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
982  auto skin = filter_true_skin(get_skin(body_ents));
983 
984  CHKERR mField.add_finite_element(skinElement, MF_ZERO);
985  CHKERR mField.add_ents_to_finite_element_by_type(skin, MBTRI, skinElement);
986  CHKERR mField.modify_finite_element_add_field_data(skinElement,
987  spatialH1Disp);
988  CHKERR mField.modify_finite_element_add_field_data(skinElement,
989  materialH1Positions);
990  CHKERR mField.modify_finite_element_add_field_data(skinElement,
991  hybridSpatialDisp);
992  CHKERR mField.modify_finite_element_add_field_data(skinElement,
993  hybridMaterialDisp);
994  CHKERR mField.modify_finite_element_add_field_data(skinElement,
995  contactDisp);
996  CHKERR mField.build_finite_elements(skinElement);
997  }
998 
999  if (!mField.check_finite_element(contactElement)) {
1000  if (contactFaces) {
1001  MOFEM_LOG("EP", Sev::inform)
1002  << "Contact elements " << contactFaces->size();
1003  CHKERR mField.add_finite_element(contactElement, MF_ZERO);
1004  CHKERR mField.add_ents_to_finite_element_by_type(*contactFaces, MBTRI,
1005  contactElement);
1006  CHKERR add_field_to_fe(contactElement, piolaStress);
1007  CHKERR add_field_to_fe(contactElement, spatialH1Disp);
1008  CHKERR add_field_to_fe(contactElement, hybridSpatialDisp);
1009  CHKERR add_field_to_fe(contactElement, contactDisp);
1010  CHKERR set_fe_adjacency(contactElement);
1011  CHKERR mField.build_finite_elements(contactElement);
1012  }
1013  }
1014 
1015  if (!mField.check_finite_element(skeletonElement)) {
1016  if (!skeletonFaces)
1017  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No skeleton faces");
1018  MOFEM_LOG("EP", Sev::inform)
1019  << "Skeleton elements " << skeletonFaces->size();
1020  CHKERR mField.add_finite_element(skeletonElement, MF_ZERO);
1021  CHKERR mField.add_ents_to_finite_element_by_type(*skeletonFaces, MBTRI,
1022  skeletonElement);
1023  CHKERR add_field_to_fe(skeletonElement, piolaStress);
1024  CHKERR add_field_to_fe(skeletonElement, hybridSpatialDisp);
1025  CHKERR add_field_to_fe(contactElement, contactDisp);
1026  CHKERR set_fe_adjacency(skeletonElement);
1027  CHKERR mField.build_finite_elements(skeletonElement);
1028  }
1029 
1030  if (!mField.check_finite_element(materialSkeletonElement)) {
1031 
1032  Range front_edges = get_range_from_block(mField, "FRONT", SPACE_DIM - 2);
1033 
1034  Range front_elements;
1035  for (auto l = 0; l != frontLayers; ++l) {
1036  Range front_elements_layer;
1037  CHKERR mField.get_moab().get_adjacencies(front_edges, SPACE_DIM, true,
1038  front_elements_layer,
1039  moab::Interface::UNION);
1040  front_elements.merge(front_elements_layer);
1041  front_edges.clear();
1042  CHKERR mField.get_moab().get_adjacencies(front_elements_layer,
1043  SPACE_DIM - 2, true, front_edges,
1044  moab::Interface::UNION);
1045  }
1046  Range body_ents;
1047  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
1048  Range front_elements_faces;
1049  CHKERR mField.get_moab().get_adjacencies(front_elements, SPACE_DIM - 1,
1050  true, front_elements_faces,
1051  moab::Interface::UNION);
1052  auto body_skin = filter_true_skin(get_skin(body_ents));
1053  auto front_elements_skin = filter_true_skin(get_skin(front_elements));
1054  Range material_skeleton_faces =
1055  subtract(front_elements_faces, front_elements_skin);
1056  material_skeleton_faces.merge(intersect(front_elements_skin, body_skin));
1057 
1058 #ifndef NDEBUG
1059  CHKERR save_range(mField.get_moab(), "front_elements.vtk", front_elements);
1060  CHKERR save_range(mField.get_moab(), "front_elements_skin.vtk",
1061  front_elements_skin);
1062  CHKERR save_range(mField.get_moab(), "material_skeleton_faces.vtk",
1063  material_skeleton_faces);
1064 #endif
1065 
1066  CHKERR mField.add_finite_element(materialSkeletonElement, MF_ZERO);
1067  CHKERR mField.add_ents_to_finite_element_by_type(
1068  material_skeleton_faces, MBTRI, materialSkeletonElement);
1069  CHKERR add_field_to_fe(materialSkeletonElement, eshelbyStress);
1070  CHKERR add_field_to_fe(materialSkeletonElement, hybridMaterialDisp);
1071  CHKERR set_fe_adjacency(materialSkeletonElement);
1072  CHKERR mField.build_finite_elements(materialSkeletonElement);
1073  }
1074 
1075  if (!mField.check_finite_element(crackElement)) {
1076  auto crack_faces = get_range_from_block(mField, "CRACK", SPACE_DIM - 1);
1077  CHKERR mField.add_finite_element(crackElement, MF_ZERO);
1078  CHKERR mField.add_ents_to_finite_element_by_type(crack_faces, MBTRI,
1079  crackElement);
1080  CHKERR add_field_to_fe(crackElement, eshelbyStress);
1081  CHKERR add_field_to_fe(crackElement, hybridMaterialDisp);
1082 
1083  CHKERR mField.modify_finite_element_add_field_data(crackElement,
1084  spatialH1Disp);
1085  CHKERR mField.modify_finite_element_add_field_data(crackElement,
1086  materialH1Positions);
1087  CHKERR mField.modify_finite_element_add_field_data(crackElement,
1088  hybridSpatialDisp);
1089  CHKERR mField.modify_finite_element_add_field_data(crackElement,
1090  hybridMaterialDisp);
1091  CHKERR mField.modify_finite_element_add_field_data(crackElement,
1092  contactDisp);
1093 
1094  CHKERR set_fe_adjacency(crackElement);
1095  CHKERR mField.build_finite_elements(crackElement);
1096  }
1097 
1098  MoFEMFunctionReturn(0);
1099 }
1100 
1101 MoFEMErrorCode EshelbianCore::addDMs(const BitRefLevel bit,
1102  const EntityHandle meshset) {
1103  MoFEMFunctionBegin;
1104 
1105  // find adjacencies between finite elements and dofs
1106  CHKERR mField.build_adjacencies(bit, QUIET);
1107 
1108  // Create coupled problem
1109  dM = createDM(mField.get_comm(), "DMMOFEM");
1110  CHKERR DMMoFEMCreateMoFEM(dM, &mField, "ESHELBY_PLASTICITY", bit,
1111  BitRefLevel().set());
1112  CHKERR DMMoFEMSetDestroyProblem(dM, PETSC_TRUE);
1113  CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
1114  CHKERR DMMoFEMAddElement(dM, elementVolumeName);
1115  CHKERR DMMoFEMAddElement(dM, naturalBcElement);
1116  CHKERR DMMoFEMAddElement(dM, skeletonElement);
1117  CHKERR DMMoFEMAddElement(dM, contactElement);
1118  CHKERR DMMoFEMAddElement(dM, crackElement);
1119  CHKERR DMMoFEMAddElement(dM, skinElement);
1120 
1121  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
1122  CHKERR DMSetUp(dM);
1123  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
1124 
1125  auto remove_dofs_on_broken_skin = [&](const std::string prb_name) {
1126  MoFEMFunctionBegin;
1127  for (int d : {0, 1, 2}) {
1128  std::vector<boost::weak_ptr<NumeredDofEntity>> dofs_to_remove;
1129  CHKERR mField.getInterface<ProblemsManager>()
1130  ->getSideDofsOnBrokenSpaceEntities(
1131  dofs_to_remove, prb_name, ROW, piolaStress,
1132  (*bcSpatialFreeTraction)[d], SPACE_DIM, d, d);
1133  // remove piola dofs, i.e. traction free boundary
1134  CHKERR mField.getInterface<ProblemsManager>()->removeDofs(prb_name, ROW,
1135  dofs_to_remove);
1136  CHKERR mField.getInterface<ProblemsManager>()->removeDofs(prb_name, COL,
1137  dofs_to_remove);
1138  }
1139  MoFEMFunctionReturn(0);
1140  };
1141  CHKERR remove_dofs_on_broken_skin("ESHELBY_PLASTICITY");
1142 
1143  // Create elastic sub-problem
1144  dmElastic = createDM(mField.get_comm(), "DMMOFEM");
1145  CHKERR DMMoFEMCreateSubDM(dmElastic, dM, "ELASTIC_PROBLEM");
1146  CHKERR DMMoFEMSetDestroyProblem(dmElastic, PETSC_TRUE);
1147  CHKERR DMMoFEMAddSubFieldRow(dmElastic, piolaStress);
1148  CHKERR DMMoFEMAddSubFieldRow(dmElastic, bubbleField);
1149  CHKERR DMMoFEMAddSubFieldRow(dmElastic, spatialL2Disp);
1150  CHKERR DMMoFEMAddSubFieldRow(dmElastic, rotAxis);
1151  if (!noStretch) {
1152  CHKERR DMMoFEMAddSubFieldRow(dmElastic, stretchTensor);
1153  }
1154  CHKERR DMMoFEMAddSubFieldRow(dmElastic, hybridSpatialDisp);
1155  CHKERR DMMoFEMAddSubFieldRow(dmElastic, contactDisp);
1156  CHKERR DMMoFEMAddElement(dmElastic, elementVolumeName);
1157  CHKERR DMMoFEMAddElement(dmElastic, naturalBcElement);
1158  CHKERR DMMoFEMAddElement(dmElastic, skeletonElement);
1159  CHKERR DMMoFEMAddElement(dmElastic, contactElement);
1160  CHKERR DMMoFEMAddElement(dmElastic, skinElement);
1161  CHKERR DMMoFEMSetSquareProblem(dmElastic, PETSC_TRUE);
1162  CHKERR DMMoFEMSetIsPartitioned(dmElastic, PETSC_TRUE);
1163  CHKERR DMSetUp(dmElastic);
1164 
1165  // dmMaterial = createDM(mField.get_comm(), "DMMOFEM");
1166  // CHKERR DMMoFEMCreateSubDM(dmMaterial, dM, "MATERIAL_PROBLEM");
1167  // CHKERR DMMoFEMSetDestroyProblem(dmMaterial, PETSC_TRUE);
1168  // CHKERR DMMoFEMAddSubFieldRow(dmMaterial, eshelbyStress);
1169  // CHKERR DMMoFEMAddSubFieldRow(dmMaterial, materialL2Disp);
1170  // CHKERR DMMoFEMAddElement(dmMaterh elementVolumeName);
1171  // CHKERR DMMoFEMAddElement(dmMaterial, naturalBcElement);
1172  // CHKERR DMMoFEMAddElement(dmMaterial, skinElement);
1173  // CHKERR DMMoFEMSetSquareProblem(dmMaterial, PETSC_TRUE);
1174  // CHKERR DMMoFEMSetIsPartitioned(dmMaterial, PETSC_TRUE);
1175  // CHKERR DMSetUp(dmMaterial);
1176 
1177  auto set_zero_block = [&]() {
1178  MoFEMFunctionBegin;
1179  if (!noStretch) {
1180  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1181  "ELASTIC_PROBLEM", spatialL2Disp, stretchTensor);
1182  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1183  "ELASTIC_PROBLEM", stretchTensor, spatialL2Disp);
1184  }
1185  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1186  "ELASTIC_PROBLEM", spatialL2Disp, rotAxis);
1187  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1188  "ELASTIC_PROBLEM", rotAxis, spatialL2Disp);
1189  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1190  "ELASTIC_PROBLEM", spatialL2Disp, bubbleField);
1191  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1192  "ELASTIC_PROBLEM", bubbleField, spatialL2Disp);
1193  if (!noStretch) {
1194  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1195  "ELASTIC_PROBLEM", bubbleField, bubbleField);
1196  CHKERR
1197  mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1198  "ELASTIC_PROBLEM", piolaStress, piolaStress);
1199  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1200  "ELASTIC_PROBLEM", bubbleField, piolaStress);
1201  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
1202  "ELASTIC_PROBLEM", piolaStress, bubbleField);
1203  }
1204  MoFEMFunctionReturn(0);
1205  };
1206 
1207  auto set_section = [&]() {
1208  MoFEMFunctionBegin;
1209  PetscSection section;
1210  CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
1211  &section);
1212  CHKERR DMSetSection(dmElastic, section);
1213  CHKERR DMSetGlobalSection(dmElastic, section);
1214  CHKERR PetscSectionDestroy(&section);
1215  MoFEMFunctionReturn(0);
1216  };
1217 
1218  CHKERR set_zero_block();
1219  CHKERR set_section();
1220 
1221  dmPrjSpatial = createDM(mField.get_comm(), "DMMOFEM");
1222  CHKERR DMMoFEMCreateSubDM(dmPrjSpatial, dM, "PROJECT_SPATIAL");
1223  CHKERR DMMoFEMSetDestroyProblem(dmPrjSpatial, PETSC_TRUE);
1224  CHKERR DMMoFEMAddSubFieldRow(dmPrjSpatial, spatialH1Disp);
1225  CHKERR DMMoFEMAddElement(dmPrjSpatial, elementVolumeName);
1226  CHKERR DMMoFEMSetSquareProblem(dmPrjSpatial, PETSC_TRUE);
1227  CHKERR DMMoFEMSetIsPartitioned(dmPrjSpatial, PETSC_TRUE);
1228  CHKERR DMSetUp(dmPrjSpatial);
1229 
1230  CHKERR mField.getInterface<BcManager>()
1231  ->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
1232  "PROJECT_SPATIAL", spatialH1Disp, true, false);
1233 
1234  MoFEMFunctionReturn(0);
1235 }
1236 
1237 BcDisp::BcDisp(std::string name, std::vector<double> &attr, Range &faces)
1238  : blockName(name), faces(faces) {
1239  vals.resize(3, false);
1240  flags.resize(3, false);
1241  for (int ii = 0; ii != 3; ++ii) {
1242  vals[ii] = attr[ii];
1243  flags[ii] = static_cast<int>(attr[ii + 3]);
1244  }
1245 
1246  MOFEM_LOG("EP", Sev::inform) << "Add BCDisp " << name;
1247  MOFEM_LOG("EP", Sev::inform)
1248  << "Add BCDisp vals " << vals[0] << " " << vals[1] << " " << vals[2];
1249  MOFEM_LOG("EP", Sev::inform)
1250  << "Add BCDisp flags " << flags[0] << " " << flags[1] << " " << flags[2];
1251  MOFEM_LOG("EP", Sev::inform) << "Add BCDisp nb. of faces " << faces.size();
1252 }
1253 
1254 BcRot::BcRot(std::string name, std::vector<double> &attr, Range &faces)
1255  : blockName(name), faces(faces) {
1256  vals.resize(3, false);
1257  for (int ii = 0; ii != 3; ++ii) {
1258  vals[ii] = attr[ii];
1259  }
1260  theta = attr[3];
1261 }
1262 
1263 TractionBc::TractionBc(std::string name, std::vector<double> &attr,
1264  Range &faces)
1265  : blockName(name), faces(faces) {
1266  vals.resize(3, false);
1267  flags.resize(3, false);
1268  for (int ii = 0; ii != 3; ++ii) {
1269  vals[ii] = attr[ii];
1270  flags[ii] = static_cast<int>(attr[ii + 3]);
1271  }
1272 
1273  MOFEM_LOG("EP", Sev::inform) << "Add BCForce " << name;
1274  MOFEM_LOG("EP", Sev::inform)
1275  << "Add BCForce vals " << vals[0] << " " << vals[1] << " " << vals[2];
1276  MOFEM_LOG("EP", Sev::inform)
1277  << "Add BCForce flags " << flags[0] << " " << flags[1] << " " << flags[2];
1278  MOFEM_LOG("EP", Sev::inform) << "Add BCForce nb. of faces " << faces.size();
1279 }
1280 
1281 MoFEMErrorCode
1282 EshelbianCore::getTractionFreeBc(const EntityHandle meshset,
1283  boost::shared_ptr<TractionFreeBc> &bc_ptr,
1284  const std::string contact_set_name) {
1285  MoFEMFunctionBegin;
1286 
1287  // get skin from all tets
1288  Range tets;
1289  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
1290  Range tets_skin_part;
1291  Skinner skin(&mField.get_moab());
1292  CHKERR skin.find_skin(0, tets, false, tets_skin_part);
1293  ParallelComm *pcomm =
1294  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
1295  Range tets_skin;
1296  CHKERR pcomm->filter_pstatus(tets_skin_part,
1297  PSTATUS_SHARED | PSTATUS_MULTISHARED,
1298  PSTATUS_NOT, -1, &tets_skin);
1299 
1300  bc_ptr->resize(3);
1301  for (int dd = 0; dd != 3; ++dd)
1302  (*bc_ptr)[dd] = tets_skin;
1303 
1304  // Do not remove dofs on which traction is applied
1305  if (bcSpatialDispVecPtr)
1306  for (auto &v : *bcSpatialDispVecPtr) {
1307  if (v.flags[0])
1308  (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
1309  if (v.flags[1])
1310  (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
1311  if (v.flags[2])
1312  (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
1313  }
1314 
1315  // Do not remove dofs on which rotation is applied
1316  if (bcSpatialRotationVecPtr)
1317  for (auto &v : *bcSpatialRotationVecPtr) {
1318  (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
1319  (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
1320  (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
1321  }
1322 
1323  if (bcSpatialTraction)
1324  for (auto &v : *bcSpatialTraction) {
1325  (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
1326  (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
1327  (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
1328  }
1329 
1330  // remove contact
1331  for (auto m : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
1332  std::regex((boost::format("%s(.*)") % contact_set_name).str()))) {
1333  Range faces;
1334  CHKERR m->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
1335  true);
1336  (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
1337  (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
1338  (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
1339  }
1340 
1341  MoFEMFunctionReturn(0);
1342 }
1343 
1344 /**
1345  * @brief Set integration rule on element
1346  * \param order on row
1347  * \param order on column
1348  * \param order on data
1349  *
1350  * Use maximal oder on data in order to determine integration rule
1351  *
1352  */
1353 struct VolRule {
1354  int operator()(int p_row, int p_col, int p_data) const {
1355  // if (EshelbianCore::gradApproximator > MODERATE_ROT)
1356  // return p_data + p_data + (p_data - 1);
1357  // else
1358  return 2 * p_data + 1;
1359  }
1360 };
1361 
1362 struct FaceRule {
1363  int operator()(int p_row, int p_col, int p_data) const {
1364  return 2 * (p_data + 1);
1365  }
1366 };
1367 
1368 CGGUserPolynomialBase::CGGUserPolynomialBase(
1369  boost::shared_ptr<CachePhi> cache_phi_otr)
1370  : TetPolynomialBase(), cachePhiPtr(cache_phi_otr) {}
1371 
1372 MoFEMErrorCode CGGUserPolynomialBase::query_interface(
1373  boost::typeindex::type_index type_index,
1374  BaseFunctionUnknownInterface **iface) const {
1375  *iface = const_cast<CGGUserPolynomialBase *>(this);
1376  return 0;
1377 }
1378 
1379 MoFEMErrorCode
1380 CGGUserPolynomialBase::getValue(MatrixDouble &pts,
1381  boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
1382  MoFEMFunctionBeginHot;
1383 
1384  this->cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
1385 
1386  int nb_gauss_pts = pts.size2();
1387  if (!nb_gauss_pts) {
1388  MoFEMFunctionReturnHot(0);
1389  }
1390 
1391  if (pts.size1() < 3) {
1392  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1393  "Wrong dimension of pts, should be at least 3 rows with "
1394  "coordinates");
1395  }
1396 
1397  const auto base = this->cTx->bAse;
1398  EntitiesFieldData &data = this->cTx->dAta;
1399 
1400  switch (this->cTx->sPace) {
1401  case HDIV:
1402  CHKERR getValueHdivForCGGBubble(pts);
1403  break;
1404  case L2:
1405  data.dataOnEntities[MBVERTEX][0].getN(base).resize(nb_gauss_pts, 4, false);
1406  CHKERR Tools::shapeFunMBTET(
1407  &*data.dataOnEntities[MBVERTEX][0].getN(base).data().begin(),
1408  &pts(0, 0), &pts(1, 0), &pts(2, 0), nb_gauss_pts);
1409  data.dataOnEntities[MBVERTEX][0].getDiffN(base).resize(4, 3, false);
1410  std::copy(Tools::diffShapeFunMBTET.begin(), Tools::diffShapeFunMBTET.end(),
1411  data.dataOnEntities[MBVERTEX][0].getDiffN(base).data().begin());
1412  this->cTx->basePolynomialsType0 = Legendre_polynomials;
1413  CHKERR getValueL2AinsworthBase(pts);
1414  break;
1415  default:
1416  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
1417  }
1418 
1419  MoFEMFunctionReturnHot(0);
1420 }
1421 
1422 MoFEMErrorCode
1423 CGGUserPolynomialBase::getValueHdivForCGGBubble(MatrixDouble &pts) {
1424  MoFEMFunctionBegin;
1425 
1426  // This should be used only in case USER_BASE is selected
1427  if (this->cTx->bAse != USER_BASE) {
1428  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1429  "Wrong base, should be USER_BASE");
1430  }
1431  // get access to data structures on element
1432  EntitiesFieldData &data = this->cTx->dAta;
1433  // Get approximation order on element. Note that bubble functions are only
1434  // on tetrahedron.
1435  const int order = data.dataOnEntities[MBTET][0].getOrder();
1436  /// number of integration points
1437  const int nb_gauss_pts = pts.size2();
1438 
1439  auto check_cache = [this](int order, int nb_gauss_pts) -> bool {
1440  if (cachePhiPtr) {
1441  return cachePhiPtr->get<0>() == order &&
1442  cachePhiPtr->get<1>() == nb_gauss_pts;
1443  } else {
1444  return false;
1445  }
1446  };
1447 
1448  if (check_cache(order, nb_gauss_pts)) {
1449  auto &mat = cachePhiPtr->get<2>();
1450  auto &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
1451  phi.resize(mat.size1(), mat.size2(), false);
1452  noalias(phi) = mat;
1453 
1454  } else {
1455  // calculate shape functions, i.e. barycentric coordinates
1456  shapeFun.resize(nb_gauss_pts, 4, false);
1457  CHKERR ShapeMBTET(&*shapeFun.data().begin(), &pts(0, 0), &pts(1, 0),
1458  &pts(2, 0), nb_gauss_pts);
1459  // derivatives of shape functions
1460  double diff_shape_fun[12];
1461  CHKERR ShapeDiffMBTET(diff_shape_fun);
1462 
1463  const int nb_base_functions = NBVOLUMETET_CCG_BUBBLE(order);
1464  // get base functions and set size
1465  MatrixDouble &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
1466  phi.resize(nb_gauss_pts, 9 * nb_base_functions, false);
1467  // finally calculate base functions
1468  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_phi(
1469  &phi(0, 0), &phi(0, 1), &phi(0, 2),
1470 
1471  &phi(0, 3), &phi(0, 4), &phi(0, 5),
1472 
1473  &phi(0, 6), &phi(0, 7), &phi(0, 8));
1474  CHKERR CGG_BubbleBase_MBTET(order, &shapeFun(0, 0), diff_shape_fun, t_phi,
1475  nb_gauss_pts);
1476 
1477  if (cachePhiPtr) {
1478  cachePhiPtr->get<0>() = order;
1479  cachePhiPtr->get<1>() = nb_gauss_pts;
1480  cachePhiPtr->get<2>().resize(phi.size1(), phi.size2(), false);
1481  noalias(cachePhiPtr->get<2>()) = phi;
1482  }
1483  }
1484 
1485  MoFEMFunctionReturn(0);
1486 }
1487 
1488 MoFEMErrorCode EshelbianCore::setBaseVolumeElementOps(
1489  const int tag, const bool do_rhs, const bool do_lhs,
1490  boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe) {
1491  MoFEMFunctionBegin;
1492  fe = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
1493 
1494  auto bubble_cache =
1495  boost::make_shared<CGGUserPolynomialBase::CachePhi>(0, 0, MatrixDouble());
1496  fe->getUserPolynomialBase() =
1497  boost::make_shared<CGGUserPolynomialBase>(bubble_cache);
1498  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
1499  fe->getOpPtrVector(), {HDIV, H1, L2}, materialH1Positions, frontAdjEdges);
1500 
1501  // set integration rule
1502  fe->getRuleHook = VolRule();
1503 
1504  if (!dataAtPts) {
1505  dataAtPts =
1506  boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
1507  dataAtPts->physicsPtr = physicalEquations;
1508  }
1509 
1510  // calculate fields values
1511  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
1512  piolaStress, dataAtPts->getApproxPAtPts()));
1513  fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
1514  bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
1515  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
1516  piolaStress, dataAtPts->getDivPAtPts()));
1517 
1518  if (noStretch) {
1519  fe->getOpPtrVector().push_back(
1520  physicalEquations->returnOpCalculateStretchFromStress(
1521  dataAtPts, physicalEquations));
1522  } else {
1523  fe->getOpPtrVector().push_back(
1524  new OpCalculateTensor2SymmetricFieldValues<3>(
1525  stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
1526  }
1527 
1528  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1529  rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
1530  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1531  spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
1532 
1533  // velocities
1534  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
1535  spatialL2Disp, dataAtPts->getSmallWL2DotAtPts(), MBTET));
1536  if (noStretch) {
1537  } else {
1538  fe->getOpPtrVector().push_back(
1539  new OpCalculateTensor2SymmetricFieldValuesDot<3>(
1540  stretchTensor, dataAtPts->getLogStretchDotTensorAtPts(), MBTET));
1541  }
1542  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
1543  rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
1544 
1545  // acceleration
1546  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
1547  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDotDot<3>(
1548  spatialL2Disp, dataAtPts->getSmallWL2DotDotAtPts(), MBTET));
1549  }
1550 
1551  // H1 displacements
1552  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1553  spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
1554  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
1555  spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
1556 
1557  // calculate other derived quantities
1558  fe->getOpPtrVector().push_back(
1559  new OpCalculateRotationAndSpatialGradient(dataAtPts));
1560 
1561  // evaluate integration points
1562  if (noStretch) {
1563  } else {
1564  fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
1565  tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
1566  }
1567 
1568  MoFEMFunctionReturn(0);
1569 }
1570 
1571 MoFEMErrorCode EshelbianCore::setVolumeElementOps(
1572  const int tag, const bool add_elastic, const bool add_material,
1573  boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_rhs,
1574  boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_lhs) {
1575  MoFEMFunctionBegin;
1576 
1577  /** Contact requires that body is marked */
1578  auto get_body_range = [this](auto name, int dim) {
1579  std::map<int, Range> map;
1580 
1581  for (auto m_ptr :
1582  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
1583 
1584  (boost::format("%s(.*)") % name).str()
1585 
1586  ))
1587 
1588  ) {
1589  Range ents;
1590  CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
1591  dim, ents, true),
1592  "by dim");
1593  map[m_ptr->getMeshsetId()] = ents;
1594  }
1595 
1596  return map;
1597  };
1598 
1599  auto rule_contact = [](int, int, int o) { return -1; };
1600  auto refine = Tools::refineTriangle(contactRefinementLevels);
1601 
1602  auto set_rule_contact = [refine](
1603 
1604  ForcesAndSourcesCore *fe_raw_ptr, int order_row,
1605  int order_col, int order_data
1606 
1607  ) {
1608  MoFEMFunctionBegin;
1609  auto rule = 2 * order_data;
1610  fe_raw_ptr->gaussPts = Tools::refineTriangleIntegrationPts(rule, refine);
1611  MoFEMFunctionReturn(0);
1612  };
1613 
1614  auto time_scale = boost::make_shared<TimeScale>();
1615 
1616  // Right hand side
1617  CHKERR setBaseVolumeElementOps(tag, true, false, fe_rhs);
1618 
1619  // elastic
1620  if (add_elastic) {
1621  fe_rhs->getOpPtrVector().push_back(
1622  new OpSpatialEquilibrium(spatialL2Disp, dataAtPts, alphaW, alphaRho));
1623  fe_rhs->getOpPtrVector().push_back(
1624  new OpSpatialRotation(rotAxis, dataAtPts));
1625  if (noStretch) {
1626  // do nothing - no stretch approximation
1627  } else {
1628  fe_rhs->getOpPtrVector().push_back(
1629  physicalEquations->returnOpSpatialPhysical(stretchTensor, dataAtPts,
1630  alphaU));
1631  }
1632  fe_rhs->getOpPtrVector().push_back(
1633  new OpSpatialConsistencyP(piolaStress, dataAtPts));
1634  fe_rhs->getOpPtrVector().push_back(
1635  new OpSpatialConsistencyBubble(bubbleField, dataAtPts));
1636  fe_rhs->getOpPtrVector().push_back(
1637  new OpSpatialConsistencyDivTerm(piolaStress, dataAtPts));
1638 
1639  auto set_hybridisation = [&](auto &pip) {
1640  MoFEMFunctionBegin;
1641 
1642  using BoundaryEle =
1643  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
1644  using EleOnSide =
1645  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
1646  using SideEleOp = EleOnSide::UserDataOperator;
1647  using BdyEleOp = BoundaryEle::UserDataOperator;
1648 
1649  // First: Iterate over skeleton FEs adjacent to Domain FEs
1650  // Note: BoundaryEle, i.e. uses skeleton interation rule
1651  auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
1652  mField, skeletonElement, SPACE_DIM - 1, Sev::noisy);
1653  op_loop_skeleton_side->getSideFEPtr()->getRuleHook = FaceRule();
1654  CHKERR EshelbianPlasticity::
1655  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
1656  op_loop_skeleton_side->getOpPtrVector(), {L2},
1657  materialH1Positions, frontAdjEdges);
1658 
1659  // Second: Iterate over domain FEs adjacent to skelton, particularly one
1660  // domain element.
1661  auto broken_data_ptr =
1662  boost::make_shared<std::vector<BrokenBaseSideData>>();
1663  // Note: EleOnSide, i.e. uses on domain projected skeleton rule
1664  auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
1665  mField, elementVolumeName, SPACE_DIM, Sev::noisy);
1666  op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
1667  boost::make_shared<CGGUserPolynomialBase>();
1668  CHKERR
1669  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
1670  op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
1671  materialH1Positions, frontAdjEdges);
1672  op_loop_domain_side->getOpPtrVector().push_back(
1673  new OpGetBrokenBaseSideData<SideEleOp>(piolaStress, broken_data_ptr));
1674  auto flux_mat_ptr = boost::make_shared<MatrixDouble>();
1675  op_loop_domain_side->getOpPtrVector().push_back(
1676  new OpCalculateHVecTensorField<SPACE_DIM, SPACE_DIM>(piolaStress,
1677  flux_mat_ptr));
1678  op_loop_domain_side->getOpPtrVector().push_back(
1679  new OpSetFlux<SideEleOp>(broken_data_ptr, flux_mat_ptr));
1680 
1681  // Assemble on skeleton
1682  op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
1683  using OpC_dHybrid = FormsIntegrators<BdyEleOp>::Assembly<A>::LinearForm<
1684  GAUSS>::OpBrokenSpaceConstrainDHybrid<SPACE_DIM>;
1685  using OpC_dBroken = FormsIntegrators<BdyEleOp>::Assembly<A>::LinearForm<
1686  GAUSS>::OpBrokenSpaceConstrainDFlux<SPACE_DIM>;
1687  op_loop_skeleton_side->getOpPtrVector().push_back(
1688  new OpC_dHybrid(hybridSpatialDisp, broken_data_ptr, 1.));
1689  auto hybrid_ptr = boost::make_shared<MatrixDouble>();
1690  op_loop_skeleton_side->getOpPtrVector().push_back(
1691  new OpCalculateVectorFieldValues<SPACE_DIM>(hybridSpatialDisp,
1692  hybrid_ptr));
1693  op_loop_skeleton_side->getOpPtrVector().push_back(
1694  new OpC_dBroken(broken_data_ptr, hybrid_ptr, 1.));
1695 
1696  // Add skeleton to domain pipeline
1697  pip.push_back(op_loop_skeleton_side);
1698 
1699  MoFEMFunctionReturn(0);
1700  };
1701 
1702  auto set_contact = [&](auto &pip) {
1703  MoFEMFunctionBegin;
1704 
1705  using BoundaryEle =
1706  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
1707  using EleOnSide =
1708  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
1709  using SideEleOp = EleOnSide::UserDataOperator;
1710  using BdyEleOp = BoundaryEle::UserDataOperator;
1711 
1712  // First: Iterate over skeleton FEs adjacent to Domain FEs
1713  // Note: BoundaryEle, i.e. uses skeleton interation rule
1714  auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
1715  mField, contactElement, SPACE_DIM - 1, Sev::noisy);
1716 
1717  op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
1718  op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
1719  CHKERR EshelbianPlasticity::
1720  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
1721  op_loop_skeleton_side->getOpPtrVector(), {L2},
1722  materialH1Positions, frontAdjEdges);
1723 
1724  // Second: Iterate over domain FEs adjacent to skelton, particularly
1725  // one domain element.
1726  auto broken_data_ptr =
1727  boost::make_shared<std::vector<BrokenBaseSideData>>();
1728 
1729  // Data storing contact fields
1730  auto contact_common_data_ptr =
1731  boost::make_shared<ContactOps::CommonData>();
1732 
1733  auto add_ops_domain_side = [&](auto &pip) {
1734  MoFEMFunctionBegin;
1735  // Note: EleOnSide, i.e. uses on domain projected skeleton rule
1736  auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
1737  mField, elementVolumeName, SPACE_DIM, Sev::noisy);
1738  op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
1739  boost::make_shared<CGGUserPolynomialBase>();
1740  CHKERR
1741  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
1742  op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
1743  materialH1Positions, frontAdjEdges);
1744  op_loop_domain_side->getOpPtrVector().push_back(
1745  new OpGetBrokenBaseSideData<SideEleOp>(piolaStress,
1746  broken_data_ptr));
1747  op_loop_domain_side->getOpPtrVector().push_back(
1748  new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
1749  piolaStress, contact_common_data_ptr->contactTractionPtr()));
1750  pip.push_back(op_loop_domain_side);
1751  MoFEMFunctionReturn(0);
1752  };
1753 
1754  auto add_ops_contact_rhs = [&](auto &pip) {
1755  MoFEMFunctionBegin;
1756  // get body id and SDF range
1757  auto contact_sfd_map_range_ptr =
1758  boost::make_shared<std::map<int, Range>>(
1759  get_body_range("CONTACT_SDF", SPACE_DIM - 1));
1760 
1761  pip.push_back(new OpCalculateVectorFieldValues<3>(
1762  contactDisp, contact_common_data_ptr->contactDispPtr()));
1763  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
1764  pip.push_back(
1765  new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
1766  pip.push_back(new OpTreeSearch(
1767  contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
1768  get_range_from_block(mField, "CONTACT", SPACE_DIM - 1), nullptr,
1769  nullptr));
1770  pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryL2Rhs(
1771  contactDisp, contact_common_data_ptr, contactTreeRhs,
1772  contact_sfd_map_range_ptr));
1773  pip.push_back(
1774  new EshelbianPlasticity::OpConstrainBoundaryHDivRhs<A, GAUSS>(
1775  broken_data_ptr, contact_common_data_ptr, contactTreeRhs));
1776 
1777  MoFEMFunctionReturn(0);
1778  };
1779 
1780  // push ops to face/side pipeline
1781  CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
1782  CHKERR add_ops_contact_rhs(op_loop_skeleton_side->getOpPtrVector());
1783 
1784  // Add skeleton to domain pipeline
1785  pip.push_back(op_loop_skeleton_side);
1786 
1787  MoFEMFunctionReturn(0);
1788  };
1789 
1790  CHKERR set_hybridisation(fe_rhs->getOpPtrVector());
1791  CHKERR set_contact(fe_rhs->getOpPtrVector());
1792 
1793  // Body forces
1794  using BodyNaturalBC =
1795  NaturalBC<VolumeElementForcesAndSourcesCore::UserDataOperator>::
1796  Assembly<PETSC>::LinearForm<GAUSS>;
1797  using OpBodyForce =
1798  BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
1799  CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
1800  fe_rhs->getOpPtrVector(), mField, spatialL2Disp, {time_scale},
1801  "BODY_FORCE", Sev::inform);
1802  }
1803 
1804  // Left hand side
1805  CHKERR setBaseVolumeElementOps(tag, true, true, fe_lhs);
1806 
1807  // elastic
1808  if (add_elastic) {
1809 
1810  const bool symmetric_system =
1811  (gradApproximator <= MODERATE_ROT && rotSelector == SMALL_ROT);
1812 
1813  if (noStretch) {
1814  fe_lhs->getOpPtrVector().push_back(
1815  new OpSpatialConsistency_dP_dP(piolaStress, piolaStress, dataAtPts));
1816  fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_dP(
1817  bubbleField, piolaStress, dataAtPts));
1818  fe_lhs->getOpPtrVector().push_back(
1819  new OpSpatialConsistency_dBubble_dBubble(bubbleField, bubbleField,
1820  dataAtPts));
1821  } else {
1822  fe_lhs->getOpPtrVector().push_back(
1823  physicalEquations->returnOpSpatialPhysical_du_du(
1824  stretchTensor, stretchTensor, dataAtPts, alphaU));
1825  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
1826  stretchTensor, piolaStress, dataAtPts, true));
1827  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
1828  stretchTensor, bubbleField, dataAtPts, true));
1829  if (!symmetric_system) {
1830  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
1831  stretchTensor, rotAxis, dataAtPts, false));
1832  }
1833  }
1834 
1835  fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
1836  spatialL2Disp, piolaStress, dataAtPts, true));
1837  fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
1838  spatialL2Disp, spatialL2Disp, dataAtPts, alphaW, alphaRho));
1839 
1840  fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dP_domega(
1841  piolaStress, rotAxis, dataAtPts, symmetric_system));
1842  fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
1843  bubbleField, rotAxis, dataAtPts, symmetric_system));
1844 
1845  if (!symmetric_system) {
1846  fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dP(
1847  rotAxis, piolaStress, dataAtPts, false));
1848  fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dBubble(
1849  rotAxis, bubbleField, dataAtPts, false));
1850  fe_lhs->getOpPtrVector().push_back(
1851  new OpSpatialRotation_domega_domega(rotAxis, rotAxis, dataAtPts));
1852  }
1853 
1854  auto set_hybridisation = [&](auto &pip) {
1855  MoFEMFunctionBegin;
1856 
1857  using BoundaryEle =
1858  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
1859  using EleOnSide =
1860  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
1861  using SideEleOp = EleOnSide::UserDataOperator;
1862  using BdyEleOp = BoundaryEle::UserDataOperator;
1863 
1864  // First: Iterate over skeleton FEs adjacent to Domain FEs
1865  // Note: BoundaryEle, i.e. uses skeleton interation rule
1866  auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
1867  mField, skeletonElement, SPACE_DIM - 1, Sev::noisy);
1868  op_loop_skeleton_side->getSideFEPtr()->getRuleHook = FaceRule();
1869  CHKERR EshelbianPlasticity::
1870  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
1871  op_loop_skeleton_side->getOpPtrVector(), {L2},
1872  materialH1Positions, frontAdjEdges);
1873 
1874  // Second: Iterate over domain FEs adjacent to skelton, particularly one
1875  // domain element.
1876  auto broken_data_ptr =
1877  boost::make_shared<std::vector<BrokenBaseSideData>>();
1878  // Note: EleOnSide, i.e. uses on domain projected skeleton rule
1879  auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
1880  mField, elementVolumeName, SPACE_DIM, Sev::noisy);
1881  op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
1882  boost::make_shared<CGGUserPolynomialBase>();
1883  CHKERR
1884  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
1885  op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
1886  materialH1Positions, frontAdjEdges);
1887  op_loop_domain_side->getOpPtrVector().push_back(
1888  new OpGetBrokenBaseSideData<SideEleOp>(piolaStress, broken_data_ptr));
1889 
1890  op_loop_skeleton_side->getOpPtrVector().push_back(op_loop_domain_side);
1891  using OpC = FormsIntegrators<BdyEleOp>::Assembly<A>::BiLinearForm<
1892  GAUSS>::OpBrokenSpaceConstrain<SPACE_DIM>;
1893  op_loop_skeleton_side->getOpPtrVector().push_back(
1894  new OpC(hybridSpatialDisp, broken_data_ptr, 1., true, false));
1895 
1896  pip.push_back(op_loop_skeleton_side);
1897 
1898  MoFEMFunctionReturn(0);
1899  };
1900 
1901  auto set_contact = [&](auto &pip) {
1902  MoFEMFunctionBegin;
1903 
1904  using BoundaryEle =
1905  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::BoundaryEle;
1906  using EleOnSide =
1907  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
1908  using SideEleOp = EleOnSide::UserDataOperator;
1909  using BdyEleOp = BoundaryEle::UserDataOperator;
1910 
1911  // First: Iterate over skeleton FEs adjacent to Domain FEs
1912  // Note: BoundaryEle, i.e. uses skeleton interation rule
1913  auto op_loop_skeleton_side = new OpLoopSide<BoundaryEle>(
1914  mField, contactElement, SPACE_DIM - 1, Sev::noisy);
1915 
1916  op_loop_skeleton_side->getSideFEPtr()->getRuleHook = rule_contact;
1917  op_loop_skeleton_side->getSideFEPtr()->setRuleHook = set_rule_contact;
1918  CHKERR EshelbianPlasticity::
1919  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
1920  op_loop_skeleton_side->getOpPtrVector(), {L2},
1921  materialH1Positions, frontAdjEdges);
1922 
1923  // Second: Iterate over domain FEs adjacent to skelton, particularly
1924  // one domain element.
1925  auto broken_data_ptr =
1926  boost::make_shared<std::vector<BrokenBaseSideData>>();
1927 
1928  // Data storing contact fields
1929  auto contact_common_data_ptr =
1930  boost::make_shared<ContactOps::CommonData>();
1931 
1932  auto add_ops_domain_side = [&](auto &pip) {
1933  MoFEMFunctionBegin;
1934  // Note: EleOnSide, i.e. uses on domain projected skeleton rule
1935  auto op_loop_domain_side = new OpBrokenLoopSide<EleOnSide>(
1936  mField, elementVolumeName, SPACE_DIM, Sev::noisy);
1937  op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
1938  boost::make_shared<CGGUserPolynomialBase>();
1939  CHKERR
1940  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
1941  op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
1942  materialH1Positions, frontAdjEdges);
1943  op_loop_domain_side->getOpPtrVector().push_back(
1944  new OpGetBrokenBaseSideData<SideEleOp>(piolaStress,
1945  broken_data_ptr));
1946  op_loop_domain_side->getOpPtrVector().push_back(
1947  new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
1948  piolaStress, contact_common_data_ptr->contactTractionPtr()));
1949  pip.push_back(op_loop_domain_side);
1950  MoFEMFunctionReturn(0);
1951  };
1952 
1953  auto add_ops_contact_lhs = [&](auto &pip) {
1954  MoFEMFunctionBegin;
1955  pip.push_back(new OpCalculateVectorFieldValues<3>(
1956  contactDisp, contact_common_data_ptr->contactDispPtr()));
1957  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
1958  pip.push_back(
1959  new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
1960  pip.push_back(new OpTreeSearch(
1961  contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
1962  get_range_from_block(mField, "CONTACT", SPACE_DIM - 1), nullptr,
1963  nullptr));
1964 
1965  // get body id and SDF range
1966  auto contact_sfd_map_range_ptr =
1967  boost::make_shared<std::map<int, Range>>(
1968  get_body_range("CONTACT_SDF", SPACE_DIM - 1));
1969 
1970  pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dU(
1971  contactDisp, contactDisp, contact_common_data_ptr, contactTreeRhs,
1972  contact_sfd_map_range_ptr));
1973  pip.push_back(
1974  new EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dP<A, GAUSS>(
1975  contactDisp, broken_data_ptr, contact_common_data_ptr,
1976  contactTreeRhs, contact_sfd_map_range_ptr));
1977  pip.push_back(
1978  new EshelbianPlasticity::OpConstrainBoundaryHDivLhs_dU<A, GAUSS>(
1979  broken_data_ptr, contactDisp, contact_common_data_ptr,
1980  contactTreeRhs));
1981 
1982  MoFEMFunctionReturn(0);
1983  };
1984 
1985  // push ops to face/side pipeline
1986  CHKERR add_ops_domain_side(op_loop_skeleton_side->getOpPtrVector());
1987  CHKERR add_ops_contact_lhs(op_loop_skeleton_side->getOpPtrVector());
1988 
1989  // Add skeleton to domain pipeline
1990  pip.push_back(op_loop_skeleton_side);
1991 
1992  MoFEMFunctionReturn(0);
1993  };
1994 
1995  CHKERR set_hybridisation(fe_lhs->getOpPtrVector());
1996  CHKERR set_contact(fe_lhs->getOpPtrVector());
1997  }
1998 
1999  if (add_material) {
2000  }
2001 
2002  MoFEMFunctionReturn(0);
2003 }
2004 
2005 MoFEMErrorCode EshelbianCore::setFaceElementOps(
2006  const bool add_elastic, const bool add_material,
2007  boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
2008  boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
2009  MoFEMFunctionBegin;
2010 
2011  fe_rhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
2012  fe_lhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
2013 
2014  // set integration rule
2015  fe_rhs->getRuleHook = [](int, int, int p) { return 2 * (p + 1); };
2016  fe_lhs->getRuleHook = [](int, int, int p) { return 2 * (p + 1); };
2017 
2018  CHKERR
2019  EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2020  fe_rhs->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
2021  CHKERR
2022  EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2023  fe_lhs->getOpPtrVector(), {L2}, materialH1Positions, frontAdjEdges);
2024 
2025  if (add_elastic) {
2026 
2027  auto get_broken_op_side = [this](auto &pip) {
2028  using EleOnSide =
2029  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
2030  using SideEleOp = EleOnSide::UserDataOperator;
2031  // Iterate over domain FEs adjacent to boundary.
2032  auto broken_data_ptr =
2033  boost::make_shared<std::vector<BrokenBaseSideData>>();
2034  // Note: EleOnSide, i.e. uses on domain projected skeleton rule
2035  auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
2036  mField, elementVolumeName, SPACE_DIM, Sev::noisy);
2037  op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2038  boost::make_shared<CGGUserPolynomialBase>();
2039  CHKERR
2040  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2041  op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2042  materialH1Positions, frontAdjEdges);
2043  op_loop_domain_side->getOpPtrVector().push_back(
2044  new OpGetBrokenBaseSideData<SideEleOp>(piolaStress, broken_data_ptr));
2045  pip.push_back(op_loop_domain_side);
2046  return broken_data_ptr;
2047  };
2048 
2049  auto broken_data_ptr = get_broken_op_side(fe_rhs->getOpPtrVector());
2050 
2051  fe_rhs->getOpPtrVector().push_back(
2052  new OpDispBc(broken_data_ptr, bcSpatialDispVecPtr,
2053  {
2054 
2055  boost::make_shared<TimeScale>("disp_history.txt")
2056 
2057  }));
2058  fe_rhs->getOpPtrVector().push_back(new OpRotationBc(
2059  broken_data_ptr, bcSpatialRotationVecPtr,
2060 
2061  {
2062 
2063  boost::make_shared<TimeScale>("rotation_history.txt")
2064 
2065  }));
2066 
2067  fe_rhs->getOpPtrVector().push_back(
2068  new OpBrokenTractionBc(hybridSpatialDisp, bcSpatialTraction));
2069  }
2070 
2071  MoFEMFunctionReturn(0);
2072 }
2073 
2074 MoFEMErrorCode EshelbianCore::setContactElementRhsOps(
2075 
2076  boost::shared_ptr<ContactTree> &fe_contact_tree
2077 
2078 ) {
2079  MoFEMFunctionBegin;
2080 
2081  /** Contact requires that body is marked */
2082  auto get_body_range = [this](auto name, int dim) {
2083  std::map<int, Range> map;
2084 
2085  for (auto m_ptr :
2086  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
2087 
2088  (boost::format("%s(.*)") % name).str()
2089 
2090  ))
2091 
2092  ) {
2093  Range ents;
2094  CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
2095  dim, ents, true),
2096  "by dim");
2097  map[m_ptr->getMeshsetId()] = ents;
2098  }
2099 
2100  return map;
2101  };
2102 
2103  auto get_map_skin = [this](auto &&map) {
2104  ParallelComm *pcomm =
2105  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
2106 
2107  Skinner skin(&mField.get_moab());
2108  for (auto &m : map) {
2109  Range skin_faces;
2110  CHKERR skin.find_skin(0, m.second, false, skin_faces);
2111  CHK_MOAB_THROW(pcomm->filter_pstatus(skin_faces,
2112  PSTATUS_SHARED | PSTATUS_MULTISHARED,
2113  PSTATUS_NOT, -1, nullptr),
2114  "filter");
2115  m.second.swap(skin_faces);
2116  }
2117  return map;
2118  };
2119 
2120  /* The above code is written in C++ and it appears to be defining and using
2121  various operations on boundary elements and side elements. */
2122  using BoundaryEle = FaceElementForcesAndSourcesCore;
2123  using BoundaryEleOp = BoundaryEle::UserDataOperator;
2124 
2125  auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
2126 
2127  auto calcs_side_traction = [&](auto &pip) {
2128  MoFEMFunctionBegin;
2129  using EleOnSide =
2130  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
2131  using SideEleOp = EleOnSide::UserDataOperator;
2132  auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
2133  mField, elementVolumeName, SPACE_DIM, Sev::noisy);
2134  op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2135  boost::make_shared<CGGUserPolynomialBase>();
2136  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2137  op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2138  materialH1Positions, frontAdjEdges);
2139  op_loop_domain_side->getOpPtrVector().push_back(
2140  new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
2141  piolaStress, contact_common_data_ptr->contactTractionPtr()));
2142  pip.push_back(op_loop_domain_side);
2143  MoFEMFunctionReturn(0);
2144  };
2145 
2146  auto add_contact_three = [&]() {
2147  MoFEMFunctionBegin;
2148  auto tree_moab_ptr = boost::make_shared<moab::Core>();
2149  fe_contact_tree = boost::make_shared<ContactTree>(
2150  mField, tree_moab_ptr, spaceOrder,
2151  get_map_skin(get_body_range("BODY", 3)));
2152  fe_contact_tree->getOpPtrVector().push_back(
2153  new OpCalculateVectorFieldValues<3>(
2154  contactDisp, contact_common_data_ptr->contactDispPtr()));
2155  CHKERR calcs_side_traction(fe_contact_tree->getOpPtrVector());
2156  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
2157  fe_contact_tree->getOpPtrVector().push_back(
2158  new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
2159  fe_contact_tree->getOpPtrVector().push_back(
2160  new OpMoveNode(fe_contact_tree, contact_common_data_ptr, u_h1_ptr));
2161  MoFEMFunctionReturn(0);
2162  };
2163 
2164  CHKERR add_contact_three();
2165 
2166  MoFEMFunctionReturn(0);
2167 }
2168 
2169 MoFEMErrorCode EshelbianCore::setElasticElementOps(const int tag) {
2170  MoFEMFunctionBegin;
2171 
2172  // Add contact operators. Note that only for rhs. THe lhs is assembled with
2173  // volume element, to enable schur complement evaluation.
2174  CHKERR setContactElementRhsOps(contactTreeRhs);
2175 
2176  CHKERR setVolumeElementOps(tag, true, false, elasticFeRhs, elasticFeLhs);
2177  CHKERR setFaceElementOps(true, false, elasticBcRhs, elasticBcLhs);
2178 
2179  auto adj_cache =
2180  boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
2181 
2182  auto get_op_contact_bc = [&]() {
2183  using SideEle = FaceElementForcesAndSourcesCore;
2184  auto op_loop_side = new OpLoopSide<SideEle>(
2185  mField, contactElement, SPACE_DIM - 1, Sev::noisy, adj_cache);
2186  return op_loop_side;
2187  };
2188 
2189  MoFEMFunctionReturn(0);
2190 }
2191 
2192 MoFEMErrorCode EshelbianCore::setElasticElementToTs(DM dm) {
2193  MoFEMFunctionBegin;
2194  boost::shared_ptr<FEMethod> null;
2195 
2196  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
2197 
2198  CHKERR DMMoFEMTSSetI2Function(dm, elementVolumeName, elasticFeRhs, null,
2199  null);
2200  CHKERR DMMoFEMTSSetI2Function(dm, naturalBcElement, elasticBcRhs, null,
2201  null);
2202  CHKERR DMMoFEMTSSetI2Jacobian(dm, elementVolumeName, elasticFeLhs, null,
2203  null);
2204  // CHKERR DMMoFEMTSSetI2Jacobian(dm, naturalBcElement, elasticBcLhs, null,
2205  // null);
2206 
2207  } else {
2208  CHKERR DMMoFEMTSSetIFunction(dm, elementVolumeName, elasticFeRhs, null,
2209  null);
2210  CHKERR DMMoFEMTSSetIFunction(dm, naturalBcElement, elasticBcRhs, null,
2211  null);
2212  CHKERR DMMoFEMTSSetIJacobian(dm, elementVolumeName, elasticFeLhs, null,
2213  null);
2214  // CHKERR DMMoFEMTSSetIJacobian(dm, naturalBcElement, elasticBcLhs, null,
2215  // null);
2216  }
2217 
2218  MoFEMFunctionReturn(0);
2219 }
2220 
2221 #include "impl/EshelbianMonitor.cpp"
2222 #include "impl/SetUpSchurImpl.cpp"
2223 #include "impl/TSElasticPostStep.cpp"
2224 
2225 MoFEMErrorCode EshelbianCore::solveElastic(TS ts, Vec x) {
2226  MoFEMFunctionBegin;
2227 
2228 #ifdef PYTHON_SDF
2229 
2230  boost::shared_ptr<ContactOps::SDFPython> sdf_python_ptr;
2231 
2232  auto file_exists = [](std::string myfile) {
2233  std::ifstream file(myfile.c_str());
2234  if (file) {
2235  return true;
2236  }
2237  return false;
2238  };
2239 
2240  char sdf_file_name[255] = "sdf.py";
2241  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-sdf_file",
2242  sdf_file_name, 255, PETSC_NULL);
2243 
2244  if (file_exists(sdf_file_name)) {
2245  MOFEM_LOG("EP", Sev::inform) << sdf_file_name << " file found";
2246  sdf_python_ptr = boost::make_shared<ContactOps::SDFPython>();
2247  CHKERR sdf_python_ptr->sdfInit(sdf_file_name);
2248  ContactOps::sdfPythonWeakPtr = sdf_python_ptr;
2249  } else {
2250  MOFEM_LOG("EP", Sev::warning) << sdf_file_name << " file NOT found";
2251  }
2252 #else
2253 #endif
2254 
2255  boost::shared_ptr<TsCtx> ts_ctx;
2256  CHKERR DMMoFEMGetTsCtx(dmElastic, ts_ctx);
2257  CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx.get(), PETSC_NULL);
2258 
2259  boost::shared_ptr<FEMethod> monitor_ptr(new EshelbianMonitor(*this));
2260  ts_ctx->getLoopsMonitor().push_back(
2261  TsCtx::PairNameFEMethodPtr(elementVolumeName, monitor_ptr));
2262 
2263  CHKERR TSAppendOptionsPrefix(ts, "elastic_");
2264  CHKERR TSSetFromOptions(ts);
2265 
2266  CHKERR TSSetDM(ts, dmElastic);
2267 
2268  SNES snes;
2269  CHKERR TSGetSNES(ts, &snes);
2270 
2271  PetscViewerAndFormat *vf;
2272  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
2273  PETSC_VIEWER_DEFAULT, &vf);
2274  CHKERR SNESMonitorSet(
2275  snes,
2276  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
2277  vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
2278 
2279  PetscSection section;
2280  CHKERR DMGetSection(dmElastic, &section);
2281  int num_fields;
2282  CHKERR PetscSectionGetNumFields(section, &num_fields);
2283  for (int ff = 0; ff != num_fields; ff++) {
2284  const char *field_name;
2285  CHKERR PetscSectionGetFieldName(section, ff, &field_name);
2286  MOFEM_LOG_C("EP", Sev::inform, "Field %d name %s", ff, field_name);
2287  }
2288 
2289  CHKERR DMoFEMMeshToLocalVector(dmElastic, x, INSERT_VALUES, SCATTER_FORWARD);
2290  CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
2291  CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
2292 
2293  // Adding field split solver
2294  boost::shared_ptr<SetUpSchur> schur_ptr;
2295  if constexpr (A == AssemblyType::BLOCK_SCHUR) {
2296  schur_ptr = SetUpSchur::createSetUpSchur(mField, &*this);
2297  CHKERR schur_ptr->setUp(ts);
2298  }
2299 
2300  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
2301  Vec xx;
2302  CHKERR VecDuplicate(x, &xx);
2303  CHKERR VecZeroEntries(xx);
2304  CHKERR TS2SetSolution(ts, x, xx);
2305  CHKERR VecDestroy(&xx);
2306  } else {
2307  CHKERR TSSetSolution(ts, x);
2308  }
2309 
2310  TetPolynomialBase::switchCacheBaseOn<HDIV>(
2311  {elasticFeLhs.get(), elasticFeRhs.get()});
2312  CHKERR TSSetUp(ts);
2313  CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
2314  CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
2315  CHKERR TSElasticPostStep::postStepInitialise(this);
2316  CHKERR TSSolve(ts, PETSC_NULL);
2317  CHKERR TSElasticPostStep::postStepDestroy();
2318  TetPolynomialBase::switchCacheBaseOff<HDIV>(
2319  {elasticFeLhs.get(), elasticFeRhs.get()});
2320 
2321  // CHKERR TSGetSNES(ts, &snes);
2322  int lin_solver_iterations;
2323  CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
2324  MOFEM_LOG("EP", Sev::inform)
2325  << "Number of linear solver iterations " << lin_solver_iterations;
2326 
2327  PetscBool test_cook_flg = PETSC_FALSE;
2328  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_cook", &test_cook_flg,
2329  PETSC_NULL);
2330  if (test_cook_flg) {
2331  constexpr int expected_lin_solver_iterations = 11;
2332  if (lin_solver_iterations > expected_lin_solver_iterations)
2333  SETERRQ2(
2334  PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
2335  "Expected number of iterations is different than expected %d > %d",
2336  lin_solver_iterations, expected_lin_solver_iterations);
2337  }
2338 
2339  CHKERR gettingNorms();
2340 
2341  MoFEMFunctionReturn(0);
2342 }
2343 
2344 MoFEMErrorCode EshelbianCore::postProcessResults(const int tag,
2345  const std::string file) {
2346  MoFEMFunctionBegin;
2347 
2348  if (!dataAtPts) {
2349  dataAtPts =
2350  boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
2351  }
2352 
2353  auto post_proc_mesh = boost::make_shared<moab::Core>();
2354  auto post_proc_begin =
2355  PostProcBrokenMeshInMoabBaseBegin(mField, post_proc_mesh);
2356  auto post_proc_end = PostProcBrokenMeshInMoabBaseEnd(mField, post_proc_mesh);
2357 
2358  auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
2359 
2360  auto get_post_proc = [&](auto sense) {
2361  using FaceEle = FaceElementForcesAndSourcesCore;
2362 
2363  auto post_proc_ptr =
2364  boost::make_shared<PostProcBrokenMeshInMoabBaseCont<FaceEle>>(
2365  mField, post_proc_mesh);
2366 
2367  EshelbianPlasticity::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
2368  post_proc_ptr->getOpPtrVector(), {L2}, materialH1Positions,
2369  frontAdjEdges);
2370 
2371  auto domain_ops = [&](auto &fe, int sense) {
2372  MoFEMFunctionBegin;
2373  fe.getUserPolynomialBase() =
2374  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
2375  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2376  fe.getOpPtrVector(), {HDIV, H1, L2}, materialH1Positions,
2377  frontAdjEdges);
2378  fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
2379  piolaStress, dataAtPts->getApproxPAtPts()));
2380  fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
2381  bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
2382  if (noStretch) {
2383  fe.getOpPtrVector().push_back(
2384  physicalEquations->returnOpCalculateStretchFromStress(
2385  dataAtPts, physicalEquations));
2386  } else {
2387  fe.getOpPtrVector().push_back(
2388  new OpCalculateTensor2SymmetricFieldValues<3>(
2389  stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
2390  }
2391 
2392  fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
2393  rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
2394  fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
2395  spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
2396  fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
2397  spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
2398  fe.getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
2399  spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
2400  // evaluate derived quantities
2401  fe.getOpPtrVector().push_back(
2402  new OpCalculateRotationAndSpatialGradient(dataAtPts));
2403 
2404  // evaluate integration points
2405  fe.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
2406  tag, true, false, dataAtPts, physicalEquations));
2407  if (auto op =
2408  physicalEquations->returnOpCalculateEnergy(dataAtPts, nullptr)) {
2409  fe.getOpPtrVector().push_back(op);
2410  fe.getOpPtrVector().push_back(new OpCalculateEshelbyStress(dataAtPts));
2411  }
2412 
2413  // contact
2414  fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
2415  spatialL2Disp, contact_common_data_ptr->contactDispPtr()));
2416 
2417  // post-proc
2418  fe.getOpPtrVector().push_back(new OpPostProcDataStructure(
2419  post_proc_ptr->getPostProcMesh(), post_proc_ptr->getMapGaussPts(),
2420  dataAtPts, sense));
2421 
2422  MoFEMFunctionReturn(0);
2423  };
2424 
2425  post_proc_ptr->getOpPtrVector().push_back(
2426  new OpCalculateVectorFieldValues<3>(contactDisp,
2427  dataAtPts->getContactL2AtPts()));
2428  auto X_h1_ptr = boost::make_shared<MatrixDouble>();
2429  // H1 material positions
2430  post_proc_ptr->getOpPtrVector().push_back(
2431  new OpCalculateVectorFieldValues<3>(materialH1Positions,
2432  dataAtPts->getLargeXH1AtPts()));
2433 
2434  // domain
2435  auto op_loop_side = new OpLoopSide<VolumeElementForcesAndSourcesCoreOnSide>(
2436  mField, elementVolumeName, SPACE_DIM);
2437  domain_ops(*(op_loop_side->getSideFEPtr()), sense);
2438  post_proc_ptr->getOpPtrVector().push_back(op_loop_side);
2439 
2440  return post_proc_ptr;
2441  };
2442 
2443  auto post_proc_ptr = get_post_proc(1);
2444  auto post_proc_negative_sense_ptr = get_post_proc(-1);
2445 
2446  // contact
2447  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
2448  auto lambda_h1_ptr = boost::make_shared<MatrixDouble>();
2449  post_proc_ptr->getOpPtrVector().push_back(
2450  new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
2451 
2452  using BoundaryEle = FaceElementForcesAndSourcesCore;
2453  using BoundaryEleOp = BoundaryEle::UserDataOperator;
2454 
2455  auto calcs_side_traction = [&](auto &pip) {
2456  MoFEMFunctionBegin;
2457  using EleOnSide =
2458  PipelineManager::ElementsAndOpsByDim<SPACE_DIM>::FaceSideEle;
2459  using SideEleOp = EleOnSide::UserDataOperator;
2460  auto op_loop_domain_side = new OpLoopSide<EleOnSide>(
2461  mField, elementVolumeName, SPACE_DIM, Sev::noisy);
2462  op_loop_domain_side->getSideFEPtr()->getUserPolynomialBase() =
2463  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
2464  EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2465  op_loop_domain_side->getOpPtrVector(), {HDIV, H1, L2},
2466  materialH1Positions, frontAdjEdges);
2467  op_loop_domain_side->getOpPtrVector().push_back(
2468  new OpCalculateHVecTensorTrace<SPACE_DIM, SideEleOp>(
2469  piolaStress, contact_common_data_ptr->contactTractionPtr()));
2470  pip.push_back(op_loop_domain_side);
2471  MoFEMFunctionReturn(0);
2472  };
2473 
2474  CHKERR calcs_side_traction(post_proc_ptr->getOpPtrVector());
2475 
2476  post_proc_ptr->getOpPtrVector().push_back(new OpTreeSearch(
2477  contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
2478  get_range_from_block(mField, "CONTACT", SPACE_DIM - 1),
2479  &post_proc_ptr->getPostProcMesh(), &post_proc_ptr->getMapGaussPts()));
2480 
2481  CHKERR DMoFEMLoopFiniteElements(dM, contactElement, contactTreeRhs);
2482 
2483  CHKERR DMoFEMPreProcessFiniteElements(dM, post_proc_begin.getFEMethod());
2484  CHKERR DMoFEMLoopFiniteElements(dM, skinElement, post_proc_ptr);
2485  CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(dM, crackElement, post_proc_ptr,
2486  0, mField.get_comm_size());
2487  CHKERR DMoFEMLoopFiniteElementsUpAndLowRank(dM, crackElement,
2488  post_proc_negative_sense_ptr, 0,
2489  mField.get_comm_size());
2490  CHKERR DMoFEMPostProcessFiniteElements(dM, post_proc_end.getFEMethod());
2491 
2492  CHKERR post_proc_end.writeFile(file.c_str());
2493  MoFEMFunctionReturn(0);
2494 }
2495 
2496 //! [Getting norms]
2497 MoFEMErrorCode EshelbianCore::gettingNorms() {
2498  MoFEMFunctionBegin;
2499 
2500  auto post_proc_norm_fe =
2501  boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
2502 
2503  auto post_proc_norm_rule_hook = [](int, int, int p) -> int {
2504  return 2 * (p);
2505  };
2506  post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
2507 
2508  post_proc_norm_fe->getUserPolynomialBase() =
2509  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
2510 
2511  CHKERR EshelbianPlasticity::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
2512  post_proc_norm_fe->getOpPtrVector(), {L2, H1, HDIV}, materialH1Positions,
2513  frontAdjEdges);
2514 
2515  enum NORMS { U_NORM_L2 = 0, U_NORM_H1, PIOLA_NORM, U_ERROR_L2, LAST_NORM };
2516  auto norms_vec =
2517  createVectorMPI(mField.get_comm(), LAST_NORM, PETSC_DETERMINE);
2518  CHKERR VecZeroEntries(norms_vec);
2519 
2520  auto u_l2_ptr = boost::make_shared<MatrixDouble>();
2521  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
2522  post_proc_norm_fe->getOpPtrVector().push_back(
2523  new OpCalculateVectorFieldValues<SPACE_DIM>(spatialL2Disp, u_l2_ptr));
2524  post_proc_norm_fe->getOpPtrVector().push_back(
2525  new OpCalculateVectorFieldValues<SPACE_DIM>(spatialH1Disp, u_h1_ptr));
2526  post_proc_norm_fe->getOpPtrVector().push_back(
2527  new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_NORM_L2));
2528  post_proc_norm_fe->getOpPtrVector().push_back(
2529  new OpCalcNormL2Tensor1<SPACE_DIM>(u_h1_ptr, norms_vec, U_NORM_H1));
2530  post_proc_norm_fe->getOpPtrVector().push_back(
2531  new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_ERROR_L2,
2532  u_h1_ptr));
2533 
2534  auto piola_ptr = boost::make_shared<MatrixDouble>();
2535  post_proc_norm_fe->getOpPtrVector().push_back(
2536  new OpCalculateHVecTensorField<3, 3>(piolaStress, piola_ptr));
2537  post_proc_norm_fe->getOpPtrVector().push_back(
2538  new OpCalcNormL2Tensor2<3, 3>(piola_ptr, norms_vec, PIOLA_NORM));
2539 
2540  TetPolynomialBase::switchCacheBaseOn<HDIV>({post_proc_norm_fe.get()});
2541  CHKERR mField.loop_finite_elements("ELASTIC_PROBLEM", elementVolumeName,
2542  *post_proc_norm_fe);
2543  TetPolynomialBase::switchCacheBaseOff<HDIV>({post_proc_norm_fe.get()});
2544 
2545  CHKERR VecAssemblyBegin(norms_vec);
2546  CHKERR VecAssemblyEnd(norms_vec);
2547  const double *norms;
2548  CHKERR VecGetArrayRead(norms_vec, &norms);
2549  MOFEM_LOG("EP", Sev::inform) << "norm_u: " << std::sqrt(norms[U_NORM_L2]);
2550  MOFEM_LOG("EP", Sev::inform) << "norm_u_h1: " << std::sqrt(norms[U_NORM_H1]);
2551  MOFEM_LOG("EP", Sev::inform)
2552  << "norm_error_u_l2: " << std::sqrt(norms[U_ERROR_L2]);
2553  MOFEM_LOG("EP", Sev::inform)
2554  << "norm_piola: " << std::sqrt(norms[PIOLA_NORM]);
2555  CHKERR VecRestoreArrayRead(norms_vec, &norms);
2556 
2557  MoFEMFunctionReturn(0);
2558 }
2559 //! [Getting norms]
2560 
2561 MoFEMErrorCode EshelbianCore::getSpatialDispBc() {
2562  MoFEMFunctionBegin;
2563 
2564  auto bc_mng = mField.getInterface<BcManager>();
2565  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
2566  "", piolaStress, false, false);
2567 
2568  bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
2569 
2570  for (auto bc : bc_mng->getBcMapByBlockName()) {
2571  if (auto disp_bc = bc.second->dispBcPtr) {
2572 
2573  MOFEM_LOG("EP", Sev::noisy) << *disp_bc;
2574 
2575  std::vector<double> block_attributes(6, 0.);
2576  if (disp_bc->data.flag1 == 1) {
2577  block_attributes[0] = disp_bc->data.value1;
2578  block_attributes[3] = 1;
2579  }
2580  if (disp_bc->data.flag2 == 1) {
2581  block_attributes[1] = disp_bc->data.value2;
2582  block_attributes[4] = 1;
2583  }
2584  if (disp_bc->data.flag3 == 1) {
2585  block_attributes[2] = disp_bc->data.value3;
2586  block_attributes[5] = 1;
2587  }
2588  auto faces = bc.second->bcEnts.subset_by_dimension(2);
2589  bcSpatialDispVecPtr->emplace_back(bc.first, block_attributes, faces);
2590  }
2591  }
2592 
2593  // old way of naming blocksets for displacement BCs
2594  CHKERR getBc(bcSpatialDispVecPtr, "SPATIAL_DISP_BC", 6);
2595 
2596  MoFEMFunctionReturn(0);
2597 }
2598 
2599 MoFEMErrorCode EshelbianCore::getSpatialTractionBc() {
2600  MoFEMFunctionBegin;
2601 
2602  auto bc_mng = mField.getInterface<BcManager>();
2603  CHKERR bc_mng->pushMarkDOFsOnEntities<ForceCubitBcData>("", piolaStress,
2604  false, false);
2605 
2606  bcSpatialTraction = boost::make_shared<TractionBcVec>();
2607 
2608  for (auto bc : bc_mng->getBcMapByBlockName()) {
2609  if (auto force_bc = bc.second->forceBcPtr) {
2610  std::vector<double> block_attributes(6, 0.);
2611  block_attributes[0] = -force_bc->data.value3 * force_bc->data.value1;
2612  block_attributes[3] = 1;
2613  block_attributes[1] = -force_bc->data.value4 * force_bc->data.value1;
2614  block_attributes[4] = 1;
2615  block_attributes[2] = -force_bc->data.value5 * force_bc->data.value1;
2616  block_attributes[5] = 1;
2617  auto faces = bc.second->bcEnts.subset_by_dimension(2);
2618  bcSpatialTraction->emplace_back(bc.first, block_attributes, faces);
2619  }
2620  }
2621 
2622  CHKERR getBc(bcSpatialTraction, "SPATIAL_TRACTION_BC", 6);
2623 
2624  MoFEMFunctionReturn(0);
2625 }
2626 
2627 } // namespace EshelbianPlasticity
2628 
2629 #include <impl/EshelbianContact.cpp>
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::CoreInterface::modify_finite_element_adjacency_table
virtual MoFEMErrorCode modify_finite_element_adjacency_table(const std::string &fe_name, const EntityType type, ElementAdjacencyFunct function)=0
modify finite element table, only for advanced user
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:589
EshelbianPlasticity::EshelbianCore::f_log_e
static double f_log_e(const double v)
Definition: EshelbianPlasticity.hpp:910
EshelbianPlasticity::EshelbianCore::f_linear
static double f_linear(const double v)
Definition: EshelbianPlasticity.hpp:935
EshelbianPlasticity::OpSpatialPhysical_du_dP
Definition: EshelbianPlasticity.hpp:712
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
EshelbianPlasticity::BcRot::BcRot
BcRot(std::string name, std::vector< double > &attr, Range &faces)
Definition: EshelbianPlasticity.cpp:1254
MoFEM::VolumeElementForcesAndSourcesCore::VolumeElementForcesAndSourcesCore
VolumeElementForcesAndSourcesCore(Interface &m_field, const EntityType type=MBTET)
Definition: VolumeElementForcesAndSourcesCore.cpp:9
EshelbianPlasticity::OpMoveNode
Definition: EshelbianContact.hpp:221
EshelbianPlasticity::OpCalculateRotationAndSpatialGradient
Definition: EshelbianPlasticity.hpp:557
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
MoFEM::CoreInterface::loop_finite_elements
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
EshelbianPlasticity::LINEAR
@ LINEAR
Definition: EshelbianPlasticity.hpp:44
MoFEM::CoreInterface::loop_dofs
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.
EshelbianPlasticity::CGGUserPolynomialBase
Definition: EshelbianPlasticity.hpp:21
EshelbianPlasticity::VolRule
Set integration rule on element.
Definition: EshelbianPlasticity.cpp:1353
H1
@ H1
continuous field
Definition: definitions.h:85
EshelbianPlasticity::EshelbianCore::alphaW
double alphaW
Definition: EshelbianPlasticity.hpp:997
EshelbianPlasticity::EshelbianCore::bcSpatialRotationVecPtr
boost::shared_ptr< BcRotVec > bcSpatialRotationVecPtr
Definition: EshelbianPlasticity.hpp:1006
EshelbianPlasticity::EshelbianCore::dd_f
static boost::function< double(const double)> dd_f
Definition: EshelbianPlasticity.hpp:905
TSElasticPostStep.cpp
EshelbianPlasticity::SetIntegrationAtFrontVolume
Definition: EshelbianPlasticity.cpp:103
EshelbianPlasticity::EshelbianCore::hybridMaterialDisp
const std::string hybridMaterialDisp
Definition: EshelbianPlasticity.hpp:975
EshelbianPlasticity::EshelbianCore::inv_f_log
static double inv_f_log(const double v)
Definition: EshelbianPlasticity.hpp:925
EshelbianPlasticity::EshelbianCore::getTractionFreeBc
MoFEMErrorCode getTractionFreeBc(const EntityHandle meshset, boost::shared_ptr< TractionFreeBc > &bc_ptr, const std::string contact_set_name)
Remove all, but entities where kinematic constrains are applied.
Definition: EshelbianPlasticity.cpp:1282
EntityHandle
EshelbianPlasticity::EshelbianCore::piolaStress
const std::string piolaStress
Definition: EshelbianPlasticity.hpp:968
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
EshelbianPlasticity::DataAtIntegrationPts
Definition: EshelbianPlasticity.hpp:55
EshelbianPlasticity::EshelbianCore::spatialH1Disp
const std::string spatialH1Disp
Definition: EshelbianPlasticity.hpp:972
EshelbianPlasticity::EshelbianCore::addDMs
MoFEMErrorCode addDMs(const BitRefLevel bit=BitRefLevel().set(0), const EntityHandle meshset=0)
Definition: EshelbianPlasticity.cpp:1101
MoFEM::EntPolynomialBaseCtx
Class used to pass element data to calculate base functions on tet,triangle,edge.
Definition: EntPolynomialBaseCtx.hpp:22
EshelbianPlasticity::EshelbianCore::setBaseVolumeElementOps
MoFEMErrorCode setBaseVolumeElementOps(const int tag, const bool do_rhs, const bool do_lhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > &fe)
Definition: EshelbianPlasticity.cpp:1488
EshelbianPlasticity::OpSpatialConsistency_dBubble_dP
Definition: EshelbianPlasticity.hpp:806
EshelbianPlasticity::EshelbianCore::projectGeometry
MoFEMErrorCode projectGeometry(const EntityHandle meshset=0)
Definition: EshelbianPlasticity.cpp:707
EshelbianPlasticity::EshelbianCore::exponentBase
static double exponentBase
Definition: EshelbianPlasticity.hpp:902
EshelbianPlasticity::EshelbianCore::inv_dd_f_linear
static double inv_dd_f_linear(const double v)
Definition: EshelbianPlasticity.hpp:941
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
EleOnSide
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::FaceSideEle EleOnSide
Definition: scalar_check_approximation.cpp:27
MoFEM::BaseFunction::DofsSideMap
multi_index_container< DofsSideMapData, indexed_by< ordered_non_unique< tag< TypeSide_mi_tag >, composite_key< DofsSideMapData, member< DofsSideMapData, EntityType, &DofsSideMapData::type >, member< DofsSideMapData, int, &DofsSideMapData::side > >>, ordered_unique< tag< EntDofIdx_mi_tag >, member< DofsSideMapData, int, &DofsSideMapData::dof > > > > DofsSideMap
Map entity stype and side to element/entity dof index.
Definition: BaseFunction.hpp:73
EshelbianPlasticity::SetIntegrationAtFrontVolume::mapRefCoords
static std::map< long int, MatrixDouble > mapRefCoords
Definition: EshelbianPlasticity.cpp:336
MoFEM::PipelineManager::ElementsAndOpsByDim
Definition: PipelineManager.hpp:38
MoFEM::OpCalculateTensor2SymmetricFieldValuesDot
Calculate symmetric tensor field rates ant integratio pts.
Definition: UserDataOperators.hpp:1071
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
EshelbianPlasticity::EshelbianCore::dataAtPts
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
Definition: EshelbianPlasticity.hpp:954
NOBASE
@ NOBASE
Definition: definitions.h:59
EshelbianPlasticity::EshelbianCore::spaceH1Order
int spaceH1Order
Definition: EshelbianPlasticity.hpp:994
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:456
EshelbianPlasticity::EshelbianCore::materialSkeletonFaces
boost::shared_ptr< Range > materialSkeletonFaces
Definition: EshelbianPlasticity.hpp:1137
MoFEM::TsCtx::getLoopsMonitor
FEMethodsSequence & getLoopsMonitor()
Get the loops to do Monitor object.
Definition: TsCtx.hpp:98
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
EshelbianPlasticity::EshelbianCore::bitAdjParentMask
BitRefLevel bitAdjParentMask
bit ref level for parent parent
Definition: EshelbianPlasticity.hpp:1143
TSElasticPostStep::postStepFun
static MoFEMErrorCode postStepFun(TS ts)
Definition: TSElasticPostStep.cpp:145
MoFEM::OpCalculateHTensorTensorField
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:2631
EshelbianPlasticity::EshelbianCore::alphaRho
double alphaRho
Definition: EshelbianPlasticity.hpp:998
MoFEM::CoreInterface::get_field_structure
virtual const Field * get_field_structure(const std::string &name, enum MoFEMTypes bh=MF_EXIST) const =0
get field structure
MoFEM::NaturalBC
Natural boundary conditions.
Definition: Natural.hpp:57
phi
static double phi
Definition: poisson_2d_dis_galerkin.cpp:30
EshelbianPlasticity::EshelbianCore::addBoundaryFiniteElement
MoFEMErrorCode addBoundaryFiniteElement(const EntityHandle meshset=0)
Definition: EshelbianPlasticity.cpp:912
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
EshelbianPlasticity::EshelbianCore::dd_f_linear
static double dd_f_linear(const double v)
Definition: EshelbianPlasticity.hpp:937
_IT_GET_DOFS_FIELD_BY_NAME_AND_ENT_FOR_LOOP_
#define _IT_GET_DOFS_FIELD_BY_NAME_AND_ENT_FOR_LOOP_(MFIELD, NAME, ENT, IT)
loop over all dofs from a moFEM field and particular field
Definition: Interface.hpp:1917
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
EshelbianPlasticity::OpRotationBc
Definition: EshelbianPlasticity.hpp:666
EshelbianPlasticity::AddHOOps
Definition: EshelbianPlasticity.hpp:1150
EshelbianPlasticity::EshelbianCore::stretchTensor
const std::string stretchTensor
Definition: EshelbianPlasticity.hpp:977
EshelbianPlasticity::FaceUserDataOperatorStabAssembly
Definition: EshelbianPlasticity.cpp:45
EshelbianPlasticity::OpSpatialConsistencyP
Definition: EshelbianPlasticity.hpp:587
HenckyOps::d_f
auto d_f
Definition: HenckyOps.hpp:16
EshelbianPlasticity::EshelbianCore::materialVolumeElement
const std::string materialVolumeElement
Definition: EshelbianPlasticity.hpp:987
EshelbianPlasticity::EshelbianCore::setElasticElementToTs
MoFEMErrorCode setElasticElementToTs(DM dm)
Definition: EshelbianPlasticity.cpp:2192
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
EshelbianPlasticity::EshelbianCore::setElasticElementOps
MoFEMErrorCode setElasticElementOps(const int tag)
Definition: EshelbianPlasticity.cpp:2169
EshelbianMonitor.cpp
Contains definition of EshelbianMonitor class.
QUAD_::order
int order
Definition: quad.h:28
EshelbianPlasticity
Definition: CGGTonsorialBubbleBase.hpp:11
EshelbianPlasticity::EshelbianCore::d_f_linear
static double d_f_linear(const double v)
Definition: EshelbianPlasticity.hpp:936
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
EshelbianPlasticity::OpSpatialConsistencyDivTerm
Definition: EshelbianPlasticity.hpp:605
MoFEM::ForceCubitBcData
Definition of the force bc data structure.
Definition: BCData.hpp:139
EshelbianPlasticity::EshelbianCore::dd_f_log
static double dd_f_log(const double v)
Definition: EshelbianPlasticity.hpp:920
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
MoFEM::CoreInterface::add_broken_field
virtual MoFEMErrorCode add_broken_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const std::vector< std::pair< EntityType, std::function< MoFEMErrorCode(BaseFunction::DofsSideMap &)>> > list_dof_side_map, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
EshelbianPlasticity::EshelbianCore::inv_f_linear
static double inv_f_linear(const double v)
Definition: EshelbianPlasticity.hpp:939
EshelbianPlasticity::EshelbianCore::d_f
static boost::function< double(const double)> d_f
Definition: EshelbianPlasticity.hpp:904
EshelbianPlasticity::EshelbianCore::dmPrjSpatial
SmartPetscObj< DM > dmPrjSpatial
Projection spatial displacement.
Definition: EshelbianPlasticity.hpp:966
MoFEM::OpSetFlux
Definition: FormsBrokenSpaceConstraintImpl.hpp:139
NBVOLUMETET_L2
#define NBVOLUMETET_L2(P)
Number of base functions on tetrahedron for L2 space.
Definition: h1_hdiv_hcurl_l2.h:27
EshelbianPlasticity::NO_H1_CONFIGURATION
@ NO_H1_CONFIGURATION
Definition: EshelbianPlasticity.hpp:43
EshelbianPlasticity::EshelbianCore::elasticFeRhs
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeRhs
Definition: EshelbianPlasticity.hpp:957
ts_ctx
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
EshelbianPlasticity::EshelbianCore::elasticBcRhs
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcRhs
Definition: EshelbianPlasticity.hpp:960
EshelbianPlasticity::EshelbianCore::elasticFeLhs
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeLhs
Definition: EshelbianPlasticity.hpp:958
EshelbianPlasticity::EshelbianCore::frontEdges
boost::shared_ptr< Range > frontEdges
Definition: EshelbianPlasticity.hpp:1133
EshelbianPlasticity::EshelbianCore::hybridSpatialDisp
const std::string hybridSpatialDisp
Definition: EshelbianPlasticity.hpp:974
sdf.r
int r
Definition: sdf.py:8
MoFEM::CoreInterface::add_ents_to_field_by_dim
virtual MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
MoFEM::CoreInterface::add_ents_to_field_by_type
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.
FTensor::Tensor2
Definition: Tensor2_value.hpp:16
EshelbianPlasticity::CGGUserPolynomialBase::shapeFun
MatrixDouble shapeFun
Definition: EshelbianPlasticity.hpp:35
MoFEM::Field
Provide data structure for (tensor) field approximation.
Definition: FieldMultiIndices.hpp:51
USER_BASE
@ USER_BASE
user implemented approximation base
Definition: definitions.h:68
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ROW
@ ROW
Definition: definitions.h:136
EshelbianPlasticity::EshelbianCore::f_log
static double f_log(const double v)
Definition: EshelbianPlasticity.hpp:914
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::TsMonitorSet
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:259
MoFEM::OpCalculateVectorFieldValuesFromPetscVecImpl
Approximate field values for given petsc vector.
Definition: UserDataOperators.hpp:595
EshelbianPlasticity::EshelbianCore::frontAdjEdges
boost::shared_ptr< Range > frontAdjEdges
Definition: EshelbianPlasticity.hpp:1134
SideEleOp
EshelbianPlasticity::BcDisp::flags
VectorInt3 flags
Definition: EshelbianPlasticity.hpp:356
MoFEM::DMMoFEMSetDestroyProblem
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMoFEM.cpp:438
EshelbianPlasticity::OpSpatialConsistency_dBubble_dBubble
Definition: EshelbianPlasticity.hpp:792
EshelbianPlasticity::EshelbianCore::d_f_log_e
static double d_f_log_e(const double v)
Definition: EshelbianPlasticity.hpp:911
EshelbianPlasticity::BcDisp::faces
Range faces
Definition: EshelbianPlasticity.hpp:354
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
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::OpCalcNormL2Tensor2
Get norm of input MatrixDouble for Tensor2.
Definition: NormsOperators.hpp:72
EshelbianPlasticity::EshelbianCore::postProcessResults
MoFEMErrorCode postProcessResults(const int tag, const std::string file)
Definition: EshelbianPlasticity.cpp:2344
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
MoFEM::OpGetBrokenBaseSideData
Definition: FormsBrokenSpaceConstraintImpl.hpp:68
EshelbianPlasticity::OpSpatialRotation
Definition: EshelbianPlasticity.hpp:580
EshelbianPlasticity::EshelbianCore::crackFaces
boost::shared_ptr< Range > crackFaces
Definition: EshelbianPlasticity.hpp:1132
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
EshelbianPlasticity::TractionBc::faces
Range faces
Definition: EshelbianPlasticity.hpp:374
EshelbianPlasticity::OpSpatialConsistency_dBubble_domega
Definition: EshelbianPlasticity.hpp:833
MoFEM::OpBrokenLoopSide
Definition: FormsBrokenSpaceConstraintImpl.hpp:15
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::DMMoFEMTSSetIJacobian
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: DMMoFEM.cpp:853
MoFEM::DMoFEMLoopFiniteElementsUpAndLowRank
PetscErrorCode DMoFEMLoopFiniteElementsUpAndLowRank(DM dm, const char fe_name[], MoFEM::FEMethod *method, int low_rank, int up_rank, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:567
EshelbianPlasticity::EshelbianCore::SetUpSchur::createSetUpSchur
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, EshelbianCore *ep_core_ptr)
Definition: SetUpSchurImpl.cpp:515
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
TSElasticPostStep::preStepFun
static MoFEMErrorCode preStepFun(TS ts)
Definition: TSElasticPostStep.cpp:95
EshelbianPlasticity::EshelbianCore::skinElement
const std::string skinElement
Definition: EshelbianPlasticity.hpp:983
EshelbianPlasticity::EshelbianCore::d_f_log
static double d_f_log(const double v)
Definition: EshelbianPlasticity.hpp:917
EshelbianPlasticity::EshelbianCore::gradApproximator
static enum RotSelector gradApproximator
Definition: EshelbianPlasticity.hpp:897
EshelbianPlasticity::EshelbianCore::spaceOrder
int spaceOrder
Definition: EshelbianPlasticity.hpp:993
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:141
save_range
static auto save_range(moab::Interface &moab, const std::string name, const Range r)
Definition: EshelbianPlasticity.cpp:74
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
BoundaryEleOp
MoFEM::DMMoFEMTSSetI2Jacobian
PetscErrorCode DMMoFEMTSSetI2Jacobian(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: DMMoFEM.cpp:1017
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::ForcesAndSourcesCore::UserDataOperator
Definition: ForcesAndSourcesCore.hpp:549
EshelbianPlasticity::EshelbianCore::addVolumeFiniteElement
MoFEMErrorCode addVolumeFiniteElement(const EntityHandle meshset=0)
Definition: EshelbianPlasticity.cpp:850
EshelbianPlasticity::BcRot::vals
VectorDouble3 vals
Definition: EshelbianPlasticity.hpp:364
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
EshelbianPlasticity::OpSpatialConsistency_dP_domega
Definition: EshelbianPlasticity.hpp:820
MoFEM::MeshRefinement
Mesh refinement interface.
Definition: MeshRefinement.hpp:26
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
EshelbianPlasticity::VolRule::operator()
int operator()(int p_row, int p_col, int p_data) const
Definition: EshelbianPlasticity.cpp:1354
MoFEM::DMMoFEMTSSetI2Function
PetscErrorCode DMMoFEMTSSetI2Function(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 implicit function evaluation function
Definition: DMMoFEM.cpp:975
MoFEM::CoreTmp
Definition: Core.hpp:36
EshelbianPlasticity::EshelbianCore::stretchSelector
static enum StretchSelector stretchSelector
Definition: EshelbianPlasticity.hpp:898
MoFEM::get_temp_meshset_ptr
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition: Templates.hpp:1886
EshelbianPlasticity::EshelbianCore::dmElastic
SmartPetscObj< DM > dmElastic
Elastic problem.
Definition: EshelbianPlasticity.hpp:964
EshelbianPlasticity::EshelbianCore::elasticBcLhs
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcLhs
Definition: EshelbianPlasticity.hpp:959
EshelbianPlasticity::EshelbianCore::contactTreeRhs
boost::shared_ptr< ContactTree > contactTreeRhs
Make a contact tree.
Definition: EshelbianPlasticity.hpp:961
MoFEM::DMMoFEMCreateSubDM
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMoFEM.cpp:215
double
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
convert.type
type
Definition: convert.py:64
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
MoFEM::ForcesAndSourcesCore::dataOnElement
const std::array< boost::shared_ptr< EntitiesFieldData >, LASTSPACE > dataOnElement
Entity data on element entity rows fields.
Definition: ForcesAndSourcesCore.hpp:461
MoFEM::PairNameFEMethodPtr
Definition: AuxPETSc.hpp:12
EshelbianPlasticity::OpBrokenTractionBc
Definition: EshelbianPlasticity.hpp:674
EshelbianPlasticity::EshelbianCore::skeletonFaces
boost::shared_ptr< Range > skeletonFaces
Definition: EshelbianPlasticity.hpp:1136
EshelbianPlasticity::EshelbianCore::rotSelector
static enum RotSelector rotSelector
Definition: EshelbianPlasticity.hpp:896
EshelbianPlasticity::OpConstrainBoundaryHDivLhs_dU
Definition: EshelbianContact.hpp:198
EshelbianPlasticity::EshelbianCore::f
static boost::function< double(const double)> f
Definition: EshelbianPlasticity.hpp:903
EshelbianPlasticity::EshelbianCore::bcSpatialDispVecPtr
boost::shared_ptr< BcDispVec > bcSpatialDispVecPtr
Definition: EshelbianPlasticity.hpp:1005
COL
@ COL
Definition: definitions.h:136
EshelbianPlasticity::OpConstrainBoundaryHDivRhs
Definition: EshelbianContact.hpp:134
EshelbianPlasticity::VolUserDataOperatorStabAssembly
Definition: EshelbianPlasticity.cpp:41
EshelbianPlasticity::CGGUserPolynomialBase::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, BaseFunctionUnknownInterface **iface) const
Definition: EshelbianPlasticity.cpp:1372
EshelbianPlasticity::EshelbianCore::inv_d_f
static boost::function< double(const double)> inv_d_f
Definition: EshelbianPlasticity.hpp:907
EshelbianPlasticity::EshelbianCore::materialSkeletonElement
const std::string materialSkeletonElement
Definition: EshelbianPlasticity.hpp:988
SideEle
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition: plastic.cpp:61
MoFEM::OpCalculateHVecTensorField
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:2572
EshelbianPlasticity::StretchSelector
StretchSelector
Definition: EshelbianPlasticity.hpp:44
MoFEM::CoreInterface::get_comm_size
virtual int get_comm_size() const =0
EshelbianPlasticity::EshelbianCore::bcSpatialTraction
boost::shared_ptr< TractionBcVec > bcSpatialTraction
Definition: EshelbianPlasticity.hpp:1007
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
EshelbianPlasticity::EshelbianCore::mField
MoFEM::Interface & mField
Definition: EshelbianPlasticity.hpp:952
QUAD_::npoints
int npoints
Definition: quad.h:29
EshelbianPlasticity::OpDispBc
Definition: EshelbianPlasticity.hpp:635
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
EshelbianPlasticity::EshelbianCore::noStretch
static PetscBool noStretch
Definition: EshelbianPlasticity.hpp:899
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dU
Definition: EshelbianContact.hpp:155
EshelbianPlasticity::EshelbianCore::solveElastic
MoFEMErrorCode solveElastic(TS ts, Vec x)
Definition: EshelbianPlasticity.cpp:2225
EshelbianPlasticity.hpp
Eshelbian plasticity interface.
EshelbianPlasticity::BcDisp::BcDisp
BcDisp(std::string name, std::vector< double > &attr, Range &faces)
Definition: EshelbianPlasticity.cpp:1237
MoFEM::PostProcBrokenMeshInMoabBaseEnd
PostProcBrokenMeshInMoabBaseEndImpl< PostProcBrokenMeshInMoabBase< ForcesAndSourcesCore > > PostProcBrokenMeshInMoabBaseEnd
Enable to run stack of post-processing elements. Use this to end stack.
Definition: PostProcBrokenMeshInMoabBase.hpp:952
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:546
EshelbianPlasticity::CGGUserPolynomialBase::getValueHdivForCGGBubble
MoFEMErrorCode getValueHdivForCGGBubble(MatrixDouble &pts)
Definition: EshelbianPlasticity.cpp:1423
MoFEM::OpCalcNormL2Tensor1
Get norm of input MatrixDouble for Tensor1.
Definition: NormsOperators.hpp:44
EshelbianPlasticity::EshelbianCore::setSingularity
static PetscBool setSingularity
Definition: EshelbianPlasticity.hpp:900
HenckyOps::dd_f
auto dd_f
Definition: HenckyOps.hpp:17
EshelbianPlasticity::FaceRule
Definition: EshelbianPlasticity.cpp:1362
EshelbianPlasticity::EshelbianCore::contactElement
const std::string contactElement
Definition: EshelbianPlasticity.hpp:985
EshelbianPlasticity::CGG_BubbleBase_MBTET
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.
Definition: CGGTonsorialBubbleBase.cpp:20
EshelbianPlasticity::EshelbianCore::setFaceElementOps
MoFEMErrorCode setFaceElementOps(const bool add_elastic, const bool add_material, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_lhs)
Definition: EshelbianPlasticity.cpp:2005
EshelbianPlasticity::EshelbianCore::dM
SmartPetscObj< DM > dM
Coupled problem all fields.
Definition: EshelbianPlasticity.hpp:963
EshelbianPlasticity::SetIntegrationAtFrontVolume::VolFe
Definition: EshelbianPlasticity.cpp:325
Legendre_polynomials
PetscErrorCode Legendre_polynomials(int p, double s, double *diff_s, double *L, double *diffL, const int dim)
Calculate Legendre approximation basis.
Definition: base_functions.c:15
MoFEM::DMoFEMPostProcessFiniteElements
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:556
EshelbianPlasticity::OpCalculateEshelbyStress
Definition: EshelbianPlasticity.hpp:545
EshelbianPlasticity::OpSpatialPhysical_du_domega
Definition: EshelbianPlasticity.hpp:736
MoFEM::BaseFunctionUnknownInterface
Definition: BaseFunction.hpp:13
EshelbianPlasticity::EshelbianCore::rotAxis
const std::string rotAxis
Definition: EshelbianPlasticity.hpp:978
MoFEM::DMMoFEMCreateMoFEM
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: DMMoFEM.cpp:114
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
EshelbianPlasticity::OpPostProcDataStructure
Definition: EshelbianPlasticity.hpp:846
EshelbianPlasticity::EshelbianCore::spatialL2Disp
const std::string spatialL2Disp
Definition: EshelbianPlasticity.hpp:970
BiLinearForm
QUAD_3D_TABLE
static QUAD *const QUAD_3D_TABLE[]
Definition: quad.h:187
EshelbianPlasticity::EshelbianCore::contactFaces
boost::shared_ptr< Range > contactFaces
Definition: EshelbianPlasticity.hpp:1131
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
FractureMechanics::mapRefCoords
static std::map< long int, MatrixDouble > mapRefCoords
Definition: CrackFrontElement.cpp:53
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
EshelbianPlasticity::EshelbianCore::bcSpatialFreeTraction
boost::shared_ptr< TractionFreeBc > bcSpatialFreeTraction
Definition: EshelbianPlasticity.hpp:1008
MoFEM::OpCalculateHVecTensorTrace
Calculate trace of vector (Hdiv/Hcurl) space.
Definition: UserDataOperators.hpp:2760
EshelbianPlasticity::SetIntegrationAtFrontVolume::operator()
MoFEMErrorCode operator()(ForcesAndSourcesCore *fe_raw_ptr, int order_row, int order_col, int order_data)
Definition: EshelbianPlasticity.cpp:107
EshelbianPlasticity::EshelbianCore::dd_f_log_e
static double dd_f_log_e(const double v)
Definition: EshelbianPlasticity.hpp:912
TSElasticPostStep::postStepDestroy
static MoFEMErrorCode postStepDestroy()
Definition: TSElasticPostStep.cpp:85
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
EshelbianPlasticity::EshelbianCore::bubbleField
const std::string bubbleField
Definition: EshelbianPlasticity.hpp:979
MoFEM::UnknownInterface
base class for all interface classes
Definition: UnknownInterface.hpp:34
cholesky.hpp
cholesky decomposition
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into EntitiesFieldData
Definition: ForcesAndSourcesCore.hpp:22
QUAD_3D_TABLE_SIZE
#define QUAD_3D_TABLE_SIZE
Definition: quad.h:186
EshelbianPlasticity::EshelbianCore
Definition: EshelbianPlasticity.hpp:894
MoFEM::OpCalculateHVecTensorDivergence
Calculate divergence of tonsorial field using vectorial base.
Definition: UserDataOperators.hpp:2692
EshelbianPlasticity::EshelbianCore::inv_f
static boost::function< double(const double)> inv_f
Definition: EshelbianPlasticity.hpp:906
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
MoFEM::Tools::tetVolume
static double tetVolume(const double *coords)
Calculate volume of tetrahedron.
Definition: Tools.cpp:30
Range
EshelbianPlasticity::CGGUserPolynomialBase::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: EshelbianPlasticity.cpp:1380
EshelbianPlasticity::SMALL_ROT
@ SMALL_ROT
Definition: EshelbianPlasticity.hpp:43
FTensor::dd
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
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
CGGTonsorialBubbleBase.hpp
Implementation of tonsorial bubble base div(v) = 0.
EshelbianPlasticity::CGGUserPolynomialBase::cachePhiPtr
boost::shared_ptr< CachePhi > cachePhiPtr
Definition: EshelbianPlasticity.hpp:36
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
EshelbianPlasticity::EshelbianCore::physicalEquations
boost::shared_ptr< PhysicalEquations > physicalEquations
Definition: EshelbianPlasticity.hpp:955
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
EshelbianPlasticity::OpConstrainBoundaryL2Rhs
Definition: EshelbianContact.hpp:118
EshelbianPlasticity::LARGE_ROT
@ LARGE_ROT
Definition: EshelbianPlasticity.hpp:43
EshelbianPlasticity::LOG
@ LOG
Definition: EshelbianPlasticity.hpp:44
EshelbianPlasticity::EshelbianCore::elementVolumeName
const std::string elementVolumeName
Definition: EshelbianPlasticity.hpp:981
EshelbianPlasticity::EshelbianCore::frontVertices
boost::shared_ptr< Range > frontVertices
Definition: EshelbianPlasticity.hpp:1135
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
EshelbianPlasticity::EshelbianCore::frontLayers
int frontLayers
Definition: EshelbianPlasticity.hpp:1001
NBVOLUMETET_CCG_BUBBLE
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
Definition: CGGTonsorialBubbleBase.hpp:19
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
EshelbianPlasticity::TractionBc::TractionBc
TractionBc(std::string name, std::vector< double > &attr, Range &faces)
Definition: EshelbianPlasticity.cpp:1263
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
eps
static const double eps
Definition: check_base_functions_derivatives_on_tet.cpp:11
EshelbianPlasticity::RotSelector
RotSelector
Definition: EshelbianPlasticity.hpp:43
EshelbianPlasticity::EshelbianCore::bitAdjParent
BitRefLevel bitAdjParent
bit ref level for parent
Definition: EshelbianPlasticity.hpp:1142
EshelbianPlasticity::TractionBc::flags
VectorInt3 flags
Definition: EshelbianPlasticity.hpp:376
DISCONTINUOUS
@ DISCONTINUOUS
Broken continuity (No effect on L2 space)
Definition: definitions.h:101
EshelbianPlasticity::EshelbianCore::bitAdjEntMask
BitRefLevel bitAdjEntMask
bit ref level for parent parent
Definition: EshelbianPlasticity.hpp:1146
EshelbianPlasticity::EshelbianCore::inv_d_f_linear
static double inv_d_f_linear(const double v)
Definition: EshelbianPlasticity.hpp:940
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
ShapeDiffMBTET
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:319
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
EshelbianContact.cpp
MoFEM::DMMoFEMTSSetIFunction
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: DMMoFEM.cpp:800
MoFEM::OpCalculateTensor2SymmetricFieldValues
Calculate symmetric tensor field values at integration pts.
Definition: UserDataOperators.hpp:977
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
EshelbianContact.hpp
MoFEM::ForcesAndSourcesCore::gaussPts
MatrixDouble gaussPts
Matrix of integration points.
Definition: ForcesAndSourcesCore.hpp:109
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
EshelbianPlasticity::EshelbianCore::alphaU
double alphaU
Definition: EshelbianPlasticity.hpp:996
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
EshelbianMonitor
Definition: EshelbianMonitor.cpp:6
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
EshelbianPlasticity::EshelbianCore::skeletonElement
const std::string skeletonElement
Definition: EshelbianPlasticity.hpp:984
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
sdf_hertz_2d_axisymm_plane.d
float d
Definition: sdf_hertz_2d_axisymm_plane.py:4
EshelbianPlasticity::EshelbianCore::inv_d_f_log
static double inv_d_f_log(const double v)
Definition: EshelbianPlasticity.hpp:928
MoFEM::DMMoFEMGetTsCtx
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMoFEM.cpp:1142
EshelbianPlasticity::EshelbianCore::addFields
MoFEMErrorCode addFields(const EntityHandle meshset=0)
Definition: EshelbianPlasticity.cpp:483
TSElasticPostStep::postStepInitialise
static MoFEMErrorCode postStepInitialise(EshelbianCore *ep_ptr)
Definition: TSElasticPostStep.cpp:11
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
BdyEleOp
BoundaryEle::UserDataOperator BdyEleOp
Definition: test_broken_space.cpp:37
MoFEM::EntitiesFieldData::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition: EntitiesFieldData.hpp:57
EshelbianPlasticity::EshelbianCore::inv_dd_f
static boost::function< double(const double)> inv_dd_f
Definition: EshelbianPlasticity.hpp:908
EshelbianPlasticity::EshelbianCore::inv_dd_f_log
static double inv_dd_f_log(const double v)
Definition: EshelbianPlasticity.hpp:931
EshelbianPlasticity::FaceRule::operator()
int operator()(int p_row, int p_col, int p_data) const
Definition: EshelbianPlasticity.cpp:1363
MoFEM::BLOCK_SCHUR
@ BLOCK_SCHUR
Definition: FormsIntegrators.hpp:108
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
EshelbianPlasticity::EshelbianCore::getOptions
MoFEMErrorCode getOptions()
Definition: EshelbianPlasticity.cpp:377
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
MoFEM::FEMethod::getFEEntityHandle
EntityHandle getFEEntityHandle() const
Definition: LoopMethods.hpp:460
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
EshelbianPlasticity::BcRot::theta
double theta
Definition: EshelbianPlasticity.hpp:365
QUIET
@ QUIET
Definition: definitions.h:221
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
EshelbianPlasticity::EshelbianCore::contactRefinementLevels
int contactRefinementLevels
Definition: EshelbianPlasticity.hpp:1000
EshelbianPlasticity::EshelbianCore::getBc
MoFEMErrorCode getBc(boost::shared_ptr< BC > &bc_vec_ptr, const std::string block_name, const int nb_attributes)
Definition: EshelbianPlasticity.hpp:1011
MoFEM::CoreInterface::check_finite_element
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:238
EshelbianPlasticity::EshelbianCore::parentAdjSkeletonFunctionDim2
boost::shared_ptr< ParentFiniteElementAdjacencyFunctionSkeleton< 2 > > parentAdjSkeletonFunctionDim2
Definition: EshelbianPlasticity.hpp:1140
EshelbianPlasticity::OpSpatialEquilibrium
Definition: EshelbianPlasticity.hpp:569
SetUpSchurImpl.cpp
EshelbianPlasticity::OpSpatialEquilibrium_dw_dP
Definition: EshelbianPlasticity.hpp:688
MoFEM::CoreInterface::set_field_order
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.
MoFEM::MeshsetsManager::getCubitMeshsetPtr
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
Definition: MeshsetsManager.cpp:578
EshelbianPlasticity::SetIntegrationAtFrontVolume::frontNodes
static boost::shared_ptr< Range > frontNodes
Definition: EshelbianPlasticity.cpp:334
MoFEM::EntitiesFieldData
data structure for finite element entity
Definition: EntitiesFieldData.hpp:40
EshelbianPlasticity::BcDisp::vals
VectorDouble3 vals
Definition: EshelbianPlasticity.hpp:355
EshelbianPlasticity::EshelbianCore::getSpatialTractionBc
MoFEMErrorCode getSpatialTractionBc()
Definition: EshelbianPlasticity.cpp:2599
EshelbianPlasticity::EshelbianCore::setVolumeElementOps
MoFEMErrorCode setVolumeElementOps(const int tag, const bool add_elastic, const bool add_material, boost::shared_ptr< VolumeElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > &fe_lhs)
Definition: EshelbianPlasticity.cpp:1571
EshelbianPlasticity::EshelbianCore::materialOrder
int materialOrder
Definition: EshelbianPlasticity.hpp:995
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dP
Definition: EshelbianContact.hpp:172
EshelbianPlasticity::EshelbianCore::crackElement
const std::string crackElement
Definition: EshelbianPlasticity.hpp:986
EshelbianPlasticity::OpSpatialRotation_domega_domega
Definition: EshelbianPlasticity.hpp:769
get_range_from_block
static auto get_range_from_block(MoFEM::Interface &m_field, const std::string block_name, int dim)
Definition: EshelbianPlasticity.cpp:52
EshelbianPlasticity::MODERATE_ROT
@ MODERATE_ROT
Definition: EshelbianPlasticity.hpp:43
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
EshelbianPlasticity::EshelbianCore::materialL2Disp
const std::string materialL2Disp
Definition: EshelbianPlasticity.hpp:971
EshelbianPlasticity::OpSpatialEquilibrium_dw_dw
Definition: EshelbianPlasticity.hpp:699
convert.int
int
Definition: convert.py:64
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
MoFEM::ForcesAndSourcesCore::UserDataOperator::getKSPB
Mat getKSPB() const
Definition: ForcesAndSourcesCore.hpp:1104
EshelbianPlasticity::OpSpatialRotation_domega_dBubble
Definition: EshelbianPlasticity.hpp:758
EshelbianPlasticity::EshelbianCore::setContactElementRhsOps
MoFEMErrorCode setContactElementRhsOps(boost::shared_ptr< ContactTree > &fe_contact_tree)
Definition: EshelbianPlasticity.cpp:2074
EshelbianPlasticity::EshelbianCore::eshelbyStress
const std::string eshelbyStress
Definition: EshelbianPlasticity.hpp:969
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:496
EshelbianPlasticity::SetIntegrationAtFrontVolume::frontEdges
static boost::shared_ptr< Range > frontEdges
Definition: EshelbianPlasticity.cpp:335
MoFEM::SCHUR
@ SCHUR
Definition: FormsIntegrators.hpp:106
MoFEM::CoreInterface::add_field
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.
EshelbianPlasticity::EshelbianCore::bitAdjEnt
BitRefLevel bitAdjEnt
bit ref level for parent
Definition: EshelbianPlasticity.hpp:1145
EshelbianPlasticity::OpSpatialRotation_domega_dP
Definition: EshelbianPlasticity.hpp:747
EshelbianPlasticity::EshelbianCore::contactDisp
const std::string contactDisp
Definition: EshelbianPlasticity.hpp:976
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1290
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
l
FTensor::Index< 'l', 3 > l
Definition: matrix_function.cpp:21
EshelbianPlasticity::EshelbianCore::getSpatialDispBc
MoFEMErrorCode getSpatialDispBc()
[Getting norms]
Definition: EshelbianPlasticity.cpp:2561
EshelbianPlasticity::OpSpatialConsistencyBubble
Definition: EshelbianPlasticity.hpp:596
MoFEM::PostProcBrokenMeshInMoabBaseBegin
PostProcBrokenMeshInMoabBaseBeginImpl< PostProcBrokenMeshInMoabBase< ForcesAndSourcesCore > > PostProcBrokenMeshInMoabBaseBegin
Enable to run stack of post-processing elements. Use this to begin stack.
Definition: PostProcBrokenMeshInMoabBase.hpp:934
MoFEM::ForcesAndSourcesCore::mField
Interface & mField
Definition: ForcesAndSourcesCore.hpp:24
EshelbianPlasticity::EshelbianCore::naturalBcElement
const std::string naturalBcElement
Definition: EshelbianPlasticity.hpp:982
EshelbianPlasticity::EshelbianCore::gettingNorms
MoFEMErrorCode gettingNorms()
[Getting norms]
Definition: EshelbianPlasticity.cpp:2497
EshelbianPlasticity::EshelbianCore::~EshelbianCore
virtual ~EshelbianCore()
EshelbianPlasticity::CGGUserPolynomialBase::CGGUserPolynomialBase
CGGUserPolynomialBase(boost::shared_ptr< CachePhi > cache_phi=nullptr)
Definition: EshelbianPlasticity.cpp:1368
EshelbianPlasticity::TractionBc::vals
VectorDouble3 vals
Definition: EshelbianPlasticity.hpp:375
MoFEM::TetPolynomialBase
Calculate base functions on tetrahedral.
Definition: TetPolynomialBase.hpp:17
MoFEM::ForcesAndSourcesCore::UserDataOperator::UserDataOperator
UserDataOperator(const FieldSpace space, const char type=OPSPACE, const bool symm=true)
Definition: ForcesAndSourcesCore.cpp:1977
EshelbianPlasticity::OpSpatialPhysical_du_dBubble
Definition: EshelbianPlasticity.hpp:724
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1123
EshelbianPlasticity::OpTreeSearch
Definition: EshelbianContact.hpp:237
EshelbianPlasticity::EshelbianCore::materialH1Positions
const std::string materialH1Positions
Definition: EshelbianPlasticity.hpp:973
ShapeMBTET
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:306
quad.h
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
EshelbianPlasticity::OpSpatialConsistency_dP_dP
Definition: EshelbianPlasticity.hpp:779
Generated by Doxygen 1.8.17 and hosted at University of Glasgow