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 #include <MoFEM.hpp>
9 using namespace MoFEM;
10 
11 #include <CGGTonsorialBubbleBase.hpp>
12 
13 #include <EshelbianPlasticity.hpp>
14 #include <boost/math/constants/constants.hpp>
15 
16 #include <cholesky.hpp>
17 #ifdef PYTHON_SDF
18 #include <boost/python.hpp>
19 #include <boost/python/def.hpp>
20 #include <boost/python/numpy.hpp>
21 namespace bp = boost::python;
22 namespace np = boost::python::numpy;
23 
24 #pragma message "With PYTHON_SDF"
25 
26 #else
27 
28 #pragma message "Without PYTHON_SDF"
29 
30 #endif
31 
32 #include <EshelbianContact.hpp>
33 
34 namespace EshelbianPlasticity {
35 struct VolUserDataOperatorStabAssembly
36  : public VolumeElementForcesAndSourcesCore::UserDataOperator {
37  using VolumeElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
38 };
39 struct FaceUserDataOperatorStabAssembly
40  : public FaceElementForcesAndSourcesCore::UserDataOperator {
41  using FaceElementForcesAndSourcesCore::UserDataOperator::UserDataOperator;
42 };
43 
44 } // namespace EshelbianPlasticity
45 
46 static auto get_range_from_block(MoFEM::Interface &m_field,
47  const std::string block_name) {
48  Range r;
49 
50  auto mesh_mng = m_field.getInterface<MeshsetsManager>();
51  auto bcs = mesh_mng->getCubitMeshsetPtr(
52 
53  std::regex((boost::format("%s(.*)") % block_name).str())
54 
55  );
56 
57  for (auto bc : bcs) {
58  Range faces;
59  CHK_MOAB_THROW(
60  bc->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 2, faces, true),
61  "get meshset ents");
62  r.merge(faces);
63  }
64 
65  return r;
66 };
67 
68 template <>
69 typename MoFEM::OpBaseImpl<
70  SCHUR,
71  EshelbianPlasticity::VolUserDataOperatorStabAssembly>::MatSetValuesHook
72  MoFEM::OpBaseImpl<SCHUR,
73  EshelbianPlasticity::VolUserDataOperatorStabAssembly>::
74  matSetValuesHook = [](ForcesAndSourcesCore::UserDataOperator *op_ptr,
75  const EntitiesFieldData::EntData &row_data,
76  const EntitiesFieldData::EntData &col_data,
77  MatrixDouble &m) {
78  return MatSetValues<AssemblyTypeSelector<SCHUR>>(
79  op_ptr->getKSPB(), row_data, col_data, m, ADD_VALUES);
80  };
81 
82 namespace EshelbianPlasticity {
83 
84 enum RotSelector EshelbianCore::rotSelector = LARGE_ROT;
85 enum RotSelector EshelbianCore::gradApperoximator = LARGE_ROT;
86 enum StretchSelector EshelbianCore::stretchSelector = LOG;
87 
88 double EshelbianCore::exponentBase = exp(1);
89 boost::function<double(const double)> EshelbianCore::f = EshelbianCore::f_log;
90 boost::function<double(const double)> EshelbianCore::d_f =
91  EshelbianCore::d_f_log;
92 boost::function<double(const double)> EshelbianCore::dd_f =
93  EshelbianCore::dd_f_log;
94 boost::function<double(const double)> EshelbianCore::inv_f =
95  EshelbianCore::inv_f_log;
96 boost::function<double(const double)> EshelbianCore::inv_d_f =
97  EshelbianCore::inv_d_f_log;
98 boost::function<double(const double)> EshelbianCore::inv_dd_f =
99  EshelbianCore::inv_dd_f_log;
100 
101 MoFEMErrorCode
102 EshelbianCore::query_interface(boost::typeindex::type_index type_index,
103  UnknownInterface **iface) const {
104  *iface = const_cast<EshelbianCore *>(this);
105  return 0;
106 }
107 
108 MoFEMErrorCode OpJacobian::doWork(int side, EntityType type, EntData &data) {
109  MoFEMFunctionBegin;
110 
111  if (evalRhs)
112  CHKERR evaluateRhs(data);
113 
114  if (evalLhs)
115  CHKERR evaluateLhs(data);
116 
117  MoFEMFunctionReturn(0);
118 }
119 
120 EshelbianCore::EshelbianCore(MoFEM::Interface &m_field)
121  : mField(m_field), piolaStress("P"), eshelbyStress("S"),
122  spatialL2Disp("wL2"), spatialH1Disp("wH1"), contactDisp("contactDisp"),
123  materialL2Disp("W"), stretchTensor("u"), rotAxis("omega"),
124  materialGradient("G"), tauField("TAU"), lambdaField("LAMBDA"),
125  bubbleField("BUBBLE"), elementVolumeName("EP"),
126  naturalBcElement("NATURAL_BC"), essentialBcElement("ESSENTIAL_BC"),
127  skinElement("SKIN_ELEMENT"), contactElement("CONTACT_ELEMENT") {
128 
129  ierr = getOptions();
130  CHKERRABORT(PETSC_COMM_WORLD, ierr);
131 }
132 
133 EshelbianCore::~EshelbianCore() = default;
134 
135 MoFEMErrorCode EshelbianCore::getOptions() {
136  MoFEMFunctionBegin;
137  const char *list_rots[] = {"small", "moderate", "large"};
138  PetscInt choice_rot = EshelbianCore::rotSelector;
139  PetscInt choice_grad = EshelbianCore::gradApperoximator;
140 
141  const char *list_stretches[] = {"linear", "log"};
142  PetscInt choice_stretch = StretchSelector::LOG;
143 
144  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
145  "none");
146 
147  spaceOrder = 2;
148  CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
149  spaceOrder, &spaceOrder, PETSC_NULL);
150  spaceH1Order = -1;
151  CHKERR PetscOptionsInt("-space_h1_order", "approximation oder for space", "",
152  spaceH1Order, &spaceH1Order, PETSC_NULL);
153 
154  materialOrder = 1;
155  CHKERR PetscOptionsInt("-material_order", "approximation oder for material",
156  "", materialOrder, &materialOrder, PETSC_NULL);
157 
158  alphaU = 0;
159  CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alphaU,
160  &alphaU, PETSC_NULL);
161 
162  alphaW = 0;
163  CHKERR PetscOptionsScalar("-viscosity_alpha_w", "viscosity", "", alphaW,
164  &alphaW, PETSC_NULL);
165 
166  alphaRho = 0;
167  CHKERR PetscOptionsScalar("-density_alpha_rho", "density", "", alphaRho,
168  &alphaRho, PETSC_NULL);
169 
170  precEps = 0;
171  CHKERR PetscOptionsScalar("-preconditioner_eps", "preconditioner_eps", "",
172  precEps, &precEps, PETSC_NULL);
173 
174  precEpsOmega = 0;
175  CHKERR PetscOptionsScalar("-preconditioner_eps_omega", "preconditioner_eps",
176  "", precEpsOmega, &precEpsOmega, PETSC_NULL);
177  precEpsW = 0;
178  CHKERR PetscOptionsScalar("-preconditioner_eps_w", "preconditioner_eps", "",
179  precEpsW, &precEpsW, PETSC_NULL);
180  precEpsContactDisp = 0;
181  CHKERR PetscOptionsScalar("-preconditioner_eps_contact_disp",
182  "preconditioner_eps", "", precEpsContactDisp,
183  &precEpsContactDisp, PETSC_NULL);
184 
185  CHKERR PetscOptionsEList("-rotations", "rotations", "", list_rots,
186  LARGE_ROT + 1, list_rots[choice_rot], &choice_rot,
187  PETSC_NULL);
188  CHKERR PetscOptionsEList("-grad", "gradient of defamation approximate", "",
189  list_rots, LARGE_ROT + 1, list_rots[choice_grad],
190  &choice_grad, PETSC_NULL);
191 
192  CHKERR PetscOptionsScalar("-exponent_base", "exponent_base", "", exponentBase,
193  &EshelbianCore::exponentBase, PETSC_NULL);
194  CHKERR PetscOptionsEList(
195  "-stretches", "stretches", "", list_stretches, StretchSelector::LOG + 1,
196  list_stretches[choice_stretch], &choice_stretch, PETSC_NULL);
197 
198  ierr = PetscOptionsEnd();
199  CHKERRG(ierr);
200 
201  EshelbianCore::rotSelector = static_cast<RotSelector>(choice_rot);
202  EshelbianCore::gradApperoximator = static_cast<RotSelector>(choice_grad);
203  EshelbianCore::stretchSelector = static_cast<StretchSelector>(choice_stretch);
204 
205  switch (EshelbianCore::stretchSelector) {
206  case StretchSelector::LINEAR:
207  EshelbianCore::f = f_linear;
208  EshelbianCore::d_f = d_f_linear;
209  EshelbianCore::dd_f = dd_f_linear;
210  EshelbianCore::inv_f = inv_f_linear;
211  EshelbianCore::inv_d_f = inv_d_f_linear;
212  EshelbianCore::inv_dd_f = inv_dd_f_linear;
213  break;
214  case StretchSelector::LOG:
215  EshelbianCore::f = f_log;
216  EshelbianCore::d_f = d_f_log;
217  EshelbianCore::dd_f = dd_f_log;
218  EshelbianCore::inv_f = inv_f_log;
219  EshelbianCore::inv_d_f = inv_d_f_log;
220  EshelbianCore::inv_dd_f = inv_dd_f_log;
221  break;
222  default:
223  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "Unknown stretch");
224  break;
225  };
226 
227  precEpsOmega += precEps;
228  precEpsW += precEps;
229 
230  MOFEM_LOG("EP", Sev::inform) << "spaceOrder " << spaceOrder;
231  MOFEM_LOG("EP", Sev::inform) << "spaceH1Order " << spaceH1Order;
232  MOFEM_LOG("EP", Sev::inform) << "materialOrder " << materialOrder;
233  MOFEM_LOG("EP", Sev::inform) << "alphaU " << alphaU;
234  MOFEM_LOG("EP", Sev::inform) << "alphaW " << alphaW;
235  MOFEM_LOG("EP", Sev::inform) << "alphaRho " << alphaRho;
236  MOFEM_LOG("EP", Sev::inform) << "precEps " << precEps;
237  MOFEM_LOG("EP", Sev::inform) << "precEpsOmega " << precEpsOmega;
238  MOFEM_LOG("EP", Sev::inform) << "precEpsW " << precEpsW;
239  MOFEM_LOG("EP", Sev::inform)
240  << "Rotations " << list_rots[EshelbianCore::rotSelector];
241  MOFEM_LOG("EP", Sev::inform) << "Gradient of deformation "
242  << list_rots[EshelbianCore::gradApperoximator];
243  if (exponentBase != exp(1))
244  MOFEM_LOG("EP", Sev::inform)
245  << "Base exponent " << EshelbianCore::exponentBase;
246  else
247  MOFEM_LOG("EP", Sev::inform) << "Base exponent e";
248  MOFEM_LOG("EP", Sev::inform) << "Stretch " << list_stretches[choice_stretch];
249 
250  if (spaceH1Order == -1)
251  spaceH1Order = spaceOrder;
252 
253  MoFEMFunctionReturn(0);
254 }
255 
256 MoFEMErrorCode EshelbianCore::addFields(const EntityHandle meshset) {
257  MoFEMFunctionBegin;
258 
259  Range tets;
260  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
261  Range tets_skin_part;
262  Skinner skin(&mField.get_moab());
263  CHKERR skin.find_skin(0, tets, false, tets_skin_part);
264  ParallelComm *pcomm =
265  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
266  Range tets_skin;
267  CHKERR pcomm->filter_pstatus(tets_skin_part,
268  PSTATUS_SHARED | PSTATUS_MULTISHARED,
269  PSTATUS_NOT, -1, &tets_skin);
270 
271  if (bcSpatialDispVecPtr)
272  for (auto &v : *bcSpatialDispVecPtr) {
273  tets_skin = subtract(tets_skin, v.faces);
274  }
275  if (bcSpatialRotationVecPtr)
276  for (auto &v : *bcSpatialRotationVecPtr) {
277  tets_skin = subtract(tets_skin, v.faces);
278  }
279  if (bcSpatialTraction)
280  for (auto &v : *bcSpatialTraction) {
281  tets_skin = subtract(tets_skin, v.faces);
282  }
283 
284  auto subtract_faces_where_displacements_are_applied =
285  [&](const std::string block_name) {
286  MoFEMFunctionBegin;
287  auto contact_range = get_range_from_block(mField, block_name);
288  tets_skin = subtract(tets_skin, contact_range);
289  MoFEMFunctionReturn(0);
290  };
291  CHKERR subtract_faces_where_displacements_are_applied("CONTACT");
292 
293  Range faces;
294  CHKERR mField.get_moab().get_adjacencies(tets, 2, true, faces,
295  moab::Interface::UNION);
296  Range faces_not_on_the_skin = subtract(faces, tets_skin);
297 
298  auto add_hdiv_field = [&](const std::string field_name, const int order,
299  const int dim) {
300  MoFEMFunctionBegin;
301  CHKERR mField.add_field(field_name, HDIV, DEMKOWICZ_JACOBI_BASE, dim,
302  MB_TAG_SPARSE, MF_ZERO);
303  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
304  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
305  CHKERR mField.set_field_order(faces_not_on_the_skin, field_name, order);
306  CHKERR mField.set_field_order(tets_skin, field_name, 0);
307  MoFEMFunctionReturn(0);
308  };
309 
310  auto add_l2_field = [this, meshset](const std::string field_name,
311  const int order, const int dim) {
312  MoFEMFunctionBegin;
313  CHKERR mField.add_field(field_name, L2, AINSWORTH_LEGENDRE_BASE, dim,
314  MB_TAG_DENSE, MF_ZERO);
315  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
316  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
317  MoFEMFunctionReturn(0);
318  };
319 
320  auto add_h1_field = [this, meshset](const std::string field_name,
321  const int order, const int dim) {
322  MoFEMFunctionBegin;
323  CHKERR mField.add_field(field_name, H1, AINSWORTH_LEGENDRE_BASE, dim,
324  MB_TAG_DENSE, MF_ZERO);
325  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
326  CHKERR mField.set_field_order(meshset, MBVERTEX, field_name, 1);
327  CHKERR mField.set_field_order(meshset, MBEDGE, field_name, order);
328  CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
329  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
330  MoFEMFunctionReturn(0);
331  };
332 
333  auto add_l2_field_by_range = [this](const std::string field_name,
334  const int order, const int dim,
335  const int field_dim, Range &&r) {
336  MoFEMFunctionBegin;
337  CHKERR mField.add_field(field_name, L2, AINSWORTH_LEGENDRE_BASE, field_dim,
338  MB_TAG_SPARSE, MF_ZERO);
339  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(r);
340  CHKERR mField.add_ents_to_field_by_dim(r, dim, field_name);
341  CHKERR mField.set_field_order(r, field_name, order);
342  MoFEMFunctionReturn(0);
343  };
344 
345  auto add_bubble_field = [this, meshset](const std::string field_name,
346  const int order, const int dim) {
347  MoFEMFunctionBegin;
348  CHKERR mField.add_field(field_name, HDIV, USER_BASE, dim, MB_TAG_DENSE,
349  MF_ZERO);
350  // Modify field
351  auto field_ptr = mField.get_field_structure(field_name);
352  auto field_order_table =
353  const_cast<Field *>(field_ptr)->getFieldOrderTable();
354  auto get_cgg_bubble_order_zero = [](int p) { return 0; };
355  auto get_cgg_bubble_order_tet = [](int p) {
356  return NBVOLUMETET_CCG_BUBBLE(p);
357  };
358  field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
359  field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
360  field_order_table[MBTRI] = get_cgg_bubble_order_zero;
361  field_order_table[MBTET] = get_cgg_bubble_order_tet;
362  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
363  CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
364  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
365  MoFEMFunctionReturn(0);
366  };
367 
368  // spatial fields
369  CHKERR add_hdiv_field(piolaStress, spaceOrder, 3);
370  CHKERR add_bubble_field(bubbleField, spaceOrder, 1);
371  CHKERR add_l2_field(spatialL2Disp, spaceOrder - 1, 3);
372  CHKERR add_l2_field(rotAxis, spaceOrder - 1, 3);
373  CHKERR add_l2_field(stretchTensor, spaceOrder, 6);
374  CHKERR add_l2_field_by_range(contactDisp, spaceOrder - 1, 2, 3,
375  get_range_from_block(mField, "CONTACT"));
376 
377  // material fields
378  // CHKERR add_hdiv_field(eshelbyStress, materialOrder, 3);
379  // CHKERR add_l2_field(materialGradient, materialOrder - 1, 9);
380  // CHKERR add_l2_field(materialL2Disp, materialOrder - 1, 3);
381  // CHKERR add_l2_field(tauField, materialOrder - 1, 1);
382  // CHKERR add_l2_field(lambdaField, materialOrder - 1, 1);
383 
384  // Add history filedes
385  // CHKERR add_l2_field(materialGradient + "0", materialOrder - 1, 9);
386 
387  // spatial displacement
388  CHKERR add_h1_field(spatialH1Disp, spaceH1Order, 3);
389 
390  CHKERR mField.build_fields();
391 
392  MoFEMFunctionReturn(0);
393 }
394 
395 MoFEMErrorCode
396 EshelbianCore::addVolumeFiniteElement(const EntityHandle meshset) {
397  MoFEMFunctionBegin;
398 
399  // set finite element fields
400  auto add_field_to_fe = [this](const std::string fe,
401  const std::string field_name) {
402  MoFEMFunctionBegin;
403  CHKERR mField.modify_finite_element_add_field_row(fe, field_name);
404  CHKERR mField.modify_finite_element_add_field_col(fe, field_name);
405  CHKERR mField.modify_finite_element_add_field_data(fe, field_name);
406  MoFEMFunctionReturn(0);
407  };
408 
409  if (!mField.check_finite_element(elementVolumeName)) {
410  CHKERR mField.add_finite_element(elementVolumeName, MF_ZERO);
411  CHKERR mField.add_ents_to_finite_element_by_type(meshset, MBTET,
412  elementVolumeName);
413 
414  CHKERR add_field_to_fe(elementVolumeName, piolaStress);
415  CHKERR add_field_to_fe(elementVolumeName, bubbleField);
416  // CHKERR add_field_to_fe(elementVolumeName, eshelbyStress);
417  CHKERR add_field_to_fe(elementVolumeName, stretchTensor);
418  CHKERR add_field_to_fe(elementVolumeName, rotAxis);
419  CHKERR add_field_to_fe(elementVolumeName, spatialL2Disp);
420  CHKERR add_field_to_fe(elementVolumeName, stretchTensor);
421  CHKERR add_field_to_fe(elementVolumeName, spatialH1Disp);
422  CHKERR add_field_to_fe(elementVolumeName, contactDisp);
423  // CHKERR add_field_to_fe(elementVolumeName, materialGradient);
424  // CHKERR mField.modify_finite_element_add_field_data(elementVolumeName,
425  // materialGradient +
426  // "0");
427  }
428 
429  // build finite elements data structures
430  CHKERR mField.build_finite_elements(elementVolumeName);
431 
432  MoFEMFunctionReturn(0);
433 }
434 
435 MoFEMErrorCode
436 EshelbianCore::addBoundaryFiniteElement(const EntityHandle meshset) {
437  MoFEMFunctionBegin;
438 
439  // set finite element fields
440  auto add_field_to_fe = [this](const std::string fe,
441  const std::string field_name) {
442  MoFEMFunctionBegin;
443  CHKERR mField.modify_finite_element_add_field_row(fe, field_name);
444  CHKERR mField.modify_finite_element_add_field_col(fe, field_name);
445  CHKERR mField.modify_finite_element_add_field_data(fe, field_name);
446  MoFEMFunctionReturn(0);
447  };
448 
449  Range natural_bc_elements;
450  if (bcSpatialDispVecPtr) {
451  for (auto &v : *bcSpatialDispVecPtr) {
452  natural_bc_elements.merge(v.faces);
453  }
454  }
455  if (bcSpatialRotationVecPtr) {
456  for (auto &v : *bcSpatialRotationVecPtr) {
457  natural_bc_elements.merge(v.faces);
458  }
459  }
460  Range essentail_bc_elements;
461  if (bcSpatialTraction) {
462  for (auto &v : *bcSpatialTraction) {
463  essentail_bc_elements.merge(v.faces);
464  }
465  }
466 
467  CHKERR mField.add_finite_element(naturalBcElement, MF_ZERO);
468  CHKERR mField.add_ents_to_finite_element_by_type(natural_bc_elements, MBTRI,
469  naturalBcElement);
470  CHKERR add_field_to_fe(naturalBcElement, piolaStress);
471  // CHKERR add_field_to_fe(naturalBcElement, eshelbyStress);
472  CHKERR mField.build_finite_elements(naturalBcElement);
473 
474  CHKERR mField.add_finite_element(essentialBcElement, MF_ZERO);
475  CHKERR mField.add_ents_to_finite_element_by_type(essentail_bc_elements, MBTRI,
476  essentialBcElement);
477  CHKERR add_field_to_fe(essentialBcElement, piolaStress);
478  // CHKERR add_field_to_fe(essentialBcElement, eshelbyStress);
479  CHKERR mField.build_finite_elements(essentialBcElement);
480 
481  auto get_skin = [&]() {
482  Range body_ents;
483  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
484  Skinner skin(&mField.get_moab());
485  Range skin_ents;
486  CHKERR skin.find_skin(0, body_ents, false, skin_ents);
487  return skin_ents;
488  };
489 
490  auto filter_true_skin = [&](auto &&skin) {
491  Range boundary_ents;
492  ParallelComm *pcomm =
493  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
494  CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
495  PSTATUS_NOT, -1, &boundary_ents);
496  return boundary_ents;
497  };
498 
499  auto skin = filter_true_skin(get_skin());
500 
501  CHKERR mField.add_finite_element(skinElement, MF_ZERO);
502  CHKERR mField.add_ents_to_finite_element_by_type(skin, MBTRI, skinElement);
503  CHKERR add_field_to_fe(skinElement, piolaStress);
504  CHKERR add_field_to_fe(skinElement, spatialH1Disp);
505  CHKERR add_field_to_fe(skinElement, contactDisp);
506  // CHKERR add_field_to_fe(skinElement, eshelbyStress);
507  CHKERR mField.build_finite_elements(skinElement);
508 
509  auto contact_range = get_range_from_block(mField, "CONTACT");
510  MOFEM_LOG("EP", Sev::inform) << "Contact elements " << contact_range.size();
511 
512  CHKERR mField.add_finite_element(contactElement, MF_ZERO);
513  CHKERR mField.add_ents_to_finite_element_by_type(contact_range, MBTRI,
514  contactElement);
515  CHKERR add_field_to_fe(contactElement, piolaStress);
516  CHKERR add_field_to_fe(contactElement, spatialH1Disp);
517  CHKERR add_field_to_fe(contactElement, contactDisp);
518  CHKERR mField.build_finite_elements(contactElement);
519 
520  MoFEMFunctionReturn(0);
521 }
522 
523 MoFEMErrorCode EshelbianCore::addDMs(const BitRefLevel bit) {
524  MoFEMFunctionBegin;
525 
526  // find adjacencies between finite elements and dofs
527  CHKERR mField.build_adjacencies(bit, QUIET);
528 
529  // Create coupled problem
530  dM = createDM(mField.get_comm(), "DMMOFEM");
531  CHKERR DMMoFEMCreateMoFEM(dM, &mField, "ESHELBY_PLASTICITY", bit,
532  BitRefLevel().set());
533  CHKERR DMMoFEMSetDestroyProblem(dM, PETSC_TRUE);
534  CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
535  CHKERR DMMoFEMAddElement(dM, elementVolumeName);
536  CHKERR DMMoFEMAddElement(dM, naturalBcElement);
537  CHKERR DMMoFEMAddElement(dM, essentialBcElement);
538  CHKERR DMMoFEMAddElement(dM, contactElement);
539  CHKERR DMMoFEMAddElement(dM, skinElement);
540  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
541  CHKERR DMSetUp(dM);
542  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
543 
544  auto remove_dofs_on_essential_spatial_stress_boundary =
545  [&](const std::string prb_name) {
546  MoFEMFunctionBegin;
547  for (int d : {0, 1, 2})
548  CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
549  prb_name, piolaStress, (*bcSpatialFreeTraction)[d], d, d, 0,
550  MAX_DOFS_ON_ENTITY, NOISY, true);
551  MoFEMFunctionReturn(0);
552  };
553  CHKERR remove_dofs_on_essential_spatial_stress_boundary("ESHELBY_PLASTICITY");
554 
555  auto contact_range = get_range_from_block(mField, "CONTACT");
556 
557  // Create elastic sub-problem
558  dmElastic = createDM(mField.get_comm(), "DMMOFEM");
559  CHKERR DMMoFEMCreateSubDM(dmElastic, dM, "ELASTIC_PROBLEM");
560  CHKERR DMMoFEMSetDestroyProblem(dmElastic, PETSC_TRUE);
561  CHKERR DMMoFEMAddSubFieldRow(dmElastic, piolaStress);
562  CHKERR DMMoFEMAddSubFieldRow(dmElastic, bubbleField);
563  CHKERR DMMoFEMAddSubFieldRow(dmElastic, stretchTensor);
564  CHKERR DMMoFEMAddSubFieldRow(dmElastic, rotAxis);
565  CHKERR DMMoFEMAddSubFieldRow(dmElastic, spatialL2Disp);
566  CHKERR DMMoFEMAddSubFieldRow(dmElastic, contactDisp);
567  CHKERR DMMoFEMAddElement(dmElastic, elementVolumeName);
568  CHKERR DMMoFEMAddElement(dmElastic, naturalBcElement);
569  CHKERR DMMoFEMAddElement(dmElastic, essentialBcElement);
570  CHKERR DMMoFEMAddElement(dmElastic, contactElement);
571  CHKERR DMMoFEMAddElement(dmElastic, skinElement);
572  CHKERR DMMoFEMSetSquareProblem(dmElastic, PETSC_TRUE);
573  CHKERR DMMoFEMSetIsPartitioned(dmElastic, PETSC_TRUE);
574  CHKERR DMSetUp(dmElastic);
575 
576  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
577  "ELASTIC_PROBLEM", spatialL2Disp, stretchTensor);
578  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
579  "ELASTIC_PROBLEM", stretchTensor, spatialL2Disp);
580  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
581  "ELASTIC_PROBLEM", spatialL2Disp, rotAxis);
582  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
583  "ELASTIC_PROBLEM", rotAxis, spatialL2Disp);
584  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
585  "ELASTIC_PROBLEM", spatialL2Disp, bubbleField);
586  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
587  "ELASTIC_PROBLEM", bubbleField, spatialL2Disp);
588  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
589  "ELASTIC_PROBLEM", bubbleField, bubbleField);
590  // CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
591  // "ELASTIC_PROBLEM", piolaStress, piolaStress);
592  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
593  "ELASTIC_PROBLEM", bubbleField, piolaStress);
594  CHKERR mField.getInterface<ProblemsManager>()->addFieldToEmptyFieldBlocks(
595  "ELASTIC_PROBLEM", piolaStress, bubbleField);
596  {
597  PetscSection section;
598  CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
599  &section);
600  CHKERR DMSetSection(dmElastic, section);
601  CHKERR DMSetGlobalSection(dmElastic, section);
602  CHKERR PetscSectionDestroy(&section);
603  }
604 
605  dmPrjSpatial = createDM(mField.get_comm(), "DMMOFEM");
606  CHKERR DMMoFEMCreateSubDM(dmPrjSpatial, dM, "PROJECT_SPATIAL");
607  CHKERR DMMoFEMSetDestroyProblem(dmPrjSpatial, PETSC_TRUE);
608  CHKERR DMMoFEMAddSubFieldRow(dmPrjSpatial, spatialH1Disp);
609  CHKERR DMMoFEMAddElement(dmPrjSpatial, elementVolumeName);
610  CHKERR DMMoFEMSetSquareProblem(dmPrjSpatial, PETSC_TRUE);
611  CHKERR DMMoFEMSetIsPartitioned(dmPrjSpatial, PETSC_TRUE);
612  CHKERR DMSetUp(dmPrjSpatial);
613 
614  CHKERR mField.getInterface<BcManager>()
615  ->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
616  "PROJECT_SPATIAL", spatialH1Disp, true, false);
617 
618  MoFEMFunctionReturn(0);
619 }
620 
621 BcDisp::BcDisp(std::string name, std::vector<double> &attr, Range &faces)
622  : blockName(name), faces(faces) {
623  vals.resize(3, false);
624  flags.resize(3, false);
625  for (int ii = 0; ii != 3; ++ii) {
626  vals[ii] = attr[ii];
627  flags[ii] = static_cast<int>(attr[ii + 3]);
628  }
629 
630  MOFEM_LOG("EP", Sev::inform) << "Add BCDisp " << name;
631  MOFEM_LOG("EP", Sev::inform)
632  << "Add BCDisp vals " << vals[0] << " " << vals[1] << " " << vals[2];
633  MOFEM_LOG("EP", Sev::inform)
634  << "Add BCDisp flags " << flags[0] << " " << flags[1] << " " << flags[2];
635  MOFEM_LOG("EP", Sev::inform) << "Add BCDisp nb. of faces " << faces.size();
636 }
637 
638 BcRot::BcRot(std::string name, std::vector<double> &attr, Range &faces)
639  : blockName(name), faces(faces) {
640  vals.resize(3, false);
641  for (int ii = 0; ii != 3; ++ii) {
642  vals[ii] = attr[ii];
643  }
644  theta = attr[3];
645 }
646 
647 TractionBc::TractionBc(std::string name, std::vector<double> &attr,
648  Range &faces)
649  : blockName(name), faces(faces) {
650  vals.resize(3, false);
651  flags.resize(3, false);
652  for (int ii = 0; ii != 3; ++ii) {
653  vals[ii] = attr[ii];
654  flags[ii] = static_cast<int>(attr[ii + 3]);
655  }
656 
657  MOFEM_LOG("EP", Sev::inform) << "Add BCForce " << name;
658  MOFEM_LOG("EP", Sev::inform)
659  << "Add BCForce vals " << vals[0] << " " << vals[1] << " " << vals[2];
660  MOFEM_LOG("EP", Sev::inform)
661  << "Add BCForce flags " << flags[0] << " " << flags[1] << " " << flags[2];
662  MOFEM_LOG("EP", Sev::inform) << "Add BCForce nb. of faces " << faces.size();
663 }
664 
665 MoFEMErrorCode
666 EshelbianCore::getTractionFreeBc(const EntityHandle meshset,
667  boost::shared_ptr<TractionFreeBc> &bc_ptr,
668  const std::string contact_set_name) {
669  MoFEMFunctionBegin;
670 
671  // get skin from all tets
672  Range tets;
673  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
674  Range tets_skin_part;
675  Skinner skin(&mField.get_moab());
676  CHKERR skin.find_skin(0, tets, false, tets_skin_part);
677  ParallelComm *pcomm =
678  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
679  Range tets_skin;
680  CHKERR pcomm->filter_pstatus(tets_skin_part,
681  PSTATUS_SHARED | PSTATUS_MULTISHARED,
682  PSTATUS_NOT, -1, &tets_skin);
683 
684  bc_ptr->resize(3);
685  for (int dd = 0; dd != 3; ++dd)
686  (*bc_ptr)[dd] = tets_skin;
687 
688  if (bcSpatialDispVecPtr)
689  for (auto &v : *bcSpatialDispVecPtr) {
690  if (v.flags[0])
691  (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
692  if (v.flags[1])
693  (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
694  if (v.flags[2])
695  (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
696  }
697 
698  if (bcSpatialRotationVecPtr)
699  for (auto &v : *bcSpatialRotationVecPtr) {
700  (*bc_ptr)[0] = subtract((*bc_ptr)[0], v.faces);
701  (*bc_ptr)[1] = subtract((*bc_ptr)[1], v.faces);
702  (*bc_ptr)[2] = subtract((*bc_ptr)[2], v.faces);
703  }
704 
705  // remove contact
706  for (auto m : mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
707  std::regex((boost::format("%s(.*)") % contact_set_name).str()))) {
708  Range faces;
709  CHKERR m->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
710  true);
711  (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
712  (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
713  (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
714  }
715 
716  // for (int dd = 0; dd != 3; ++dd) {
717  // EntityHandle meshset;
718  // CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
719  // CHKERR mField.get_moab().add_entities(meshset, (*bc_ptr)[dd]);
720  // std::string file_name = disp_block_set_name +
721  // "_traction_free_bc_" + boost::lexical_cast<std::string>(dd) + ".vtk";
722  // CHKERR mField.get_moab().write_file(file_name.c_str(), " VTK ", "",
723  // &meshset, 1);
724  // CHKERR mField.get_moab().delete_entities(&meshset, 1);
725  // }
726 
727  MoFEMFunctionReturn(0);
728 }
729 
730 /**
731  * @brief Set integration rule on element
732  * \param order on row
733  * \param order on column
734  * \param order on data
735  *
736  * Use maximal oder on data in order to determine integration rule
737  *
738  */
739 struct VolRule {
740  int operator()(int p_row, int p_col, int p_data) const {
741  if (EshelbianCore::gradApperoximator > MODERATE_ROT)
742  return p_data + p_data + (p_data - 1);
743  else
744  return 2 * p_data;
745  }
746 };
747 
748 struct FaceRule {
749  int operator()(int p_row, int p_col, int p_data) const { return 2 * p_data; }
750 };
751 
752 struct CGGUserPolynomialBase : public BaseFunction {
753 
754  CGGUserPolynomialBase() {}
755  ~CGGUserPolynomialBase() {}
756 
757  MoFEMErrorCode query_interface(boost::typeindex::type_index type_index,
758  BaseFunctionUnknownInterface **iface) const {
759  *iface = const_cast<CGGUserPolynomialBase *>(this);
760  return 0;
761  }
762 
763  MoFEMErrorCode getValue(MatrixDouble &pts,
764  boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
765  MoFEMFunctionBeginHot;
766 
767  cTx = ctx_ptr->getInterface<EntPolynomialBaseCtx>();
768 
769  int nb_gauss_pts = pts.size2();
770  if (!nb_gauss_pts) {
771  MoFEMFunctionReturnHot(0);
772  }
773 
774  if (pts.size1() < 3) {
775  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
776  "Wrong dimension of pts, should be at least 3 rows with "
777  "coordinates");
778  }
779 
780  switch (cTx->sPace) {
781  case HDIV:
782  CHKERR getValueHdivForCGGBubble(pts);
783  break;
784  default:
785  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
786  }
787 
788  MoFEMFunctionReturnHot(0);
789  }
790 
791 private:
792  EntPolynomialBaseCtx *cTx;
793 
794  MatrixDouble shapeFun;
795  boost::tuple<int, int, MatrixDouble> cachePhi = {0, 0, MatrixDouble()};
796 
797  MoFEMErrorCode getValueHdivForCGGBubble(MatrixDouble &pts) {
798  MoFEMFunctionBegin;
799 
800  // This should be used only in case USER_BASE is selected
801  if (cTx->bAse != USER_BASE) {
802  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
803  "Wrong base, should be USER_BASE");
804  }
805  // get access to data structures on element
806  EntitiesFieldData &data = cTx->dAta;
807  // Get approximation order on element. Note that bubble functions are only
808  // on tetrahedron.
809  const int order = data.dataOnEntities[MBTET][0].getOrder();
810  /// number of integration points
811  const int nb_gauss_pts = pts.size2();
812 
813  if(cachePhi.get<0>() == order && cachePhi.get<1>() == nb_gauss_pts) {
814  auto &mat = cachePhi.get<2>();
815  auto &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
816  phi.resize(mat.size1(), mat.size2(), false);
817  noalias(phi) = mat;
818  } else {
819  // calculate shape functions, i.e. barycentric coordinates
820  shapeFun.resize(nb_gauss_pts, 4, false);
821  CHKERR ShapeMBTET(&*shapeFun.data().begin(), &pts(0, 0), &pts(1, 0),
822  &pts(2, 0), nb_gauss_pts);
823  // direvatives of shape functions
824  double diff_shape_fun[12];
825  CHKERR ShapeDiffMBTET(diff_shape_fun);
826 
827  const int nb_base_functions = NBVOLUMETET_CCG_BUBBLE(order);
828  // get base functions and set size
829  MatrixDouble &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
830  phi.resize(nb_gauss_pts, 9 * nb_base_functions, false);
831  // finally calculate base functions
832  FTensor::Tensor2<FTensor::PackPtr<double *, 9>, 3, 3> t_phi(
833  &phi(0, 0), &phi(0, 1), &phi(0, 2),
834 
835  &phi(0, 3), &phi(0, 4), &phi(0, 5),
836 
837  &phi(0, 6), &phi(0, 7), &phi(0, 8));
838  CHKERR CGG_BubbleBase_MBTET(order, &shapeFun(0, 0), diff_shape_fun, t_phi,
839  nb_gauss_pts);
840 
841  cachePhi.get<0>() = order;
842  cachePhi.get<1>() = nb_gauss_pts;
843  cachePhi.get<2>().resize(phi.size1(), phi.size2(), false);
844  noalias(cachePhi.get<2>()) = phi;
845  }
846 
847  MoFEMFunctionReturn(0);
848  }
849 
850 };
851 
852 MoFEMErrorCode EshelbianCore::setBaseVolumeElementOps(
853  const int tag, const bool do_rhs, const bool do_lhs,
854  boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe) {
855  MoFEMFunctionBegin;
856  fe = boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
857 
858  fe->getUserPolynomialBase() =
859  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
860  AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(fe->getOpPtrVector(),
861  {HDIV, H1, L2});
862 
863  // set integration rule
864  fe->getRuleHook = VolRule();
865 
866  if (!dataAtPts) {
867  dataAtPts =
868  boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
869  dataAtPts->physicsPtr = physicalEquations;
870  }
871 
872  // calculate fields values
873  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
874  piolaStress, dataAtPts->getApproxPAtPts()));
875  fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
876  bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
877  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
878  piolaStress, dataAtPts->getDivPAtPts()));
879  fe->getOpPtrVector().push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
880  stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
881  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
882  rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
883  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
884  spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
885 
886  // velocities
887  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
888  spatialL2Disp, dataAtPts->getSmallWL2DotAtPts(), MBTET));
889  fe->getOpPtrVector().push_back(
890  new OpCalculateTensor2SymmetricFieldValuesDot<3>(
891  stretchTensor, dataAtPts->getLogStretchDotTensorAtPts(), MBTET));
892  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
893  rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
894 
895  // acceleration
896  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
897  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDotDot<3>(
898  spatialL2Disp, dataAtPts->getSmallWL2DotDotAtPts(), MBTET));
899  }
900 
901  // H1 displacements
902  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
903  spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
904  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
905  spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
906 
907  // calculate other derived quantities
908  fe->getOpPtrVector().push_back(
909  new OpCalculateRotationAndSpatialGradient(dataAtPts));
910 
911  // evaluate integration points
912  fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
913  tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
914 
915  MoFEMFunctionReturn(0);
916 }
917 
918 MoFEMErrorCode EshelbianCore::setVolumeElementOps(
919  const int tag, const bool add_elastic, const bool add_material,
920  boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_rhs,
921  boost::shared_ptr<VolumeElementForcesAndSourcesCore> &fe_lhs) {
922  MoFEMFunctionBegin;
923 
924  auto time_scale = boost::make_shared<TimeScale>();
925 
926  // Right hand side
927  CHKERR setBaseVolumeElementOps(tag, true, false, fe_rhs);
928 
929  // elastic
930  if (add_elastic) {
931  fe_rhs->getOpPtrVector().push_back(
932  new OpSpatialEquilibrium(spatialL2Disp, dataAtPts, alphaW, alphaRho));
933  fe_rhs->getOpPtrVector().push_back(
934  new OpSpatialRotation(rotAxis, dataAtPts));
935  fe_rhs->getOpPtrVector().push_back(
936  new OpSpatialPhysical(stretchTensor, dataAtPts, alphaU));
937  fe_rhs->getOpPtrVector().push_back(
938  new OpSpatialConsistencyP(piolaStress, dataAtPts));
939  fe_rhs->getOpPtrVector().push_back(
940  new OpSpatialConsistencyBubble(bubbleField, dataAtPts));
941  fe_rhs->getOpPtrVector().push_back(
942  new OpSpatialConsistencyDivTerm(piolaStress, dataAtPts));
943 
944  // Body forces
945  using BodyNaturalBC =
946  NaturalBC<VolumeElementForcesAndSourcesCore::UserDataOperator>::
947  Assembly<PETSC>::LinearForm<GAUSS>;
948  using OpBodyForce =
949  BodyNaturalBC::OpFlux<NaturalMeshsetType<BLOCKSET>, 1, 3>;
950  CHKERR BodyNaturalBC::AddFluxToPipeline<OpBodyForce>::add(
951  fe_rhs->getOpPtrVector(), mField, "w", {time_scale}, "BODY_FORCE",
952  Sev::inform);
953  }
954 
955  // Left hand side
956  CHKERR setBaseVolumeElementOps(tag, true, true, fe_lhs);
957 
958  // elastic
959  if (add_elastic) {
960 
961  const bool symmetric_system =
962  (gradApperoximator <= MODERATE_ROT && rotSelector == SMALL_ROT);
963 
964  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_du(
965  stretchTensor, stretchTensor, dataAtPts, alphaU));
966  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
967  stretchTensor, piolaStress, dataAtPts, true));
968  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
969  stretchTensor, bubbleField, dataAtPts, true));
970  if (!symmetric_system) {
971  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
972  stretchTensor, rotAxis, dataAtPts, false));
973  }
974 
975  fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
976  spatialL2Disp, piolaStress, dataAtPts, true));
977  fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
978  spatialL2Disp, spatialL2Disp, dataAtPts, alphaW, alphaRho));
979 
980  fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dP_domega(
981  piolaStress, rotAxis, dataAtPts, symmetric_system));
982  fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
983  bubbleField, rotAxis, dataAtPts, symmetric_system));
984 
985  if (!symmetric_system) {
986  fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dP(
987  rotAxis, piolaStress, dataAtPts, false));
988  fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dBubble(
989  rotAxis, bubbleField, dataAtPts, false));
990  fe_lhs->getOpPtrVector().push_back(
991  new OpSpatialRotation_domega_domega(rotAxis, rotAxis, dataAtPts));
992  }
993 
994  // Stabilisation
995  if constexpr (A == AssemblyType::SCHUR) {
996  // Note that we assemble to AMat, however Assembly<SCHUR> assemble by
997  // default to PMat
998  using OpMassStab =
999  FormsIntegrators<VolUserDataOperatorStabAssembly>::Assembly<
1000  SCHUR>::BiLinearForm<GAUSS>::OpMass<1, SPACE_DIM>;
1001  if (precEpsOmega > std::numeric_limits<double>::epsilon()) {
1002  fe_lhs->getOpPtrVector().push_back(
1003  new OpMassStab(rotAxis, rotAxis, [this](double, double, double) {
1004  return precEpsOmega;
1005  }));
1006  }
1007  if (std::abs(alphaRho + alphaW) <
1008  std::numeric_limits<double>::epsilon()) {
1009 
1010  fe_lhs->getOpPtrVector().push_back(new OpMassStab(
1011  spatialL2Disp, spatialL2Disp,
1012  [this](double, double, double) { return precEpsW; }));
1013  }
1014  }
1015 
1016  if (add_material) {
1017  }
1018  }
1019 
1020  MoFEMFunctionReturn(0);
1021 }
1022 
1023 MoFEMErrorCode EshelbianCore::setFaceElementOps(
1024  const bool add_elastic, const bool add_material,
1025  boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_rhs,
1026  boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_lhs) {
1027  MoFEMFunctionBegin;
1028 
1029  fe_rhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
1030  fe_lhs = boost::make_shared<FaceElementForcesAndSourcesCore>(mField);
1031 
1032  // set integration rule
1033  fe_rhs->getRuleHook = FaceRule();
1034  fe_lhs->getRuleHook = FaceRule();
1035 
1036  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(fe_rhs->getOpPtrVector(),
1037  {HDIV});
1038  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(fe_lhs->getOpPtrVector(),
1039  {HDIV});
1040 
1041  if (add_elastic) {
1042 
1043  fe_rhs->getOpPtrVector().push_back(
1044  new OpDispBc(piolaStress, dataAtPts, bcSpatialDispVecPtr,
1045  {
1046 
1047  boost::make_shared<TimeScale>("disp_history.txt")
1048 
1049  }));
1050  fe_rhs->getOpPtrVector().push_back(
1051  new OpRotationBc(piolaStress, dataAtPts, bcSpatialRotationVecPtr));
1052  }
1053 
1054  MoFEMFunctionReturn(0);
1055 }
1056 
1057 MoFEMErrorCode EshelbianCore::setContactElementOps(
1058 
1059  boost::shared_ptr<ContactTree> &fe_contact_tree,
1060  boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_face_rhs,
1061  boost::shared_ptr<FaceElementForcesAndSourcesCore> &fe_face_lhs
1062 
1063 ) {
1064  MoFEMFunctionBegin;
1065 
1066  /** Contact requires that body is marked */
1067  auto get_body_range = [this](auto name, int dim) {
1068  std::map<int, Range> map;
1069 
1070  for (auto m_ptr :
1071  mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
1072 
1073  (boost::format("%s(.*)") % name).str()
1074 
1075  ))
1076 
1077  ) {
1078  Range ents;
1079  CHK_MOAB_THROW(m_ptr->getMeshsetIdEntitiesByDimension(mField.get_moab(),
1080  dim, ents, true),
1081  "by dim");
1082  map[m_ptr->getMeshsetId()] = ents;
1083  }
1084 
1085  return map;
1086  };
1087 
1088  auto get_map_skin = [this](auto &&map) {
1089  ParallelComm *pcomm =
1090  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
1091 
1092  Skinner skin(&mField.get_moab());
1093  for(auto &m : map) {
1094  Range skin_faces;
1095  CHKERR skin.find_skin(0, m.second, false, skin_faces);
1096  CHK_MOAB_THROW(pcomm->filter_pstatus(skin_faces,
1097  PSTATUS_SHARED | PSTATUS_MULTISHARED,
1098  PSTATUS_NOT, -1, nullptr),
1099  "filter");
1100  m.second.swap(skin_faces);
1101  }
1102  return map;
1103  };
1104 
1105  /* The above code is written in C++ and it appears to be defining and using
1106  various operations on boundary elements and side elements. */
1107  using BoundaryEle = FaceElementForcesAndSourcesCore;
1108  using BoundaryEleOp = BoundaryEle::UserDataOperator;
1109 
1110  fe_face_rhs = boost::make_shared<BoundaryEle>(mField);
1111  fe_face_lhs = boost::make_shared<BoundaryEle>(mField);
1112  auto rule = [](int, int, int o) { return -1; };
1113 
1114  int levels = 0;
1115  CHKERR PetscOptionsGetInt(PETSC_NULL, "-contact_max_post_proc_ref_level",
1116  &levels, PETSC_NULL);
1117 
1118  auto refine = Tools::refineTriangle(levels);
1119 
1120  auto set_rule = [levels, refine](
1121 
1122  ForcesAndSourcesCore *fe_raw_ptr, int order_row,
1123  int order_col, int order_data
1124 
1125  ) {
1126  MoFEMFunctionBegin;
1127  auto rule = 2 * order_data;
1128  fe_raw_ptr->gaussPts = Tools::refineTriangleIntegrationPts(rule, refine);
1129  MoFEMFunctionReturn(0);
1130  };
1131 
1132  fe_face_rhs->getRuleHook = rule;
1133  fe_face_lhs->getRuleHook = rule;
1134  fe_face_rhs->setRuleHook = set_rule;
1135  fe_face_lhs->setRuleHook = set_rule;
1136 
1137  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
1138  fe_face_rhs->getOpPtrVector(), {HDIV});
1139  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
1140  fe_face_lhs->getOpPtrVector(), {HDIV});
1141 
1142  auto add_contact_three = [&]() {
1143  MoFEMFunctionBegin;
1144  auto tree_moab_ptr = boost::make_shared<moab::Core>();
1145  fe_contact_tree = boost::make_shared<ContactTree>(
1146  mField, tree_moab_ptr, spaceOrder,
1147  get_map_skin(get_body_range("BODY", 3)));
1148  auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1149  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
1150  fe_contact_tree->getOpPtrVector(), {HDIV});
1151  fe_contact_tree->getOpPtrVector().push_back(
1152  new OpCalculateVectorFieldValues<3>(
1153  contactDisp, contact_common_data_ptr->contactDispPtr()));
1154  fe_contact_tree->getOpPtrVector().push_back(
1155  new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
1156  piolaStress, contact_common_data_ptr->contactTractionPtr()));
1157 
1158  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
1159  fe_contact_tree->getOpPtrVector().push_back(
1160  new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
1161  fe_contact_tree->getOpPtrVector().push_back(
1162  new OpMoveNode(fe_contact_tree, contact_common_data_ptr, u_h1_ptr));
1163  MoFEMFunctionReturn(0);
1164  };
1165 
1166  auto add_ops_face_lhs = [&](auto &pip) {
1167  MoFEMFunctionBegin;
1168  auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1169  pip.push_back(new OpCalculateVectorFieldValues<3>(
1170  contactDisp, contact_common_data_ptr->contactDispPtr()));
1171  pip.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
1172  piolaStress, contact_common_data_ptr->contactTractionPtr()));
1173  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
1174  pip.push_back(new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
1175  pip.push_back(new OpTreeSearch(
1176  fe_contact_tree, contact_common_data_ptr, u_h1_ptr,
1177  get_range_from_block(mField, "CONTACT"), nullptr, nullptr));
1178 
1179  auto contact_sfd_map_range_ptr = boost::make_shared<std::map<int, Range>>(
1180  get_body_range("CONTACT_SDF", 2));
1181 
1182  pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dU(
1183  contactDisp, contactDisp, contact_common_data_ptr, fe_contact_tree,
1184  contact_sfd_map_range_ptr));
1185  pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dP(
1186  contactDisp, piolaStress, contact_common_data_ptr, fe_contact_tree,
1187  contact_sfd_map_range_ptr));
1188  pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryHDivLhs_dU(
1189  piolaStress, contactDisp, contact_common_data_ptr, fe_contact_tree));
1190 
1191  // Stabilisation
1192  if constexpr (A == AssemblyType::SCHUR) {
1193  // Note that we assemble to AMat, however Assembly<SCHUR> assemble by
1194  // default to PMat
1195  using OpMassStab =
1196  FormsIntegrators<FaceUserDataOperatorStabAssembly>::Assembly<
1197  SCHUR>::BiLinearForm<GAUSS>::OpMass<1, SPACE_DIM>;
1198  if (precEpsContactDisp > std::numeric_limits<double>::epsilon()) {
1199  pip.push_back(new OpMassStab(
1200  contactDisp, contactDisp,
1201  [this](double, double, double) { return precEpsContactDisp; }));
1202  }
1203  }
1204 
1205  MoFEMFunctionReturn(0);
1206  };
1207 
1208  auto add_ops_face_rhs = [&](auto &pip) {
1209  MoFEMFunctionBegin;
1210  auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1211  pip.push_back(new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
1212  piolaStress, contact_common_data_ptr->contactTractionPtr()));
1213  pip.push_back(new OpCalculateVectorFieldValues<3>(
1214  contactDisp, contact_common_data_ptr->contactDispPtr()));
1215  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
1216  pip.push_back(new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
1217  pip.push_back(new OpTreeSearch(
1218  fe_contact_tree, contact_common_data_ptr, u_h1_ptr,
1219  get_range_from_block(mField, "CONTACT"), nullptr, nullptr));
1220  auto contact_sfd_map_range_ptr = boost::make_shared<std::map<int, Range>>(
1221  get_body_range("CONTACT_SDF", 2));
1222  pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryL2Rhs(
1223  contactDisp, contact_common_data_ptr, fe_contact_tree,
1224  contact_sfd_map_range_ptr));
1225  pip.push_back(new EshelbianPlasticity::OpConstrainBoundaryHDivRhs(
1226  piolaStress, contact_common_data_ptr, fe_contact_tree));
1227  MoFEMFunctionReturn(0);
1228  };
1229 
1230 
1231 
1232  CHKERR add_contact_three();
1233  CHKERR add_ops_face_lhs(fe_face_lhs->getOpPtrVector());
1234  CHKERR add_ops_face_rhs(fe_face_rhs->getOpPtrVector());
1235 
1236  MoFEMFunctionReturn(0);
1237 }
1238 
1239 MoFEMErrorCode EshelbianCore::setElasticElementOps(const int tag) {
1240  MoFEMFunctionBegin;
1241  CHKERR setVolumeElementOps(tag, true, false, elasticFeRhs, elasticFeLhs);
1242  CHKERR setFaceElementOps(true, false, elasticBcRhs, elasticBcLhs);
1243 
1244  auto adj_cache =
1245  boost::make_shared<ForcesAndSourcesCore::UserDataOperator::AdjCache>();
1246 
1247  auto get_op_contact_bc = [&]() {
1248  using SideEle = FaceElementForcesAndSourcesCore;
1249  auto op_loop_side = new OpLoopSide<SideEle>(
1250  mField, contactElement, SPACE_DIM - 1, Sev::noisy, adj_cache);
1251  return op_loop_side;
1252  };
1253 
1254  auto op_contact_bc = get_op_contact_bc();
1255  elasticFeLhs->getOpPtrVector().push_back(op_contact_bc);
1256 
1257  CHKERR setContactElementOps(
1258 
1259  contactTreeRhs, contactRhs, op_contact_bc->getSideFEPtr()
1260 
1261  );
1262  MoFEMFunctionReturn(0);
1263 }
1264 
1265 MoFEMErrorCode EshelbianCore::setElasticElementToTs(DM dm) {
1266  MoFEMFunctionBegin;
1267  boost::shared_ptr<FEMethod> null;
1268  boost::shared_ptr<FeTractionBc> spatial_traction_bc(
1269  new FeTractionBc(mField, piolaStress, bcSpatialTraction));
1270 
1271  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
1272 
1273  CHKERR DMMoFEMTSSetI2Function(dm, DM_NO_ELEMENT, null, spatial_traction_bc,
1274  null);
1275  CHKERR DMMoFEMTSSetI2Function(dm, elementVolumeName, elasticFeRhs, null,
1276  null);
1277  CHKERR DMMoFEMTSSetI2Function(dm, naturalBcElement, elasticBcRhs, null,
1278  null);
1279  CHKERR DMMoFEMTSSetI2Function(dm, contactElement, contactRhs, null, null);
1280  CHKERR DMMoFEMTSSetI2Function(dm, DM_NO_ELEMENT, null, null,
1281  spatial_traction_bc);
1282  CHKERR DMMoFEMTSSetI2Jacobian(dm, elementVolumeName, elasticFeLhs, null,
1283  null);
1284  CHKERR DMMoFEMTSSetI2Jacobian(dm, naturalBcElement, elasticBcLhs, null,
1285  null);
1286 
1287  } else {
1288  CHKERR DMMoFEMTSSetIFunction(dm, DM_NO_ELEMENT, null, spatial_traction_bc,
1289  null);
1290  CHKERR DMMoFEMTSSetIFunction(dm, elementVolumeName, elasticFeRhs, null,
1291  null);
1292  CHKERR DMMoFEMTSSetIFunction(dm, naturalBcElement, elasticBcRhs, null,
1293  null);
1294  CHKERR DMMoFEMTSSetIFunction(dm, contactElement, contactRhs, null, null);
1295  CHKERR DMMoFEMTSSetIFunction(dm, DM_NO_ELEMENT, null, null,
1296  spatial_traction_bc);
1297  CHKERR DMMoFEMTSSetIJacobian(dm, elementVolumeName, elasticFeLhs, null,
1298  null);
1299  CHKERR DMMoFEMTSSetIJacobian(dm, naturalBcElement, elasticBcLhs, null,
1300  null);
1301  }
1302  MoFEMFunctionReturn(0);
1303 }
1304 
1305 #include "impl/EshelbianMonitor.cpp"
1306 #include "impl/TSElasticPostStep.cpp"
1307 #include "impl/SetUpSchurImpl.cpp"
1308 
1309 MoFEMErrorCode EshelbianCore::solveElastic(TS ts, Mat m, Vec f, Vec x) {
1310  MoFEMFunctionBegin;
1311 
1312 #ifdef PYTHON_SDF
1313 
1314  boost::shared_ptr<ContactOps::SDFPython> sdf_python_ptr;
1315 
1316  auto file_exists = [](std::string myfile) {
1317  std::ifstream file(myfile.c_str());
1318  if (file) {
1319  return true;
1320  }
1321  return false;
1322  };
1323 
1324  if (file_exists("sdf.py")) {
1325  MOFEM_LOG("EP", Sev::inform) << "sdf.py file found";
1326  sdf_python_ptr = boost::make_shared<ContactOps::SDFPython>();
1327  CHKERR sdf_python_ptr->sdfInit("sdf.py");
1328  ContactOps::sdfPythonWeakPtr = sdf_python_ptr;
1329  } else {
1330  MOFEM_LOG("EP", Sev::warning) << "sdf.py file NOT found";
1331  }
1332 #else
1333 #endif
1334 
1335  boost::shared_ptr<TsCtx> ts_ctx;
1336  CHKERR DMMoFEMGetTsCtx(dmElastic, ts_ctx);
1337 
1338  CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx.get(), PETSC_NULL);
1339 
1340  boost::shared_ptr<FEMethod> monitor_ptr(new EshelbianMonitor(*this));
1341  ts_ctx->getLoopsMonitor().push_back(
1342  TsCtx::PairNameFEMethodPtr(elementVolumeName, monitor_ptr));
1343 
1344  CHKERR TSAppendOptionsPrefix(ts, "elastic_");
1345  CHKERR TSSetFromOptions(ts);
1346 
1347  CHKERR TSSetDM(ts, dmElastic);
1348 
1349  SNES snes;
1350  CHKERR TSGetSNES(ts, &snes);
1351 
1352  PetscViewerAndFormat *vf;
1353  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
1354  PETSC_VIEWER_DEFAULT, &vf);
1355  CHKERR SNESMonitorSet(
1356  snes,
1357  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
1358  vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
1359 
1360  PetscSection section;
1361  CHKERR DMGetSection(dmElastic, &section);
1362  int num_fields;
1363  CHKERR PetscSectionGetNumFields(section, &num_fields);
1364  for (int ff = 0; ff != num_fields; ff++) {
1365  const char *field_name;
1366  CHKERR PetscSectionGetFieldName(section, ff, &field_name);
1367  MOFEM_LOG_C("EP", Sev::inform, "Field %d name %s", ff, field_name);
1368  }
1369 
1370  CHKERR DMoFEMMeshToLocalVector(dmElastic, x, INSERT_VALUES, SCATTER_FORWARD);
1371 
1372  CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
1373  CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
1374 
1375  // Adding field split solver
1376  boost::shared_ptr<SetUpSchur> schur_ptr;
1377 
1378  if constexpr (A == AssemblyType::SCHUR) {
1379 
1380  auto p = matDuplicate(m, MAT_DO_NOT_COPY_VALUES);
1381 
1382  auto ts_ctx_ptr = getDMTsCtx(dmElastic);
1383  ts_ctx_ptr->zeroMatrix = false;
1384 
1385  // If density is larger than zero, use dynamic time solver
1386  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
1387  CHKERR TSSetI2Function(ts, f, PETSC_NULL, PETSC_NULL);
1388  CHKERR TSSetI2Jacobian(ts, m, p, PETSC_NULL, PETSC_NULL);
1389  } else {
1390  CHKERR TSSetIFunction(ts, f, PETSC_NULL, PETSC_NULL);
1391  CHKERR TSSetIJacobian(ts, m, p, PETSC_NULL, PETSC_NULL);
1392  }
1393 
1394  KSP ksp;
1395  CHKERR SNESGetKSP(snes, &ksp);
1396  PC pc;
1397  CHKERR KSPGetPC(ksp, &pc);
1398  schur_ptr = SetUpSchur::createSetUpSchur(
1399  mField, SmartPetscObj<Mat>(m, true), p, &*this);
1400  CHKERR schur_ptr->setUp(ksp);
1401 
1402  elasticFeLhs->preProcessHook = [&]() {
1403  MoFEMFunctionBegin;
1404  *(elasticFeLhs->matAssembleSwitch) = false;
1405  CHKERR schur_ptr->preProc();
1406  MoFEMFunctionReturn(0);
1407  };
1408  elasticFeLhs->postProcessHook = [&]() {
1409  MoFEMFunctionBegin;
1410  MoFEMFunctionReturn(0);
1411  };
1412  elasticBcLhs->preProcessHook = [&]() {
1413  MoFEMFunctionBegin;
1414  MoFEMFunctionReturn(0);
1415  };
1416  elasticBcLhs->preProcessHook = [&]() {
1417  MoFEMFunctionBegin;
1418  *(elasticBcLhs->matAssembleSwitch) = false;
1419  CHKERR schur_ptr->postProc();
1420  MoFEMFunctionReturn(0);
1421  };
1422  }
1423 
1424  if (std::abs(alphaRho) > std::numeric_limits<double>::epsilon()) {
1425  Vec xx;
1426  CHKERR VecDuplicate(x, &xx);
1427  CHKERR VecZeroEntries(xx);
1428  CHKERR TS2SetSolution(ts, x, xx);
1429  CHKERR VecDestroy(&xx);
1430  } else {
1431  CHKERR TSSetSolution(ts, x);
1432  }
1433 
1434  CHKERR TSSetUp(ts);
1435  CHKERR TSSetPreStep(ts, TSElasticPostStep::preStepFun);
1436  CHKERR TSSetPostStep(ts, TSElasticPostStep::postStepFun);
1437  CHKERR TSElasticPostStep::postStepInitialise(this);
1438  CHKERR TSSolve(ts, PETSC_NULL);
1439  CHKERR TSElasticPostStep::postStepDestroy();
1440 
1441  // CHKERR TSGetSNES(ts, &snes);
1442  int lin_solver_iterations;
1443  CHKERR SNESGetLinearSolveIterations(snes, &lin_solver_iterations);
1444  MOFEM_LOG("EP", Sev::inform)
1445  << "Number of linear solver iterations " << lin_solver_iterations;
1446 
1447  PetscBool test_cook_flg = PETSC_FALSE;
1448  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_cook", &test_cook_flg,
1449  PETSC_NULL);
1450  if (test_cook_flg) {
1451  constexpr int expected_lin_solver_iterations = 11;
1452  // if (lin_solver_iterations != expected_lin_solver_iterations)
1453  // SETERRQ2(
1454  // PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
1455  // "Expected number of iterations is different than expected %d !=
1456  // %d", lin_solver_iterations, expected_lin_solver_iterations);
1457  }
1458 
1459  CHKERR gettingNorms();
1460 
1461  MoFEMFunctionReturn(0);
1462 }
1463 
1464 
1465 MoFEMErrorCode EshelbianCore::postProcessResults(const int tag,
1466  const std::string file) {
1467  MoFEMFunctionBegin;
1468 
1469  if (!dataAtPts) {
1470  dataAtPts =
1471  boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
1472  }
1473 
1474  PostProcBrokenMeshInMoab<FaceElementForcesAndSourcesCore> post_proc(mField);
1475  AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(post_proc.getOpPtrVector(),
1476  {HDIV});
1477 
1478  auto contact_common_data_ptr = boost::make_shared<ContactOps::CommonData>();
1479 
1480  auto domain_ops = [&](auto &fe) {
1481  MoFEMFunctionBegin;
1482  fe.getUserPolynomialBase() =
1483  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1484  AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(fe.getOpPtrVector(),
1485  {HDIV, H1, L2});
1486  fe.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
1487  piolaStress, dataAtPts->getApproxPAtPts()));
1488  fe.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
1489  bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
1490  fe.getOpPtrVector().push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
1491  stretchTensor, dataAtPts->getLogStretchTensorAtPts(), MBTET));
1492  fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1493  rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
1494  fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1495  spatialL2Disp, dataAtPts->getSmallWL2AtPts(), MBTET));
1496  fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1497  spatialH1Disp, dataAtPts->getSmallWH1AtPts()));
1498  fe.getOpPtrVector().push_back(new OpCalculateVectorFieldGradient<3, 3>(
1499  spatialH1Disp, dataAtPts->getSmallWGradH1AtPts()));
1500  // evaluate derived quantities
1501  fe.getOpPtrVector().push_back(
1502  new OpCalculateRotationAndSpatialGradient(dataAtPts));
1503 
1504  // evaluate integration points
1505  fe.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
1506  tag, true, false, dataAtPts, physicalEquations));
1507 
1508  // contact
1509  fe.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
1510  spatialL2Disp, contact_common_data_ptr->contactDispPtr()));
1511 
1512  MoFEMFunctionReturn(0);
1513  };
1514 
1515  auto op_loop_side = new OpLoopSide<VolumeElementForcesAndSourcesCoreOnSide>(
1516  mField, elementVolumeName, SPACE_DIM);
1517  CHKERR domain_ops(*(op_loop_side->getSideFEPtr()));
1518  post_proc.getOpPtrVector().push_back(op_loop_side);
1519  post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1520  contactDisp, dataAtPts->getContactL2AtPts()));
1521  post_proc.getOpPtrVector().push_back(new OpPostProcDataStructure(
1522  post_proc.getPostProcMesh(), post_proc.getMapGaussPts(), dataAtPts));
1523 
1524  // contact
1525  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
1526  auto lambda_h1_ptr = boost::make_shared<MatrixDouble>();
1527  post_proc.getOpPtrVector().push_back(
1528  new OpCalculateVectorFieldValues<3>(spatialH1Disp, u_h1_ptr));
1529 
1530  using BoundaryEle = FaceElementForcesAndSourcesCore;
1531  using BoundaryEleOp = BoundaryEle::UserDataOperator;
1532 
1533  post_proc.getOpPtrVector().push_back(
1534  new OpCalculateHVecTensorTrace<SPACE_DIM, BoundaryEleOp>(
1535  piolaStress, contact_common_data_ptr->contactTractionPtr()));
1536 
1537  post_proc.getOpPtrVector().push_back(new OpTreeSearch(
1538  contactTreeRhs, contact_common_data_ptr, u_h1_ptr,
1539  get_range_from_block(mField, "CONTACT"), &post_proc.getPostProcMesh(),
1540  &post_proc.getMapGaussPts()));
1541 
1542  CHKERR DMoFEMLoopFiniteElements(dM, contactElement, contactTreeRhs);
1543  CHKERR DMoFEMLoopFiniteElements(dM, skinElement.c_str(), &post_proc);
1544  CHKERR post_proc.writeFile(file.c_str());
1545  MoFEMFunctionReturn(0);
1546 }
1547 
1548 //! [Getting norms]
1549 MoFEMErrorCode EshelbianCore::gettingNorms() {
1550  MoFEMFunctionBegin;
1551 
1552  auto post_proc_norm_fe =
1553  boost::make_shared<VolumeElementForcesAndSourcesCore>(mField);
1554 
1555  auto post_proc_norm_rule_hook = [](int, int, int p) -> int {
1556  return 2 * (p);
1557  };
1558  post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
1559 
1560  post_proc_norm_fe->getUserPolynomialBase() =
1561  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1562 
1563  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
1564  post_proc_norm_fe->getOpPtrVector(), {L2, H1, HDIV});
1565 
1566  enum NORMS { U_NORM_L2 = 0, U_NORM_H1, PIOLA_NORM, U_ERROR_L2, LAST_NORM };
1567  auto norms_vec =
1568  createVectorMPI(mField.get_comm(), LAST_NORM, PETSC_DETERMINE);
1569  CHKERR VecZeroEntries(norms_vec);
1570 
1571  auto u_l2_ptr = boost::make_shared<MatrixDouble>();
1572  auto u_h1_ptr = boost::make_shared<MatrixDouble>();
1573  post_proc_norm_fe->getOpPtrVector().push_back(
1574  new OpCalculateVectorFieldValues<SPACE_DIM>(spatialL2Disp, u_l2_ptr));
1575  post_proc_norm_fe->getOpPtrVector().push_back(
1576  new OpCalculateVectorFieldValues<SPACE_DIM>(spatialH1Disp, u_h1_ptr));
1577  post_proc_norm_fe->getOpPtrVector().push_back(
1578  new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_NORM_L2));
1579  post_proc_norm_fe->getOpPtrVector().push_back(
1580  new OpCalcNormL2Tensor1<SPACE_DIM>(u_h1_ptr, norms_vec, U_NORM_H1));
1581  post_proc_norm_fe->getOpPtrVector().push_back(
1582  new OpCalcNormL2Tensor1<SPACE_DIM>(u_l2_ptr, norms_vec, U_ERROR_L2,
1583  u_h1_ptr));
1584 
1585  auto piola_ptr = boost::make_shared<MatrixDouble>();
1586  post_proc_norm_fe->getOpPtrVector().push_back(
1587  new OpCalculateHVecTensorField<3, 3>(piolaStress, piola_ptr));
1588  post_proc_norm_fe->getOpPtrVector().push_back(
1589  new OpCalcNormL2Tensor2<3, 3>(piola_ptr, norms_vec, PIOLA_NORM));
1590 
1591  CHKERR mField.loop_finite_elements("ELASTIC_PROBLEM", elementVolumeName,
1592  *post_proc_norm_fe);
1593 
1594  CHKERR VecAssemblyBegin(norms_vec);
1595  CHKERR VecAssemblyEnd(norms_vec);
1596  const double *norms;
1597  CHKERR VecGetArrayRead(norms_vec, &norms);
1598  MOFEM_LOG("EP", Sev::inform) << "norm_u: " << std::sqrt(norms[U_NORM_L2]);
1599  MOFEM_LOG("EP", Sev::inform) << "norm_u_h1: " << std::sqrt(norms[U_NORM_H1]);
1600  MOFEM_LOG("EP", Sev::inform)
1601  << "norm_error_u_l2: " << std::sqrt(norms[U_ERROR_L2]);
1602  MOFEM_LOG("EP", Sev::inform)
1603  << "norm_piola: " << std::sqrt(norms[PIOLA_NORM]);
1604  CHKERR VecRestoreArrayRead(norms_vec, &norms);
1605 
1606  MoFEMFunctionReturn(0);
1607 }
1608 //! [Getting norms]
1609 
1610 MoFEMErrorCode EshelbianCore::getSpatialDispBc() {
1611  MoFEMFunctionBegin;
1612 
1613  auto bc_mng = mField.getInterface<BcManager>();
1614  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
1615  "", piolaStress, false, false);
1616 
1617  bcSpatialDispVecPtr = boost::make_shared<BcDispVec>();
1618 
1619  for (auto bc : bc_mng->getBcMapByBlockName()) {
1620  if (auto disp_bc = bc.second->dispBcPtr) {
1621 
1622  MOFEM_LOG("EP", Sev::noisy) << *disp_bc;
1623 
1624  std::vector<double> block_attributes(6, 0.);
1625  if (disp_bc->data.flag1 == 1) {
1626  block_attributes[0] = disp_bc->data.value1;
1627  block_attributes[3] = 1;
1628  }
1629  if (disp_bc->data.flag2 == 1) {
1630  block_attributes[1] = disp_bc->data.value2;
1631  block_attributes[4] = 1;
1632  }
1633  if (disp_bc->data.flag3 == 1) {
1634  block_attributes[2] = disp_bc->data.value3;
1635  block_attributes[5] = 1;
1636  }
1637  auto faces = bc.second->bcEnts.subset_by_dimension(2);
1638  bcSpatialDispVecPtr->emplace_back(bc.first, block_attributes, faces);
1639  }
1640  }
1641 
1642  // old way of naming blocksets for displacement BCs
1643  CHKERR getBc(bcSpatialDispVecPtr, "SPATIAL_DISP_BC", 6);
1644 
1645  MoFEMFunctionReturn(0);
1646 }
1647 
1648 MoFEMErrorCode EshelbianCore::getSpatialTractionBc() {
1649  MoFEMFunctionBegin;
1650 
1651  auto bc_mng = mField.getInterface<BcManager>();
1652  CHKERR bc_mng->pushMarkDOFsOnEntities<ForceCubitBcData>("", piolaStress,
1653  false, false);
1654 
1655  bcSpatialTraction = boost::make_shared<TractionBcVec>();
1656 
1657  for (auto bc : bc_mng->getBcMapByBlockName()) {
1658  if (auto force_bc = bc.second->forceBcPtr) {
1659  std::vector<double> block_attributes(6, 0.);
1660  block_attributes[0] = -force_bc->data.value3 * force_bc->data.value1;
1661  block_attributes[3] = 1;
1662  block_attributes[1] = -force_bc->data.value4 * force_bc->data.value1;
1663  block_attributes[4] = 1;
1664  block_attributes[2] = -force_bc->data.value5 * force_bc->data.value1;
1665  block_attributes[5] = 1;
1666  auto faces = bc.second->bcEnts.subset_by_dimension(2);
1667  bcSpatialTraction->emplace_back(bc.first, block_attributes, faces);
1668  }
1669  }
1670 
1671  CHKERR getBc(bcSpatialTraction, "SPATIAL_TRACTION_BC", 6);
1672 
1673  MoFEMFunctionReturn(0);
1674 }
1675 
1676 } // namespace EshelbianPlasticity
1677 
1678 #include <impl/EshelbianContact.cpp>
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Definition: UnknownInterface.hpp:93
CHK_MOAB_THROW
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
EshelbianPlasticity::EshelbianCore::f_linear
static double f_linear(const double v)
Definition: EshelbianPlasticity.hpp:897
EshelbianPlasticity::EshelbianCore::contactRhs
boost::shared_ptr< FaceElementForcesAndSourcesCore > contactRhs
Definition: EshelbianPlasticity.hpp:923
EshelbianPlasticity::OpSpatialPhysical_du_dP
Definition: EshelbianPlasticity.hpp:692
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
EshelbianPlasticity::BcRot::BcRot
BcRot(std::string name, std::vector< double > &attr, Range &faces)
Definition: EshelbianPlasticity.cpp:638
EshelbianPlasticity::OpMoveNode
Definition: EshelbianContact.hpp:193
EshelbianPlasticity::OpCalculateRotationAndSpatialGradient
Definition: EshelbianPlasticity.hpp:488
BoundaryEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Definition: child_and_parent.cpp:39
EshelbianPlasticity::EshelbianCore::essentialBcElement
const std::string essentialBcElement
Definition: EshelbianPlasticity.hpp:945
EshelbianPlasticity::EshelbianCore::addDMs
MoFEMErrorCode addDMs(const BitRefLevel bit=BitRefLevel().set(0))
Definition: EshelbianPlasticity.cpp:523
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:215
EshelbianPlasticity::LINEAR
@ LINEAR
Definition: EshelbianPlasticity.hpp:21
EshelbianPlasticity::CGGUserPolynomialBase
Definition: EshelbianPlasticity.cpp:752
EshelbianPlasticity::VolRule
Set integration rule on element.
Definition: EshelbianPlasticity.cpp:739
H1
@ H1
continuous field
Definition: definitions.h:85
EshelbianPlasticity::EshelbianCore::alphaW
double alphaW
Definition: EshelbianPlasticity.hpp:956
EshelbianPlasticity::EshelbianCore::bcSpatialRotationVecPtr
boost::shared_ptr< BcRotVec > bcSpatialRotationVecPtr
Definition: EshelbianPlasticity.hpp:966
EshelbianPlasticity::EshelbianCore::dd_f
static boost::function< double(const double)> dd_f
Definition: EshelbianPlasticity.hpp:871
EshelbianPlasticity::CGGUserPolynomialBase::CGGUserPolynomialBase
CGGUserPolynomialBase()
Definition: EshelbianPlasticity.cpp:754
TSElasticPostStep.cpp
sdf_hertz.d
float d
Definition: sdf_hertz.py:5
EshelbianPlasticity::EshelbianCore::inv_f_log
static double inv_f_log(const double v)
Definition: EshelbianPlasticity.hpp:887
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:666
EntityHandle
EshelbianPlasticity::EshelbianCore::piolaStress
const std::string piolaStress
Definition: EshelbianPlasticity.hpp:930
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
EshelbianPlasticity::DataAtIntegrationPts
Definition: EshelbianPlasticity.hpp:32
EshelbianPlasticity::EshelbianCore::spatialH1Disp
const std::string spatialH1Disp
Definition: EshelbianPlasticity.hpp:933
NOISY
@ NOISY
Definition: definitions.h:211
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:852
EshelbianPlasticity::EshelbianCore::exponentBase
static double exponentBase
Definition: EshelbianPlasticity.hpp:868
EshelbianPlasticity::EshelbianCore::inv_dd_f_linear
static double inv_dd_f_linear(const double v)
Definition: EshelbianPlasticity.hpp:903
MoFEM::OpFluxRhsImpl
Definition: Natural.hpp:39
MoFEM::OpCalculateTensor2SymmetricFieldValuesDot
Calculate symmetric tensor field rates ant integratio pts.
Definition: UserDataOperators.hpp:1072
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:916
EshelbianPlasticity::EshelbianCore::spaceH1Order
int spaceH1Order
Definition: EshelbianPlasticity.hpp:953
MoFEM::PostProcBrokenMeshInMoabBase::getMapGaussPts
auto & getMapGaussPts()
Get vector of vectors associated to integration points.
Definition: PostProcBrokenMeshInMoabBase.hpp:637
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:460
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
TSElasticPostStep::postStepFun
static MoFEMErrorCode postStepFun(TS ts)
Definition: TSElasticPostStep.cpp:135
MoFEM::OpCalculateHTensorTensorField
Calculate tenor field using tensor base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:2635
EshelbianPlasticity::EshelbianCore::alphaRho
double alphaRho
Definition: EshelbianPlasticity.hpp:957
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:29
EshelbianPlasticity::EshelbianCore::addBoundaryFiniteElement
MoFEMErrorCode addBoundaryFiniteElement(const EntityHandle meshset=0)
Definition: EshelbianPlasticity.cpp:436
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:899
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
EshelbianPlasticity::OpRotationBc
Definition: EshelbianPlasticity.hpp:589
EshelbianPlasticity::EshelbianCore::stretchTensor
const std::string stretchTensor
Definition: EshelbianPlasticity.hpp:936
EshelbianPlasticity::FaceUserDataOperatorStabAssembly
Definition: EshelbianPlasticity.cpp:39
EshelbianPlasticity::OpSpatialConsistencyP
Definition: EshelbianPlasticity.hpp:538
HenckyOps::d_f
auto d_f
Definition: HenckyOps.hpp:16
MoFEM::EntPolynomialBaseCtx::dAta
EntitiesFieldData & dAta
Definition: EntPolynomialBaseCtx.hpp:36
EshelbianPlasticity::EshelbianCore::setElasticElementToTs
MoFEMErrorCode setElasticElementToTs(DM dm)
Definition: EshelbianPlasticity.cpp:1265
EshelbianPlasticity::EshelbianCore::setElasticElementOps
MoFEMErrorCode setElasticElementOps(const int tag)
Definition: EshelbianPlasticity.cpp:1239
EshelbianMonitor.cpp
Contains definition of EshelbianMonitor class.
EshelbianPlasticity
Definition: CGGTonsorialBubbleBase.hpp:11
EshelbianPlasticity::EshelbianCore::d_f_linear
static double d_f_linear(const double v)
Definition: EshelbianPlasticity.hpp:898
get_range_from_block
static auto get_range_from_block(MoFEM::Interface &m_field, const std::string block_name)
Definition: EshelbianPlasticity.cpp:46
MoFEM.hpp
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
EshelbianPlasticity::OpSpatialConsistencyDivTerm
Definition: EshelbianPlasticity.hpp:556
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:882
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:527
MoFEM::DisplacementCubitBcData
Definition of the displacement bc data structure.
Definition: BCData.hpp:76
MoFEM::BaseFunction
Base class if inherited used to calculate base functions.
Definition: BaseFunction.hpp:40
SideEle
PostProcEleByDim< SPACE_DIM >::SideEle SideEle
Definition: adolc_plasticity.cpp:98
EshelbianPlasticity::EshelbianCore::inv_f_linear
static double inv_f_linear(const double v)
Definition: EshelbianPlasticity.hpp:901
MoFEM::PostProcBrokenMeshInMoabBase::writeFile
MoFEMErrorCode writeFile(const std::string file_name)
wrote results in (MOAB) format, use "file_name.h5m"
Definition: PostProcBrokenMeshInMoabBase.hpp:655
EshelbianPlasticity::EshelbianCore::d_f
static boost::function< double(const double)> d_f
Definition: EshelbianPlasticity.hpp:870
EshelbianPlasticity::EshelbianCore::dmPrjSpatial
SmartPetscObj< DM > dmPrjSpatial
Projection spatial displacement.
Definition: EshelbianPlasticity.hpp:928
EshelbianPlasticity::EshelbianCore::elasticFeRhs
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeRhs
Definition: EshelbianPlasticity.hpp:919
ts_ctx
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
EshelbianPlasticity::EshelbianCore::elasticBcRhs
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcRhs
Definition: EshelbianPlasticity.hpp:922
EshelbianPlasticity::EshelbianCore::elasticFeLhs
boost::shared_ptr< VolumeElementForcesAndSourcesCore > elasticFeLhs
Definition: EshelbianPlasticity.hpp:920
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.cpp:794
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:501
order
constexpr int order
Definition: dg_projection.cpp:18
EshelbianPlasticity::EshelbianCore::gradApperoximator
static enum RotSelector gradApperoximator
Definition: EshelbianPlasticity.hpp:865
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
EshelbianPlasticity::EshelbianCore::f_log
static double f_log(const double v)
Definition: EshelbianPlasticity.hpp:876
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::BcDisp::flags
VectorInt3 flags
Definition: EshelbianPlasticity.hpp:301
MoFEM::DMMoFEMSetDestroyProblem
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMoFEM.cpp:442
EshelbianPlasticity::BcDisp::faces
Range faces
Definition: EshelbianPlasticity.hpp:299
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:69
EshelbianPlasticity::EshelbianCore::postProcessResults
MoFEMErrorCode postProcessResults(const int tag, const std::string file)
Definition: EshelbianPlasticity.cpp:1465
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:170
EshelbianPlasticity::OpSpatialRotation
Definition: EshelbianPlasticity.hpp:508
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
EshelbianPlasticity::TractionBc::faces
Range faces
Definition: EshelbianPlasticity.hpp:319
EshelbianPlasticity::OpSpatialConsistency_dBubble_domega
Definition: EshelbianPlasticity.hpp:790
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:857
EshelbianPlasticity::EshelbianCore::precEpsW
double precEpsW
Definition: EshelbianPlasticity.hpp:960
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:87
EshelbianPlasticity::EshelbianCore::skinElement
const std::string skinElement
Definition: EshelbianPlasticity.hpp:946
EshelbianPlasticity::EshelbianCore::d_f_log
static double d_f_log(const double v)
Definition: EshelbianPlasticity.hpp:879
EshelbianPlasticity::EshelbianCore::spaceOrder
int spaceOrder
Definition: EshelbianPlasticity.hpp:952
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:137
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
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:1021
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:396
EshelbianPlasticity::BcRot::vals
VectorDouble3 vals
Definition: EshelbianPlasticity.hpp:309
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
EshelbianPlasticity::OpSpatialConsistency_dP_domega
Definition: EshelbianPlasticity.hpp:776
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:740
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:979
EshelbianPlasticity::EshelbianCore::stretchSelector
static enum StretchSelector stretchSelector
Definition: EshelbianPlasticity.hpp:866
EshelbianPlasticity::EshelbianCore::dmElastic
SmartPetscObj< DM > dmElastic
Elastic problem.
Definition: EshelbianPlasticity.hpp:927
EshelbianPlasticity::EshelbianCore::elasticBcLhs
boost::shared_ptr< FaceElementForcesAndSourcesCore > elasticBcLhs
Definition: EshelbianPlasticity.hpp:921
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::contactTreeRhs
boost::shared_ptr< ContactTree > contactTreeRhs
Make a contact tree.
Definition: EshelbianPlasticity.hpp:924
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:219
double
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
MoFEM::ProblemsManager::removeDofsOnEntities
MoFEMErrorCode removeDofsOnEntities(const std::string problem_name, const std::string field_name, const Range ents, const int lo_coeff=0, const int hi_coeff=MAX_DOFS_ON_ENTITY, const int lo_order=0, const int hi_order=100, int verb=VERBOSE, const bool debug=false)
Remove DOFs from problem.
Definition: ProblemsManager.cpp:2605
convert.type
type
Definition: convert.py:64
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:302
EshelbianPlasticity::EshelbianCore::precEpsContactDisp
double precEpsContactDisp
Definition: EshelbianPlasticity.hpp:961
DM_NO_ELEMENT
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:10
MoFEM::PairNameFEMethodPtr
Definition: AuxPETSc.hpp:12
EshelbianPlasticity::EshelbianCore::rotSelector
static enum RotSelector rotSelector
Definition: EshelbianPlasticity.hpp:864
EshelbianPlasticity::OpConstrainBoundaryHDivLhs_dU
Definition: EshelbianContact.hpp:177
EshelbianPlasticity::EshelbianCore::f
static boost::function< double(const double)> f
Definition: EshelbianPlasticity.hpp:869
EshelbianPlasticity::EshelbianCore::bcSpatialDispVecPtr
boost::shared_ptr< BcDispVec > bcSpatialDispVecPtr
Definition: EshelbianPlasticity.hpp:965
EshelbianPlasticity::OpConstrainBoundaryHDivRhs
Definition: EshelbianContact.hpp:129
EshelbianPlasticity::VolUserDataOperatorStabAssembly
Definition: EshelbianPlasticity.cpp:35
EshelbianPlasticity::CGGUserPolynomialBase::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, BaseFunctionUnknownInterface **iface) const
Definition: EshelbianPlasticity.cpp:757
EshelbianPlasticity::EshelbianCore::inv_d_f
static boost::function< double(const double)> inv_d_f
Definition: EshelbianPlasticity.hpp:873
MoFEM::EntPolynomialBaseCtx::sPace
const FieldSpace sPace
Definition: EntPolynomialBaseCtx.hpp:37
MoFEM::OpCalculateHVecTensorField
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Definition: UserDataOperators.hpp:2576
EshelbianPlasticity::StretchSelector
StretchSelector
Definition: EshelbianPlasticity.hpp:21
EshelbianPlasticity::EshelbianCore::bcSpatialTraction
boost::shared_ptr< TractionBcVec > bcSpatialTraction
Definition: EshelbianPlasticity.hpp:967
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:128
EshelbianPlasticity::EshelbianCore::mField
MoFEM::Interface & mField
Definition: EshelbianPlasticity.hpp:914
EshelbianPlasticity::OpDispBc
Definition: EshelbianPlasticity.hpp:563
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dU
Definition: EshelbianContact.hpp:143
EshelbianPlasticity::FeTractionBc
Definition: EshelbianPlasticity.hpp:624
EshelbianPlasticity.hpp
Eshelbian plasticity interface.
EshelbianPlasticity::BcDisp::BcDisp
BcDisp(std::string name, std::vector< double > &attr, Range &faces)
Definition: EshelbianPlasticity.cpp:621
EshelbianPlasticity::CGGUserPolynomialBase::getValueHdivForCGGBubble
MoFEMErrorCode getValueHdivForCGGBubble(MatrixDouble &pts)
Definition: EshelbianPlasticity.cpp:797
MoFEM::OpCalcNormL2Tensor1
Get norm of input MatrixDouble for Tensor1.
Definition: NormsOperators.hpp:42
HenckyOps::dd_f
auto dd_f
Definition: HenckyOps.hpp:17
EshelbianPlasticity::FaceRule
Definition: EshelbianPlasticity.cpp:748
EshelbianPlasticity::EshelbianCore::contactElement
const std::string contactElement
Definition: EshelbianPlasticity.hpp:947
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:1023
EshelbianPlasticity::EshelbianCore::dM
SmartPetscObj< DM > dM
Coupled problem all fields.
Definition: EshelbianPlasticity.hpp:926
EshelbianPlasticity::CGGUserPolynomialBase::~CGGUserPolynomialBase
~CGGUserPolynomialBase()
Definition: EshelbianPlasticity.cpp:755
MoFEM::matDuplicate
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
Definition: PetscSmartObj.hpp:230
EshelbianPlasticity::OpSpatialPhysical_du_domega
Definition: EshelbianPlasticity.hpp:718
EshelbianPlasticity::EshelbianCore
enum RotSelector EshelbianCore
Definition: EshelbianPlasticity.cpp:84
MoFEM::BaseFunctionUnknownInterface
Definition: BaseFunction.hpp:13
EshelbianPlasticity::EshelbianCore::rotAxis
const std::string rotAxis
Definition: EshelbianPlasticity.hpp:937
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:118
EshelbianPlasticity::CGGUserPolynomialBase::cTx
EntPolynomialBaseCtx * cTx
Definition: EshelbianPlasticity.cpp:792
EshelbianPlasticity::OpPostProcDataStructure
Definition: EshelbianPlasticity.hpp:804
EshelbianPlasticity::EshelbianCore::spatialL2Disp
const std::string spatialL2Disp
Definition: EshelbianPlasticity.hpp:932
BiLinearForm
EshelbianPlasticity::OpSpatialPhysical_du_du
Definition: EshelbianPlasticity.hpp:667
MoFEM::EntPolynomialBaseCtx::bAse
const FieldApproximationBase bAse
Definition: EntPolynomialBaseCtx.hpp:38
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
field_name
constexpr auto field_name
Definition: poisson_2d_homogeneous.cpp:13
EshelbianPlasticity::EshelbianCore::bcSpatialFreeTraction
boost::shared_ptr< TractionFreeBc > bcSpatialFreeTraction
Definition: EshelbianPlasticity.hpp:968
MoFEM::OpCalculateHVecTensorTrace
Calculate trace of vector (Hdiv/Hcurl) space.
Definition: UserDataOperators.hpp:2764
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:503
TSElasticPostStep::postStepDestroy
static MoFEMErrorCode postStepDestroy()
Definition: TSElasticPostStep.cpp:77
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1536
EshelbianPlasticity::EshelbianCore::bubbleField
const std::string bubbleField
Definition: EshelbianPlasticity.hpp:941
MoFEM::UnknownInterface
base class for all interface classes
Definition: UnknownInterface.hpp:34
cholesky.hpp
cholesky decomposition
MAX_DOFS_ON_ENTITY
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
MoFEM::ForcesAndSourcesCore
structure to get information form mofem into EntitiesFieldData
Definition: ForcesAndSourcesCore.hpp:22
EshelbianPlasticity::EshelbianCore
Definition: EshelbianPlasticity.hpp:862
MoFEM::OpCalculateHVecTensorDivergence
Calculate divergence of tonsorial field using vectorial base.
Definition: UserDataOperators.hpp:2696
EshelbianPlasticity::EshelbianCore::inv_f
static boost::function< double(const double)> inv_f
Definition: EshelbianPlasticity.hpp:872
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
Range
EshelbianPlasticity::CGGUserPolynomialBase::getValue
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Definition: EshelbianPlasticity.cpp:763
EshelbianPlasticity::EshelbianCore::SetUpSchur::createSetUpSchur
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, SmartPetscObj< Mat > m, SmartPetscObj< Mat > p, EshelbianCore *ep_core_ptr)
Definition: SetUpSchurImpl.cpp:230
EshelbianPlasticity::SMALL_ROT
@ SMALL_ROT
Definition: EshelbianPlasticity.hpp:20
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:98
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
CGGTonsorialBubbleBase.hpp
Implementation of tonsorial bubble base div(v) = 0.
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
EshelbianPlasticity::EshelbianCore::physicalEquations
boost::shared_ptr< PhysicalEquations > physicalEquations
Definition: EshelbianPlasticity.hpp:917
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
EshelbianPlasticity::OpConstrainBoundaryL2Rhs
Definition: EshelbianContact.hpp:114
EshelbianPlasticity::LARGE_ROT
@ LARGE_ROT
Definition: EshelbianPlasticity.hpp:20
EshelbianPlasticity::LOG
@ LOG
Definition: EshelbianPlasticity.hpp:21
EshelbianPlasticity::EshelbianCore::elementVolumeName
const std::string elementVolumeName
Definition: EshelbianPlasticity.hpp:943
NBVOLUMETET_CCG_BUBBLE
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
Definition: CGGTonsorialBubbleBase.hpp:19
EshelbianPlasticity::EshelbianCore::precEps
double precEps
Definition: EshelbianPlasticity.hpp:958
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
EshelbianPlasticity::TractionBc::TractionBc
TractionBc(std::string name, std::vector< double > &attr, Range &faces)
Definition: EshelbianPlasticity.cpp:647
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
EshelbianPlasticity::RotSelector
RotSelector
Definition: EshelbianPlasticity.hpp:20
EshelbianPlasticity::TractionBc::flags
VectorInt3 flags
Definition: EshelbianPlasticity.hpp:321
EshelbianPlasticity::EshelbianCore::inv_d_f_linear
static double inv_d_f_linear(const double v)
Definition: EshelbianPlasticity.hpp:902
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:804
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
MoFEM::OpCalculateTensor2SymmetricFieldValues
Calculate symmetric tensor field values at integration pts.
Definition: UserDataOperators.hpp:978
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
EshelbianPlasticity::EshelbianCore::alphaU
double alphaU
Definition: EshelbianPlasticity.hpp:955
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:198
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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_filed)=0
set finite element field data
EshelbianPlasticity::EshelbianCore::inv_d_f_log
static double inv_d_f_log(const double v)
Definition: EshelbianPlasticity.hpp:890
MoFEM::DMMoFEMGetTsCtx
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMoFEM.cpp:1146
EshelbianPlasticity::EshelbianCore::addFields
MoFEMErrorCode addFields(const EntityHandle meshset=0)
Definition: EshelbianPlasticity.cpp:256
TSElasticPostStep::postStepInitialise
static MoFEMErrorCode postStepInitialise(EshelbianCore *ep_ptr)
Definition: TSElasticPostStep.cpp:11
MoFEM::EntitiesFieldData::dataOnEntities
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
Definition: EntitiesFieldData.hpp:56
EshelbianPlasticity::EshelbianCore::inv_dd_f
static boost::function< double(const double)> inv_dd_f
Definition: EshelbianPlasticity.hpp:874
EshelbianPlasticity::EshelbianCore::inv_dd_f_log
static double inv_dd_f_log(const double v)
Definition: EshelbianPlasticity.hpp:893
EshelbianPlasticity::FaceRule::operator()
int operator()(int p_row, int p_col, int p_data) const
Definition: EshelbianPlasticity.cpp:749
MoFEM::PostProcBrokenMeshInMoabBase::getPostProcMesh
auto & getPostProcMesh()
Get postprocessing mesh.
Definition: PostProcBrokenMeshInMoabBase.hpp:641
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
EshelbianPlasticity::EshelbianCore::getOptions
MoFEMErrorCode getOptions()
Definition: EshelbianPlasticity.cpp:135
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::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:440
EshelbianPlasticity::BcRot::theta
double theta
Definition: EshelbianPlasticity.hpp:310
EshelbianPlasticity::EshelbianCore::setContactElementOps
MoFEMErrorCode setContactElementOps(boost::shared_ptr< ContactTree > &fe_contact_tree, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_rhs, boost::shared_ptr< FaceElementForcesAndSourcesCore > &fe_lhs)
Definition: EshelbianPlasticity.cpp:1057
QUIET
@ QUIET
Definition: definitions.h:208
EshelbianPlasticity::EshelbianCore::getBc
MoFEMErrorCode getBc(boost::shared_ptr< BC > &bc_vec_ptr, const std::string block_name, const int nb_attributes)
Definition: EshelbianPlasticity.hpp:971
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:242
EshelbianPlasticity::OpSpatialEquilibrium
Definition: EshelbianPlasticity.hpp:497
SetUpSchurImpl.cpp
EshelbianPlasticity::OpSpatialEquilibrium_dw_dP
Definition: EshelbianPlasticity.hpp:641
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:575
MoFEM::EntitiesFieldData
data structure for finite element entity
Definition: EntitiesFieldData.hpp:40
EshelbianPlasticity::BcDisp::vals
VectorDouble3 vals
Definition: EshelbianPlasticity.hpp:300
MoFEM::SmartPetscObj< Mat >
EshelbianPlasticity::EshelbianCore::getSpatialTractionBc
MoFEMErrorCode getSpatialTractionBc()
Definition: EshelbianPlasticity.cpp:1648
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:918
EshelbianPlasticity::EshelbianCore::materialOrder
int materialOrder
Definition: EshelbianPlasticity.hpp:954
EshelbianPlasticity::OpConstrainBoundaryL2Lhs_dP
Definition: EshelbianContact.hpp:160
EshelbianPlasticity::OpSpatialRotation_domega_domega
Definition: EshelbianPlasticity.hpp:766
EshelbianPlasticity::MODERATE_ROT
@ MODERATE_ROT
Definition: EshelbianPlasticity.hpp:20
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:590
EshelbianPlasticity::OpSpatialEquilibrium_dw_dw
Definition: EshelbianPlasticity.hpp:653
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
EshelbianPlasticity::EshelbianCore::solveElastic
MoFEMErrorCode solveElastic(TS ts, Mat m, Vec f, Vec x)
Definition: EshelbianPlasticity.cpp:1309
HDIV
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
MOFEM_NOT_IMPLEMENTED
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
EshelbianPlasticity::CGGUserPolynomialBase::cachePhi
boost::tuple< int, int, MatrixDouble > cachePhi
Definition: EshelbianPlasticity.cpp:795
MoFEM::ForcesAndSourcesCore::UserDataOperator::getKSPB
Mat getKSPB() const
Definition: ForcesAndSourcesCore.hpp:1103
EshelbianPlasticity::OpSpatialRotation_domega_dBubble
Definition: EshelbianPlasticity.hpp:754
MoFEM::getDMTsCtx
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1060
CHKERRG
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:483
MoFEM::SCHUR
@ SCHUR
Definition: FormsIntegrators.hpp:104
EshelbianPlasticity::OpSpatialPhysical
Definition: EshelbianPlasticity.hpp:515
EshelbianPlasticity::OpSpatialRotation_domega_dP
Definition: EshelbianPlasticity.hpp:742
EshelbianPlasticity::EshelbianCore::contactDisp
const std::string contactDisp
Definition: EshelbianPlasticity.hpp:934
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1289
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
EshelbianPlasticity::EshelbianCore::getSpatialDispBc
MoFEMErrorCode getSpatialDispBc()
[Getting norms]
Definition: EshelbianPlasticity.cpp:1610
EshelbianPlasticity::OpSpatialConsistencyBubble
Definition: EshelbianPlasticity.hpp:547
EshelbianPlasticity::EshelbianCore::naturalBcElement
const std::string naturalBcElement
Definition: EshelbianPlasticity.hpp:944
EshelbianPlasticity::EshelbianCore::precEpsOmega
double precEpsOmega
Definition: EshelbianPlasticity.hpp:959
EshelbianPlasticity::EshelbianCore::gettingNorms
MoFEMErrorCode gettingNorms()
[Getting norms]
Definition: EshelbianPlasticity.cpp:1549
EshelbianPlasticity::EshelbianCore::~EshelbianCore
virtual ~EshelbianCore()
EshelbianPlasticity::TractionBc::vals
VectorDouble3 vals
Definition: EshelbianPlasticity.hpp:320
MoFEM::PostProcBrokenMeshInMoab< FaceElementForcesAndSourcesCore >
Definition: PostProcBrokenMeshInMoabBase.hpp:677
MoFEM::ForcesAndSourcesCore::UserDataOperator::UserDataOperator
UserDataOperator(const FieldSpace space, const char type=OPSPACE, const bool symm=true)
Definition: ForcesAndSourcesCore.cpp:1942
EshelbianPlasticity::OpSpatialPhysical_du_dBubble
Definition: EshelbianPlasticity.hpp:705
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1127
EshelbianPlasticity::OpTreeSearch
Definition: EshelbianContact.hpp:209
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
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
Generated by Doxygen 1.8.17 and hosted at University of Glasgow