v0.9.2
EshelbianPlasticity.cpp
Go to the documentation of this file.
1 /**
2  * \file EshelbianPlasticity.cpp
3  * \example EshelbianPlasticity.cpp
4  *
5  * \brief Eshelbian plasticity implementation
6  */
7 
8 #include <MoFEM.hpp>
9 using namespace MoFEM;
10 
11 #include <BasicFiniteElements.hpp>
13 
14 #include <EshelbianPlasticity.hpp>
15 #include <boost/math/constants/constants.hpp>
16 
17 namespace EshelbianPlasticity {
18 
19 MoFEMErrorCode EshelbianCore::query_interface(const MOFEMuuid &uuid,
20  UnknownInterface **iface) const {
22  *iface = NULL;
24  *iface = const_cast<EshelbianCore *>(this);
26  }
27  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
29 }
30 
31 MoFEMErrorCode OpJacobian::doWork(int side, EntityType type, EntData &data) {
33  if (type != MBVERTEX)
35 
36  if (evalRhs)
37  CHKERR evaluateRhs(data);
38 
39  if (evalLhs)
40  CHKERR evaluateLhs(data);
41 
43 }
44 
45 EshelbianCore::EshelbianCore(MoFEM::Interface &m_field)
46  : mField(m_field), piolaStress("P"), eshelbyStress("S"),
47  spatialDisp("w"), materialDisp("W"),
48  streachTensor("u"), rotAxis("omega"), materialGradient("G"),
49  tauField("TAU"), lambdaField("LAMBDA"), bubbleField("BUBBLE"),
50  elementVolumeName("EP"), naturalBcElement("NATURAL_BC"),
51  essentialBcElement("ESSENTIAL_BC") {
52 
53  ierr = getOptions();
54  CHKERRABORT(PETSC_COMM_WORLD, ierr);
55 }
56 
58 }
59 
62  CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Eshelbian plasticity",
63  "none");
64 
65  spaceOrder = 1;
66  CHKERR PetscOptionsInt("-space_order", "approximation oder for space", "",
67  spaceOrder, &spaceOrder, PETSC_NULL);
68  materialOrder = 1;
69  CHKERR PetscOptionsInt("-material_order", "approximation oder for material",
70  "", materialOrder, &materialOrder, PETSC_NULL);
71 
72  alpha_u = 0;
73  CHKERR PetscOptionsScalar("-viscosity_alpha_u", "viscosity", "", alpha_u,
74  &alpha_u, PETSC_NULL);
75 
76  alpha_w = 0;
77  CHKERR PetscOptionsScalar("-viscosity_alpha_w", "viscosity", "", alpha_w,
78  &alpha_w, PETSC_NULL);
79 
80  preconditioner_eps = 1e-3;
81  CHKERR PetscOptionsScalar("-preconditioner_eps", "preconditioner_eps", "",
83  PETSC_NULL);
84 
85  ierr = PetscOptionsEnd();
86  CHKERRG(ierr);
88 }
89 
92 
93  Range tets;
94  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
95  Range tets_skin_part;
96  Skinner skin(&mField.get_moab());
97  CHKERR skin.find_skin(0, tets, false, tets_skin_part);
98  ParallelComm *pcomm =
99  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
100  Range tets_skin;
101  CHKERR pcomm->filter_pstatus(tets_skin_part,
102  PSTATUS_SHARED | PSTATUS_MULTISHARED,
103  PSTATUS_NOT, -1, &tets_skin);
104 
105  auto subtract_faces_where_displacements_are_applied =
106  [&](const std::string disp_block_set_name) {
109  if (it->getName().compare(0, disp_block_set_name.length(),
110  disp_block_set_name) == 0) {
111  Range faces;
112  CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2,
113  faces, true);
114  tets_skin = subtract(tets_skin, faces);
115  }
116  }
118  };
119  CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_DISP_BC");
120  CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_ROTATION_BC");
121  CHKERR subtract_faces_where_displacements_are_applied("SPATIAL_TRACTION_BC");
122 
123  Range faces;
124  CHKERR mField.get_moab().get_adjacencies(tets, 2, true, faces,
125  moab::Interface::UNION);
126  Range faces_not_on_the_skin = subtract(faces, tets_skin);
127 
128  auto add_hdiv_field = [&](const std::string field_name, const int order,
129  const int dim) {
132  MB_TAG_SPARSE, MF_ZERO);
133  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
134  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
135  CHKERR mField.set_field_order(faces_not_on_the_skin, field_name, order);
136  CHKERR mField.set_field_order(tets_skin, field_name, 0);
138  };
139 
140  auto add_hdiv_rt_field = [&](const std::string field_name, const int order,
141  const int dim) {
144  MB_TAG_DENSE, MF_ZERO);
145  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
146  CHKERR mField.set_field_order(meshset, MBTET, field_name, 0);
147  CHKERR mField.set_field_order(tets_skin, field_name, order);
149  };
150 
151  auto add_l2_field = [this, meshset](const std::string field_name,
152  const int order, const int dim) {
155  MB_TAG_DENSE, MF_ZERO);
156  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
157  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
159  };
160 
161  auto add_h1_field = [this, meshset](const std::string field_name,
162  const int order, const int dim) {
165  MB_TAG_DENSE, MF_ZERO);
166  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
167  CHKERR mField.set_field_order(meshset, MBVERTEX, field_name, 1);
168  CHKERR mField.set_field_order(meshset, MBEDGE, field_name, order);
169  CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
170  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
172  };
173 
174  auto add_bubble_field = [this, meshset](const std::string field_name,
175  const int order, const int dim) {
177  CHKERR mField.add_field(field_name, HDIV, USER_BASE, dim, MB_TAG_DENSE,
178  MF_ZERO);
179  // Modify field
180  auto field_ptr = mField.get_field_structure(field_name);
181  auto field_order_table =
182  const_cast<Field *>(field_ptr)->getFieldOrderTable();
183  auto get_cgg_bubble_order_zero = [](int p) { return 0; };
184  auto get_cgg_bubble_order_tet = [](int p) {
185  return NBVOLUMETET_CCG_BUBBLE(p);
186  };
187  field_order_table[MBVERTEX] = get_cgg_bubble_order_zero;
188  field_order_table[MBEDGE] = get_cgg_bubble_order_zero;
189  field_order_table[MBTRI] = get_cgg_bubble_order_zero;
190  field_order_table[MBTET] = get_cgg_bubble_order_tet;
191  CHKERR mField.add_ents_to_field_by_type(meshset, MBTET, field_name);
192  CHKERR mField.set_field_order(meshset, MBTRI, field_name, order);
193  CHKERR mField.set_field_order(meshset, MBTET, field_name, order);
195  };
196 
197  // spatial fields
198  CHKERR add_hdiv_field(piolaStress, spaceOrder, 3);
199  CHKERR add_bubble_field(bubbleField, spaceOrder, 1);
200  CHKERR add_l2_field(spatialDisp, spaceOrder - 1, 3);
201  CHKERR add_l2_field(rotAxis, spaceOrder - 1, 3);
202  CHKERR add_l2_field(streachTensor, spaceOrder, 6);
203 
204  // material fields
205  CHKERR add_hdiv_field(eshelbyStress, materialOrder, 3);
206  CHKERR add_l2_field(materialGradient, materialOrder - 1, 9);
207  // CHKERR add_l2_field(materialDisp, materialOrder - 1, 3);
208  // CHKERR add_l2_field(tauField, materialOrder - 1, 1);
209  // CHKERR add_l2_field(lambdaField, materialOrder - 1, 1);
210 
211  // Add history filedes
212  CHKERR add_l2_field(materialGradient + "0", materialOrder - 1, 9);
213 
215 
217 }
218 
222 
223  // set finite element fields
224  auto add_field_to_fe = [this](const std::string fe,
225  const std::string field_name) {
231  };
232 
237 
238  CHKERR add_field_to_fe(elementVolumeName, piolaStress);
239  CHKERR add_field_to_fe(elementVolumeName, bubbleField);
240  CHKERR add_field_to_fe(elementVolumeName, eshelbyStress);
241  CHKERR add_field_to_fe(elementVolumeName, streachTensor);
242  CHKERR add_field_to_fe(elementVolumeName, rotAxis);
243  CHKERR add_field_to_fe(elementVolumeName, spatialDisp);
244  CHKERR add_field_to_fe(elementVolumeName, streachTensor);
245  CHKERR add_field_to_fe(elementVolumeName, materialGradient);
246 
248  materialGradient + "0");
249  }
250 
251  // build finite elements data structures
253 
255 }
256 
260 
261  auto bc_elements_add_to_range = [&](const std::string disp_block_set_name,
262  Range &r) {
265  if (it->getName().compare(0, disp_block_set_name.length(),
266  disp_block_set_name) == 0) {
267  Range faces;
268  CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
269  true);
270  r.merge(faces);
271  }
272  }
274  };
275 
276  // set finite element fields
277  auto add_field_to_fe = [this](const std::string fe,
278  const std::string field_name) {
284  };
285 
286  Range natural_bc_elements;
287  CHKERR bc_elements_add_to_range("SPATIAL_DISP_BC", natural_bc_elements);
288  CHKERR bc_elements_add_to_range("SPATIAL_ROTATION_BC", natural_bc_elements);
289  Range essentail_bc_elements;
290  CHKERR bc_elements_add_to_range("SPATIAL_TRACTION_BC", essentail_bc_elements);
291 
293  CHKERR mField.add_ents_to_finite_element_by_type(natural_bc_elements, MBTRI,
295  CHKERR add_field_to_fe(naturalBcElement, piolaStress);
296  CHKERR add_field_to_fe(naturalBcElement, eshelbyStress);
298 
300  CHKERR mField.add_ents_to_finite_element_by_type(essentail_bc_elements, MBTRI,
302  CHKERR add_field_to_fe(essentialBcElement, piolaStress);
303  CHKERR add_field_to_fe(essentialBcElement, eshelbyStress);
305 
307 }
308 
311 
312  // find adjacencies between finite elements and dofs
314 
315  // Create coupled problem
316  dM = createSmartDM(mField.get_comm(), "DMMOFEM");
317  CHKERR DMMoFEMCreateMoFEM(dM, &mField, "ESHELBY_PLASTICITY", bit,
318  BitRefLevel().set());
319  CHKERR DMMoFEMSetDestroyProblem(dM, PETSC_TRUE);
320  CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
324  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
325  CHKERR DMSetUp(dM);
326  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
327 
328  auto remove_dofs_on_essential_spatial_stress_boundary =
329  [&](const std::string prb_name) {
331  for (int d : {0, 1, 2})
333  prb_name, piolaStress, (*bcSpatialFreeTraction)[d], d, d, NOISY,
334  true);
336  };
337  CHKERR remove_dofs_on_essential_spatial_stress_boundary("ESHELBY_PLASTICITY");
338 
339  // Create elastic sub-problem
340  dmElastic = createSmartDM(mField.get_comm(), "DMMOFEM");
341  CHKERR DMMoFEMCreateSubDM(dmElastic, dM, "ELASTIC_PROBLEM");
353  CHKERR DMSetUp(dmElastic);
354 
355  // Create elastic streach-problem
358  "ELASTIC_PROBLEM_STREACH_SCHUR");
369  CHKERR DMSetUp(dmElasticSchurStreach);
370 
371  // Create elastic bubble-problem
374  "ELASTIC_PROBLEM_BUBBLE_SCHUR");
384  CHKERR DMSetUp(dmElasticSchurBubble);
385 
386  // Create elastic omega-problem
389  "ELASTIC_PROBLEM_OMEGA_SCHUR");
398  CHKERR DMSetUp(dmElasticSchurOmega);
399 
400  // Create elastic tet_stress-problem
403  "ELASTIC_PROBLEM_SPATIAL_DISP_SCHUR");
407  elementVolumeName.c_str());
410  essentialBcElement.c_str());
414 
415  {
416  PetscSection section;
417  CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
418  &section);
419  CHKERR DMSetDefaultSection(dmElastic, section);
420  CHKERR DMSetDefaultGlobalSection(dmElastic, section);
421  CHKERR PetscSectionDestroy(&section);
422  }
423 
425 }
426 
427 BcDisp::BcDisp(std::string name, std::vector<double> &attr, Range &faces)
428  : blockName(name), faces(faces) {
429  vals.resize(3, false);
430  flags.resize(3, false);
431  for (int ii = 0; ii != 3; ++ii) {
432  vals[ii] = attr[ii];
433  flags[ii] = static_cast<int>(attr[ii + 3]);
434  }
435 }
436 
437 BcRot::BcRot(std::string name, std::vector<double> &attr, Range &faces)
438  : blockName(name), faces(faces) {
439  vals.resize(3, false);
440  for (int ii = 0; ii != 3; ++ii) {
441  vals[ii] = attr[ii];
442  }
443  theta = attr[3];
444 }
445 
446 TractionBc::TractionBc(std::string name, std::vector<double> &attr,
447  Range &faces)
448  : blockName(name), faces(faces) {
449  vals.resize(3, false);
450  flags.resize(3, false);
451  for (int ii = 0; ii != 3; ++ii) {
452  vals[ii] = attr[ii];
453  flags[ii] = static_cast<int>(attr[ii + 3]);
454  }
455 }
456 
459  boost::shared_ptr<TractionFreeBc> &bc_ptr,
460  const std::string disp_block_set_name,
461  const std::string rot_block_set_name) {
463  Range tets;
464  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
465  Range tets_skin_part;
466  Skinner skin(&mField.get_moab());
467  CHKERR skin.find_skin(0, tets, false, tets_skin_part);
468  ParallelComm *pcomm =
469  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
470  Range tets_skin;
471  CHKERR pcomm->filter_pstatus(tets_skin_part,
472  PSTATUS_SHARED | PSTATUS_MULTISHARED,
473  PSTATUS_NOT, -1, &tets_skin);
474 
475  bc_ptr->resize(3);
476  for (int dd = 0; dd != 3; ++dd)
477  (*bc_ptr)[dd] = tets_skin;
478 
480  if (it->getName().compare(0, disp_block_set_name.length(),
481  disp_block_set_name) == 0) {
482  std::vector<double> block_attributes;
483  CHKERR it->getAttributes(block_attributes);
484  if (block_attributes.size() != 6) {
485  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
486  "In block %s six attributes are required for given BC "
487  "blockset (3 values + "
488  "3 flags) != %d",
489  it->getName().c_str(), block_attributes.size());
490  }
491  Range faces;
492  CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
493  true);
494  if (block_attributes[3] != 0)
495  (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
496  if (block_attributes[4] != 0)
497  (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
498  if (block_attributes[5] != 0)
499  (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
500  }
501  if (it->getName().compare(0, rot_block_set_name.length(),
502  rot_block_set_name) == 0) {
503  Range faces;
504  CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
505  true);
506  (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
507  (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
508  (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
509  }
510  }
511 
512  // for (int dd = 0; dd != 3; ++dd) {
513  // EntityHandle meshset;
514  // CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
515  // CHKERR mField.get_moab().add_entities(meshset, (*bc_ptr)[dd]);
516  // std::string file_name = disp_block_set_name +
517  // "_traction_free_bc_" + boost::lexical_cast<std::string>(dd) + ".vtk";
518  // CHKERR mField.get_moab().write_file(file_name.c_str(), " VTK ", "",
519  // &meshset, 1);
520  // CHKERR mField.get_moab().delete_entities(&meshset, 1);
521  // }
522 
524 }
525 
526 /**
527  * @brief Set integration rule on element
528  * \param order on row
529  * \param order on column
530  * \param order on data
531  *
532  * Use maximal oder on data in order to determine integration rule
533  *
534  */
535 struct VolRule {
536  int operator()(int p_row, int p_col, int p_data) const {
537  return 2 * (p_data + 1);
538  }
539 };
540 
541 struct FaceRule {
542  int operator()(int p_row, int p_col, int p_data) const {
543  return 2 * (p_data);
544  }
545 };
546 
548 
551 
553  BaseFunctionUnknownInterface **iface) const {
555  *iface = NULL;
556  if (uuid == IDD_TET_BASE_FUNCTION) {
557  *iface = const_cast<CGGUserPolynomialBase *>(this);
559  } else {
560  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "wrong interference");
561  }
562  CHKERR BaseFunction::query_interface(uuid, iface);
564  }
565 
567  boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
569 
571  CHKERR ctx_ptr->query_interface(IDD_TET_BASE_FUNCTION, &iface);
572  cTx = reinterpret_cast<EntPolynomialBaseCtx *>(iface);
573 
574  int nb_gauss_pts = pts.size2();
575  if (!nb_gauss_pts) {
577  }
578 
579  if (pts.size1() < 3) {
580  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
581  "Wrong dimension of pts, should be at least 3 rows with "
582  "coordinates");
583  }
584 
585  switch (cTx->sPace) {
586  case HDIV:
588  break;
589  default:
590  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
591  }
592 
594  }
595 
596 private:
598 
600 
603 
604  // This should be used only in case USER_BASE is selected
605  if (cTx->bAse != USER_BASE) {
606  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
607  "Wrong base, should be USER_BASE");
608  }
609  // get access to data structures on element
611  // Get approximation order on element. Note that bubble functions are only
612  // on tetrahedron.
613  const int order = data.dataOnEntities[MBTET][0].getDataOrder();
614  /// number of integration points
615  const int nb_gauss_pts = pts.size2();
616  // number of base functions
617 
618  // calculate shape functions, i.e. barycentric coordinates
619  shapeFun.resize(nb_gauss_pts, 4, false);
620  CHKERR ShapeMBTET(&*shapeFun.data().begin(), &pts(0, 0), &pts(1, 0),
621  &pts(2, 0), nb_gauss_pts);
622  // direvatives of shape functions
623  double diff_shape_fun[12];
624  CHKERR ShapeDiffMBTET(diff_shape_fun);
625 
626  const int nb_base_functions = NBVOLUMETET_CCG_BUBBLE(order);
627  // get base functions and set size
628  MatrixDouble &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
629  phi.resize(nb_gauss_pts, 9 * nb_base_functions, false);
630  // finally calculate base functions
632  &phi(0, 0), &phi(0, 1), &phi(0, 2),
633 
634  &phi(0, 3), &phi(0, 4), &phi(0, 5),
635 
636  &phi(0, 6), &phi(0, 7), &phi(0, 8));
637  CHKERR CGG_BubbleBase_MBTET(order, &shapeFun(0, 0), diff_shape_fun, t_phi,
638  nb_gauss_pts);
639 
641  }
642 };
643 
645  const int tag, const bool do_rhs, const bool do_lhs,
646  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe) {
648  fe = boost::make_shared<EpElement<VolumeElementForcesAndSourcesCore>>(mField);
649 
650  fe->getUserPolynomialBase() =
651  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
652  fe->getOpPtrVector().push_back(new OpL2Transform());
653 
654  // set integration rule
655  fe->getRuleHook = VolRule();
656 
657  if (!dataAtPts) {
658  dataAtPts =
659  boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
660  dataAtPts->physicsPtr = physicalEquations;
661  }
662 
663  // calculate fields values
664  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
665  piolaStress, dataAtPts->getApproxPAtPts()));
666  fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
667  bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
668  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
669  piolaStress, dataAtPts->getDivPAtPts()));
670  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
671  eshelbyStress, dataAtPts->getApproxSigmaAtPts()));
672  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
673  eshelbyStress, dataAtPts->getDivSigmaAtPts()));
674  fe->getOpPtrVector().push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
675  streachTensor, dataAtPts->getLogStreachTensorAtPts(), MBTET));
676  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
677  rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
678  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
679  rotAxis, dataAtPts->getRotAxisDotAtPts(), MBTET));
680  fe->getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
681  materialGradient, dataAtPts->getBigGAtPts(), MBTET));
682  fe->getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
683  materialGradient + "0", dataAtPts->getBigG0AtPts(), MBTET));
684  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
685  spatialDisp, dataAtPts->getSmallWAtPts(), MBTET));
686 
687  // velocities
688  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
689  spatialDisp, dataAtPts->getSmallWDotAtPts(), MBTET));
690  fe->getOpPtrVector().push_back(
692  streachTensor, dataAtPts->getLogStreachDotTensorAtPts(), MBTET));
693 
694  // calculate other derived quantities
695  fe->getOpPtrVector().push_back(
697 
698  // evaluate integration points
699  fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
700  spatialDisp, tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
701 
703 }
704 
706  const int tag, const bool add_elastic, const bool add_material,
707  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe_rhs,
708  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe_lhs) {
710 
711  // Right hand side
712  CHKERR setBaseVolumeElementOps(tag, true, false, fe_rhs);
713 
714  // elastic
715  if (add_elastic) {
716  fe_rhs->getOpPtrVector().push_back(
718  fe_rhs->getOpPtrVector().push_back(
720  fe_rhs->getOpPtrVector().push_back(
722  fe_rhs->getOpPtrVector().push_back(
724  fe_rhs->getOpPtrVector().push_back(
726  fe_rhs->getOpPtrVector().push_back(
728  }
729 
730  // Left hand side
731  CHKERR setBaseVolumeElementOps(tag, true, true, fe_lhs);
732 
733  // elastic
734  if (add_elastic) {
735 
736  // Schur
737  fe_lhs->getOpPtrVector().push_back(
739 
740  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_du(
742  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
744 
745  fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
747  fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
749 
750  fe_lhs->getOpPtrVector().push_back(
752  fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
754 
755  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
757  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
758  streachTensor, rotAxis, dataAtPts, false));
759 
760  fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dP(
761  rotAxis, piolaStress, dataAtPts, false));
762  fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dBubble(
763  rotAxis, bubbleField, dataAtPts, false));
764  fe_lhs->getOpPtrVector().push_back(
766 
767  // Schur
768  dataAtPts->ooMatPtr = boost::make_shared<MatrixDouble>();
769  fe_lhs->getOpPtrVector().push_back(
771  if (alpha_w < std::numeric_limits<double>::epsilon()) {
772  dataAtPts->wwMatPtr = boost::make_shared<MatrixDouble>();
773  fe_lhs->getOpPtrVector().push_back(
775  } else {
776  dataAtPts->wwMatPtr.reset();
777  }
778  fe_lhs->getOpPtrVector().push_back(
780 
781  if (add_material) {
782  }
783  }
784 
786 }
787 
789  const bool add_elastic, const bool add_material,
790  boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> &fe_rhs,
791  boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> &fe_lhs) {
793 
794  fe_rhs =
795  boost::make_shared<EpElement<FaceElementForcesAndSourcesCore>>(mField);
796  fe_lhs =
797  boost::make_shared<EpElement<FaceElementForcesAndSourcesCore>>(mField);
798 
799  // set integration rule
800  fe_rhs->getRuleHook = FaceRule();
801  fe_lhs->getRuleHook = FaceRule();
802 
803  if (add_elastic) {
804  fe_rhs->getOpPtrVector().push_back(
806  fe_rhs->getOpPtrVector().push_back(
808 
809  }
810 
812 }
813 
817  elasticFeLhs);
820 }
821 
824  boost::shared_ptr<FEMethod> null;
825  boost::shared_ptr<EpElement<FeTractionBc>> spatial_traction_bc(
827 
828  CHKERR DMMoFEMTSSetIFunction(dm, DM_NO_ELEMENT, null, spatial_traction_bc,
829  null);
833  spatial_traction_bc);
834 
835  schurAssembly = boost::make_shared<EpFEMethod>();
841 }
842 
845  boost::shared_ptr<TsCtx> ts_ctx;
847  CHKERR TSSetIFunction(ts, f, PETSC_NULL, PETSC_NULL);
848  CHKERR TSSetIJacobian(ts, m, m, PETSC_NULL, PETSC_NULL);
849  CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx.get(), PETSC_NULL);
850 
851  auto add_schur_streach_op = [this](auto &list, Mat S, AO aoS) {
853  for (auto &fe : list)
854  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
855  fe_cast->addStreachSchurMatrix(S, aoS);
856  else
857  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
859  };
860 
861  auto add_schur_streach_pre = [this](auto &list, Mat S, AO aoS) {
863  for (auto &fe : list)
864  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
865  fe_cast->addStreachSchurMatrix(S, aoS);
866  else
867  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
869  };
870 
871  auto add_schur_bubble_op = [this](auto &list, Mat S, AO aoS) {
873  for (auto &fe : list)
874  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
875  fe_cast->addBubbleSchurMatrix(S, aoS);
876  else
877  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
879  };
880 
881  auto add_schur_bubble_pre = [this](auto &list, Mat S, AO aoS) {
883  for (auto &fe : list)
884  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
885  fe_cast->addBubbleSchurMatrix(S, aoS);
886  else
887  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
889  };
890 
891  auto add_schur_spatial_disp_op = [this](auto &list, Mat S, AO aoS) {
893  for (auto &fe : list)
894  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
895  fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
896  else
897  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
899  };
900 
901  auto add_schur_spatial_disp_pre = [this](auto &list, Mat S, AO aoS) {
903  for (auto &fe : list)
904  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
905  fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
906  else
907  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
909  };
910 
911  auto add_schur_omega_op = [this](auto &list, Mat S, AO aoS) {
913  for (auto &fe : list)
914  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
915  fe_cast->addOmegaSchurMatrix(S, aoS);
916  else
917  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
919  };
920 
921  auto add_schur_omega_pre = [this](auto &list, Mat S, AO ao) {
923  for (auto &fe : list)
924  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
925  fe_cast->addOmegaSchurMatrix(S, ao);
926  else
927  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
929  };
930 
931  const MoFEM::Problem *schur_streach_prb_ptr;
932  CHKERR DMMoFEMGetProblemPtr(dmElasticSchurStreach, &schur_streach_prb_ptr);
933  if (auto sub_data = schur_streach_prb_ptr->subProblemData) {
934  Mat Suu;
935  CHKERR DMCreateMatrix(dmElasticSchurStreach, &Suu);
936  AO aoSuu;
937  CHKERR sub_data->getRowMap(&aoSuu);
938 
939  CHKERR add_schur_streach_op(ts_ctx->loops_to_do_IJacobian, Suu, aoSuu);
940  CHKERR add_schur_streach_pre(ts_ctx->preProcess_IJacobian, Suu, aoSuu);
941  CHKERR add_schur_streach_pre(ts_ctx->postProcess_IJacobian, Suu, aoSuu);
942 
943  CHKERR MatDestroy(&Suu);
944  CHKERR AODestroy(&aoSuu);
945 
946  const MoFEM::Problem *schur_bubble_prb_ptr;
947  CHKERR DMMoFEMGetProblemPtr(dmElasticSchurBubble, &schur_bubble_prb_ptr);
948  if (auto bubble_data = schur_bubble_prb_ptr->subProblemData) {
949  Mat SBubble;
950  CHKERR DMCreateMatrix(dmElasticSchurBubble, &SBubble);
951  AO aoSBubble;
952  CHKERR bubble_data->getRowMap(&aoSBubble);
953 
954  CHKERR add_schur_bubble_op(ts_ctx->loops_to_do_IJacobian, SBubble,
955  aoSBubble);
956  CHKERR add_schur_bubble_pre(ts_ctx->preProcess_IJacobian, SBubble,
957  aoSBubble);
958  CHKERR add_schur_bubble_pre(ts_ctx->postProcess_IJacobian, SBubble,
959  aoSBubble);
960 
961  CHKERR MatDestroy(&SBubble);
962  CHKERR AODestroy(&aoSBubble);
963 
964  const MoFEM::Problem *schur_omega_prb_ptr;
965  CHKERR DMMoFEMGetProblemPtr(dmElasticSchurOmega, &schur_omega_prb_ptr);
966  if (auto tet_stress_data = schur_omega_prb_ptr->subProblemData) {
967  Mat SOmega;
968  CHKERR DMCreateMatrix(dmElasticSchurOmega, &SOmega);
969  AO aoSOmega;
970  CHKERR tet_stress_data->getRowMap(&aoSOmega);
971 
972  CHKERR add_schur_omega_op(ts_ctx->loops_to_do_IJacobian, SOmega,
973  aoSOmega);
974  CHKERR add_schur_omega_pre(ts_ctx->preProcess_IJacobian, SOmega,
975  aoSOmega);
976  CHKERR add_schur_omega_pre(ts_ctx->postProcess_IJacobian, SOmega,
977  aoSOmega);
978 
979  const MoFEM::Problem *schur_spatial_disp_prb_ptr;
981  &schur_spatial_disp_prb_ptr);
982  if (auto spatial_disp_data =
983  schur_spatial_disp_prb_ptr->subProblemData) {
984 
985  Mat Sw;
986  CHKERR DMCreateMatrix(dmElasticSchurSpatialDisp, &Sw);
987  AO aoSw;
988  CHKERR spatial_disp_data->getRowMap(&aoSw);
989 
990  CHKERR add_schur_spatial_disp_op(ts_ctx->loops_to_do_IJacobian, Sw,
991  aoSw);
992  CHKERR add_schur_spatial_disp_pre(ts_ctx->preProcess_IJacobian, Sw,
993  aoSw);
994  CHKERR add_schur_spatial_disp_pre(ts_ctx->postProcess_IJacobian, Sw,
995  aoSw);
996 
997  CHKERR MatDestroy(&Sw);
998  CHKERR AODestroy(&aoSw);
999  } else
1001  "Problem does not have sub-problem data");
1002 
1003  CHKERR MatDestroy(&SOmega);
1004  CHKERR AODestroy(&aoSOmega);
1005 
1006  } else
1008  "Problem does not have sub-problem data");
1009 
1010  } else
1012  "Problem does not have sub-problem data");
1013 
1014  } else
1016  "Problem does not have sub-problem data");
1017 
1018  struct Monitor : public FEMethod {
1019 
1020  using Ele = ForcesAndSourcesCore;
1023  using SetPtsData = FieldEvaluatorInterface::SetPtsData;
1024 
1025  EshelbianCore &eP;
1026  boost::shared_ptr<SetPtsData> dataFieldEval;
1027  boost::shared_ptr<VolEle> volPostProcEle;
1028  boost::shared_ptr<double> gEnergy;
1029 
1030  Monitor(EshelbianCore &ep)
1031  : eP(ep),
1032  dataFieldEval(ep.mField.getInterface<FieldEvaluatorInterface>()
1033  ->getData<VolEle>()),
1034  volPostProcEle(new VolEle(ep.mField)), gEnergy(new double) {
1035  ierr = ep.mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
1036  dataFieldEval, "EP");
1037  CHKERRABORT(PETSC_COMM_WORLD, ierr);
1038 
1039  auto no_rule = [](int, int, int) { return -1; };
1040 
1041  auto set_element_for_field_eval = [&]() {
1042  boost::shared_ptr<Ele> vol_ele(dataFieldEval->feMethodPtr.lock());
1043  vol_ele->getRuleHook = no_rule;
1044  vol_ele->getUserPolynomialBase() =
1045  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1046  vol_ele->getOpPtrVector().push_back(new OpL2Transform());
1047 
1048  vol_ele->getOpPtrVector().push_back(
1050  eP.piolaStress, eP.dataAtPts->getApproxPAtPts()));
1051  vol_ele->getOpPtrVector().push_back(
1053  eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
1054  vol_ele->getOpPtrVector().push_back(
1056  eP.streachTensor, eP.dataAtPts->getLogStreachTensorAtPts(), MBTET));
1057  vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1058  eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
1059  vol_ele->getOpPtrVector().push_back(
1061  eP.materialGradient, eP.dataAtPts->getBigGAtPts(), MBTET));
1062  vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1063  eP.spatialDisp, eP.dataAtPts->getSmallWAtPts(), MBTET));
1064  vol_ele->getOpPtrVector().push_back(
1066  eP.dataAtPts));
1067  };
1068 
1069 
1070  auto set_element_for_post_process = [&]() {
1071  volPostProcEle->getRuleHook = VolRule();
1072  volPostProcEle->getUserPolynomialBase() =
1073  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1074  volPostProcEle->getOpPtrVector().push_back(new OpL2Transform());
1075 
1076  volPostProcEle->getOpPtrVector().push_back(
1078  eP.piolaStress, eP.dataAtPts->getApproxPAtPts()));
1079  volPostProcEle->getOpPtrVector().push_back(
1081  eP.bubbleField, eP.dataAtPts->getApproxPAtPts(), MBMAXTYPE));
1082  volPostProcEle->getOpPtrVector().push_back(
1084  eP.streachTensor, eP.dataAtPts->getLogStreachTensorAtPts(), MBTET));
1085  volPostProcEle->getOpPtrVector().push_back(
1087  eP.rotAxis, eP.dataAtPts->getRotAxisAtPts(), MBTET));
1088  volPostProcEle->getOpPtrVector().push_back(
1090  eP.materialGradient, eP.dataAtPts->getBigGAtPts(), MBTET));
1091  volPostProcEle->getOpPtrVector().push_back(
1093  eP.spatialDisp, eP.dataAtPts->getSmallWAtPts(), MBTET));
1094  volPostProcEle->getOpPtrVector().push_back(
1096  eP.dataAtPts));
1097  volPostProcEle->getOpPtrVector().push_back(
1098  new OpCalculateStrainEnergy(eP.spatialDisp, eP.dataAtPts, gEnergy));
1099  };
1100 
1101  set_element_for_field_eval();
1102  set_element_for_post_process();
1103  }
1104 
1105  MoFEMErrorCode preProcess() { return 0; }
1106 
1107  MoFEMErrorCode operator()() { return 0; }
1108 
1109  MoFEMErrorCode postProcess() {
1111 
1112  auto get_str_time = [](auto ts_t) {
1113  std::ostringstream ss;
1114  ss << boost::str(boost::format("%d") %
1115  static_cast<int>(std::ceil(ts_t * 1e6)));
1116  std::string s = ss.str();
1117  return s;
1118  };
1119 
1120  PetscViewer viewer;
1121  CHKERR PetscViewerBinaryOpen(
1122  PETSC_COMM_WORLD, ("restart_" + get_str_time(ts_t) + ".dat").c_str(),
1123  FILE_MODE_WRITE, &viewer);
1124  CHKERR VecView(ts_u, viewer);
1125  CHKERR PetscViewerDestroy(&viewer);
1126 
1127  CHKERR eP.postProcessResults(1, "out_sol_elastic_" + get_str_time(ts_t) +
1128  ".h5m");
1129 
1130  // Loop boundary elements with traction boundary conditions
1131  *gEnergy = 0;
1132  CHKERR eP.mField.loop_finite_elements(problemPtr->getName(), "EP",
1133  *volPostProcEle,nullptr);
1134 
1135  double body_energy;
1136  MPI_Allreduce(gEnergy.get(), &body_energy, 1, MPI_DOUBLE, MPI_SUM,
1137  eP.mField.get_comm());
1138  CHKERR PetscPrintf(eP.mField.get_comm(),
1139  "Step %d time %3.4g strain energy %3.6e\n", ts_step,
1140  ts_t, body_energy);
1141 
1142  auto post_proc_at_points = [&](std::array<double, 3> point,
1143  std::string str) {
1145 
1146  dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
1147 
1148  struct OpPrint : public VolOp {
1149 
1150  EshelbianCore &eP;
1151  std::array<double, 3> point;
1152  std::string str;
1153 
1154  OpPrint(EshelbianCore &ep, std::array<double, 3> &point,
1155  std::string &str)
1156  : VolOp(ep.spatialDisp, VolOp::OPROW), eP(ep), point(point),
1157  str(str) {}
1158 
1159  MoFEMErrorCode doWork(int side, EntityType type,
1162  if (type == MBTET) {
1163  if (getGaussPts().size2()) {
1164 
1165  auto t_h = getFTensor2FromMat<3, 3>(eP.dataAtPts->hAtPts);
1166  auto t_approx_P =
1167  getFTensor2FromMat<3, 3>(eP.dataAtPts->approxPAtPts);
1168 
1172  const double jac = dEterminant(t_h);
1174  t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
1175 
1176  auto add = [&]() {
1177  std::ostringstream s;
1178  s << str << " elem " << getFEEntityHandle() << " ";
1179  return s.str();
1180  };
1181 
1182  std::ostringstream print;
1183  print << add() << "comm rank " << eP.mField.get_comm_rank()
1184  << std::endl;
1185  print << add() << "point " << getVectorAdaptor(point.data(), 3)
1186  << std::endl;
1187  print << add() << "coords at gauss pts " << getCoordsAtGaussPts()
1188  << std::endl;
1189  print << add() << "w " << eP.dataAtPts->wAtPts << std::endl;
1190  print << add() << "Piola " << eP.dataAtPts->approxPAtPts
1191  << std::endl;
1192  print << add() << "Cauchy " << t_cauchy << std::endl;
1193  print << std::endl;
1194  CHKERR PetscSynchronizedPrintf(eP.mField.get_comm(), "%s",
1195  print.str().c_str());
1196 
1197  }
1198  }
1200  }
1201  };
1202 
1203  if (auto fe_ptr = dataFieldEval->feMethodPtr.lock()) {
1204 
1205  fe_ptr->getOpPtrVector().push_back(new OpPrint(eP, point, str));
1207  ->evalFEAtThePoint3D(point.data(), 1e-12, problemPtr->getName(),
1208  "EP", dataFieldEval,
1209  eP.mField.get_comm_rank(),
1211  fe_ptr->getOpPtrVector().pop_back();
1212  }
1213 
1215  };
1216 
1217  // Points for Cook beam
1218  std::array<double, 3> pointA = {48., 60., 5.};
1219  CHKERR post_proc_at_points(pointA, "Point A");
1220  PetscSynchronizedFlush(eP.mField.get_comm(), PETSC_STDOUT);
1221 
1222  std::array<double, 3> pointB = {48. / 2., 44. + (60. - 44.) / 2., 0.};
1223  CHKERR post_proc_at_points(pointB, "Point B");
1224  PetscSynchronizedFlush(eP.mField.get_comm(), PETSC_STDOUT);
1225 
1226  std::array<double, 3> pointC = {48. / 2., (44. - 0.) / 2., 0.};
1227  CHKERR post_proc_at_points(pointC, "Point C");
1228  PetscSynchronizedFlush(eP.mField.get_comm(), PETSC_STDOUT);
1229 
1231  }
1232  };
1233 
1234  boost::shared_ptr<FEMethod> monitor_ptr(new Monitor(*this));
1235  ts_ctx->get_loops_to_do_Monitor().push_back(
1237 
1238  CHKERR TSAppendOptionsPrefix(ts, "elastic_");
1239  CHKERR TSSetFromOptions(ts);
1241 }
1242 
1245 
1246  CHKERR TSSetDM(ts, dmElastic);
1247 
1248  SNES snes;
1249  CHKERR TSGetSNES(ts, &snes);
1250 
1251  PetscViewerAndFormat *vf;
1252  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
1253  PETSC_VIEWER_DEFAULT, &vf);
1254  CHKERR SNESMonitorSet(
1255  snes,
1256  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
1257  vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
1258 
1259  PetscSection section;
1260  CHKERR DMGetDefaultSection(dmElastic, &section);
1261  int num_fields;
1262  CHKERR PetscSectionGetNumFields(section, &num_fields);
1263  for (int ff = 0; ff != num_fields; ff++) {
1264  const char *field_name;
1265  CHKERR PetscSectionGetFieldName(section, ff, &field_name);
1266  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Field %d name %s\n", ff, field_name);
1267  }
1268 
1269  CHKERR DMoFEMMeshToLocalVector(dmElastic, x, INSERT_VALUES, SCATTER_FORWARD);
1270 
1271  // PetscRandom rctx;
1272  // PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
1273  // VecSetRandom(x, rctx);
1274  // VecScale(x, -1e-1);
1275  // PetscRandomDestroy(&rctx);
1276  // CHKERR VecView(x, PETSC_VIEWER_STDOUT_WORLD);
1277 
1278  CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
1279  CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
1280 
1281  // CHKERR VecView(x, PETSC_VIEWER_STDOUT_WORLD);
1282 
1283  // Adding field split solver
1284 
1285  KSP ksp;
1286  CHKERR SNESGetKSP(snes, &ksp);
1287  PC pc;
1288  CHKERR KSPGetPC(ksp, &pc);
1289  PetscBool is_uu_field_split;
1290  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_uu_field_split);
1291  if (is_uu_field_split) {
1292 
1293  const MoFEM::Problem *schur_uu_ptr;
1295  if (auto uu_data = schur_uu_ptr->subProblemData) {
1296 
1297  const MoFEM::Problem *prb_ptr;
1299  map<std::string, IS> is_map;
1300  for (int ff = 0; ff != num_fields; ff++) {
1301  const char *field_name;
1302  CHKERR PetscSectionGetFieldName(section, ff, &field_name);
1303  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1304  prb_ptr->getName(), ROW, field_name, 0, MAX_DOFS_ON_ENTITY,
1305  &is_map[field_name]);
1306  }
1307  // CHKERR mField.getInterface<ISManager>()
1308  // ->isCreateProblemFieldAndEntityType(
1309  // prb_ptr->getName(), ROW, piolaStress, MBTET, MBTET, 0,
1310  // MAX_DOFS_ON_ENTITY, &is_map["T_STRESS_ON_TETS"]);
1311  // CHKERR mField.getInterface<ISManager>()
1312  // ->isCreateProblemFieldAndEntityType(
1313  // prb_ptr->getName(), ROW, piolaStress, MBTRI, MBTRI, 0,
1314  // MAX_DOFS_ON_ENTITY, &is_map["T_STRESS_ON_TRIS"]);
1315 
1316  CHKERR uu_data->getRowIs(&is_map["E_IS_SUU"]);
1317 
1318  CHKERR PCFieldSplitSetIS(pc, NULL, is_map[streachTensor]);
1319  CHKERR PCFieldSplitSetIS(pc, NULL, is_map["E_IS_SUU"]);
1320 
1321  CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER,
1322  schurAssembly->Suu);
1323 
1324  CHKERR PCSetUp(pc);
1325  PetscInt n;
1326  KSP *uu_ksp;
1327  CHKERR PCFieldSplitGetSubKSP(pc, &n, &uu_ksp);
1328  PC bubble_pc;
1329  CHKERR KSPGetPC(uu_ksp[1], &bubble_pc);
1330  PetscBool is_bubble_field_split;
1331  PetscObjectTypeCompare((PetscObject)bubble_pc, PCFIELDSPLIT,
1332  &is_bubble_field_split);
1333  if (is_bubble_field_split) {
1334 
1335  const MoFEM::Problem *schur_bubble_ptr;
1336  CHKERR DMMoFEMGetProblemPtr(dmElasticSchurBubble, &schur_bubble_ptr);
1337  if (auto bubble_data = schur_bubble_ptr->subProblemData) {
1338 
1339  CHKERR bubble_data->getRowIs(&is_map["E_IS_BUBBLE"]);
1340 
1341  AO uu_ao;
1342  CHKERR uu_data->getRowMap(&uu_ao);
1343 
1344  CHKERR AOApplicationToPetscIS(uu_ao, is_map[bubbleField]);
1345  CHKERR PCFieldSplitSetIS(bubble_pc, NULL, is_map[bubbleField]);
1346  CHKERR PCFieldSplitSetIS(bubble_pc, NULL, is_map["E_IS_BUBBLE"]);
1347  CHKERR PCFieldSplitSetSchurPre(
1348  bubble_pc, PC_FIELDSPLIT_SCHUR_PRE_USER, schurAssembly->SBubble);
1349 
1350  CHKERR PCSetUp(bubble_pc);
1351  PetscInt bubble_n;
1352  KSP *bubble_ksp;
1353  CHKERR PCFieldSplitGetSubKSP(bubble_pc, &bubble_n, &bubble_ksp);
1354  PC omega_pc;
1355  CHKERR KSPGetPC(bubble_ksp[1], &omega_pc);
1356  PetscBool is_omega_field_split;
1357  PetscObjectTypeCompare((PetscObject)omega_pc, PCFIELDSPLIT,
1358  &is_omega_field_split);
1359 
1360  if (is_omega_field_split) {
1361 
1362  const MoFEM::Problem *schur_omega_ptr;
1363  CHKERR DMMoFEMGetProblemPtr(dmElasticSchurOmega, &schur_omega_ptr);
1364  if (auto omega_data = schur_omega_ptr->subProblemData) {
1365 
1366  AO bubble_ao;
1367  CHKERR bubble_data->getRowMap(&bubble_ao);
1368 
1369  CHKERR AOApplicationToPetscIS(uu_ao, is_map[rotAxis]);
1370  CHKERR AOApplicationToPetscIS(bubble_ao, is_map[rotAxis]);
1371  CHKERR omega_data->getRowIs(&is_map["E_IS_OMEGA"]);
1372 
1373  CHKERR PCFieldSplitSetIS(omega_pc, NULL, is_map[rotAxis]);
1374  CHKERR PCFieldSplitSetIS(omega_pc, NULL, is_map["E_IS_OMEGA"]);
1375  CHKERR PCFieldSplitSetSchurPre(omega_pc,
1376  PC_FIELDSPLIT_SCHUR_PRE_USER,
1377  schurAssembly->SOmega);
1378 
1379  CHKERR PCSetUp(omega_pc);
1380  PetscInt omega_n;
1381  KSP *omega_ksp;
1382  CHKERR PCFieldSplitGetSubKSP(omega_pc, &omega_n, &omega_ksp);
1383  PC w_pc;
1384  CHKERR KSPGetPC(omega_ksp[1], &w_pc);
1385  PetscBool is_w_field_split;
1386  PetscObjectTypeCompare((PetscObject)w_pc, PCFIELDSPLIT,
1387  &is_w_field_split);
1388  if (is_w_field_split) {
1389 
1390  const MoFEM::Problem *schur_w_ptr;
1392  &schur_w_ptr);
1393  if (auto w_data = schur_w_ptr->subProblemData) {
1394 
1395  AO omega_ao;
1396  CHKERR omega_data->getRowMap(&omega_ao);
1397 
1398  CHKERR AOApplicationToPetscIS(uu_ao, is_map[spatialDisp]);
1399  CHKERR AOApplicationToPetscIS(bubble_ao, is_map[spatialDisp]);
1400  CHKERR AOApplicationToPetscIS(omega_ao, is_map[spatialDisp]);
1401  CHKERR w_data->getRowIs(&is_map["E_IS_W"]);
1402 
1403  CHKERR PCFieldSplitSetIS(w_pc, NULL, is_map[spatialDisp]);
1404  CHKERR PCFieldSplitSetIS(w_pc, NULL, is_map["E_IS_W"]);
1405  CHKERR PCFieldSplitSetSchurPre(
1406  w_pc, PC_FIELDSPLIT_SCHUR_PRE_USER, schurAssembly->Sw);
1407 
1408  CHKERR AODestroy(&omega_ao);
1409  }
1410  }
1411 
1412  CHKERR PetscFree(omega_ksp);
1413  }
1414  }
1415  CHKERR PetscFree(bubble_ksp);
1416  CHKERR AODestroy(&uu_ao);
1417  }
1418  }
1419  CHKERR PetscFree(uu_ksp);
1420 
1421  for (auto &m : is_map)
1422  CHKERR ISDestroy(&m.second);
1423  }
1424  }
1425 
1426  CHKERR TSSolve(ts, x);
1428 }
1429 
1431  const std::string file) {
1433 
1434  if (!dataAtPts) {
1435  dataAtPts =
1436  boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
1437  }
1438 
1440  post_proc.getUserPolynomialBase() =
1441  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1442  post_proc.getOpPtrVector().push_back(new OpL2Transform());
1443 
1445  post_proc.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
1446  piolaStress, dataAtPts->getApproxPAtPts()));
1447  post_proc.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
1448  bubbleField, dataAtPts->getApproxPAtPts(), MBMAXTYPE));
1449  post_proc.getOpPtrVector().push_back(
1451  streachTensor, dataAtPts->getLogStreachTensorAtPts(), MBTET));
1452  post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1453  rotAxis, dataAtPts->getRotAxisAtPts(), MBTET));
1454  post_proc.getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
1455  materialGradient, dataAtPts->getBigGAtPts(), MBTET));
1456  post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1457  spatialDisp, dataAtPts->getSmallWAtPts(), MBTET));
1458 
1459  // evaluate derived quantities
1460  post_proc.getOpPtrVector().push_back(
1462 
1463  // evaluate integration points
1464  post_proc.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
1465  spatialDisp, tag, true, false, dataAtPts, physicalEquations));
1466  post_proc.getOpPtrVector().push_back(new OpPostProcDataStructure(
1467  post_proc.postProcMesh, post_proc.mapGaussPts, spatialDisp, dataAtPts));
1468 
1469  CHKERR DMoFEMLoopFiniteElements(dM, elementVolumeName.c_str(), &post_proc);
1470  CHKERR post_proc.writeFile(file.c_str());
1472 }
1473 
1474 } // namespace EshelbianPlasticity
SmartPetscObj< DM > dmElastic
Elastic problem.
MoFEMErrorCode setGenericFaceElementOps(const bool add_elastic, const bool add_material, boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore >> &fe_rhs, boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore >> &fe_lhs)
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, BaseFunctionUnknownInterface **iface) const
structure for User Loop Methods on finite elementsIt can be used to calculate stiffness matrices,...
Get time direvatives of values at integration pts for tensor filed rank 1, i.e. vector field.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:22
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
Deprecated interface functions.
MoFEM interface unique ID.
Problem manager is used to build and partition problems.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEMErrorCode setBaseVolumeElementOps(const int tag, const bool do_rhs, const bool do_lhs, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore >> &fe)
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
field with continuous normal traction
Definition: definitions.h:179
user implemented approximation base
Definition: definitions.h:160
virtual moab::Interface & get_moab()=0
MoFEMErrorCode setElasticElementOps(const int tag)
boost::shared_ptr< EpFEMethod > schurAssembly
data structure for finite element entityIt keeps that about indices of degrees of freedom,...
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:318
MoFEM::VolumeElementForcesAndSourcesCore VolEle
SmartPetscObj< DM > dmElasticSchurOmega
Sub problem of dmElastic Schur.
boost::shared_ptr< BcRotVec > bcSpatialRotationVecPtr
virtual bool check_finite_element(const std::string &name) const =0
Check if finite element is in database.
Section manager is used to create indexes and sectionsFIXME: ISManager is not properly testsed by ato...
Definition: ISManager.hpp:36
moab::Interface & postProcMesh
std::array< boost::ptr_vector< EntData >, MBMAXTYPE > dataOnEntities
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:507
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:445
boost::shared_ptr< SubProblemData > subProblemData
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problemIt if true is assumed that matrix has the same indexing on rows and columns....
Definition: DMMMoFEM.cpp:378
MoFEMErrorCode getValueHdivForCGGBubble(MatrixDouble &pts)
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
MoFEMErrorCode addDMs(const BitRefLevel bit=BitRefLevel().set(0))
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:425
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
auto & getUserPolynomialBase()
Get the User Polynomial Base object.
Class used to pass element data to calculate base functions on tet,triangle,edge.
const FieldApproximationBase bAse
base class for all interface classes
FTensor::Index< 'n', 2 > n
Definition: PlasticOps.hpp:68
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:550
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:202
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:514
keeps basic data about problemThis is low level structure with information about problem,...
MoFEMErrorCode removeDofsOnEntities(const std::string problem_name, const std::string field_name, const Range ents, const int lo_coeff=0, const int hi_coeff=MAX_DOFS_ON_ENTITY, int verb=VERBOSE, const bool debug=false)
Remove DOFs from problemRemove DOFs from problem which are on entities on the given range and given f...
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition: DMMMoFEM.cpp:772
Get values at integration pts for tensor filed rank 1, i.e. vector field.
static const MOFEMuuid IDD_TET_BASE_FUNCTION
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Field evaluator interface.
MoFEMErrorCode addFields(const EntityHandle meshset=0)
static const MOFEMuuid IDD_MOFEMEshelbianCrackInterface
Implementation of tonsorial bubble base div(v) = 0.
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMMoFEM.cpp:177
MoFEMErrorCode writeFile(const std::string file_name)
wrote results in (MOAB) format, use "file_name.h5m"
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
Definition: DMMMoFEM.cpp:198
const int dim
int operator()(int p_row, int p_col, int p_data) const
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
virtual int get_comm_rank() const =0
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:331
ForcesAndSourcesCore::UserDataOperator UserDataOperator
boost::shared_ptr< BcDispVec > bcSpatialDispVecPtr
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
Base class if inherited used to calculate base functions.
SmartPetscObj< DM > dmElasticSchurStreach
Sub problem of dmElastic Schur.
MoFEMErrorCode postProcessResults(const int tag, const std::string file)
MoFEMErrorCode addBoundaryFiniteElement(const EntityHandle meshset=0)
MoFEMErrorCode addVolumeFiniteElement(const EntityHandle meshset=0)
static double dEterminant(T &t)
Calculate the determinant of a 3x3 matrix or a tensor of rank 2.
Definition: Templates.hpp:452
DataForcesAndSourcesCore & dAta
MoFEMErrorCode setGenericVolumeElementOps(const int tag, const bool add_elastic, const bool add_material, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore >> &fe_rhs, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore >> &fe_lhs)
FTensor::Index< 'm', 2 > m
Definition: PlasticOps.hpp:67
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.
SmartPetscObj< DM > dmElasticSchurSpatialDisp
Sub problem of dmElastic Schur.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
boost::shared_ptr< TractionFreeBc > bcSpatialFreeTraction
TractionBc(std::string name, std::vector< double > &attr, Range &faces)
MoFEMErrorCode setUpTSElastic(TS ts, Mat m, Vec f)
MoFEMErrorCode CGG_BubbleBase_MBTET(const int p, const double *N, const double *diffN, FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > &phi, const int gdim)
Calculate CGGT tonsorial bubble base.
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
BcRot(std::string name, std::vector< double > &attr, Range &faces)
boost::shared_ptr< DataAtIntegrationPts > dataAtPts
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:42
Post processing.
MoFEMErrorCode getTractionFreeBc(const EntityHandle meshset, boost::shared_ptr< TractionFreeBc > &bc_ptr, const std::string disp_block_set_name, const std::string rot_block_set_name)
boost::shared_ptr< PhysicalEquations > physicalEquations
MoFEMErrorCode solveElastic(TS ts, Vec x)
Calculate divergence of tonsorial field using vectorial base.
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition: DMMMoFEM.cpp:105
FieldOrderTable & getFieldOrderTable()
Get the Field Order Table.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > elasticBcRhs
#define CHKERR
Inline error check.
Definition: definitions.h:602
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
std::vector< EntityHandle > mapGaussPts
constexpr int order
FTensor::Index< 'k', 2 > k
Definition: ContactOps.hpp:28
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:291
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:50
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Set integration rule on element.
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMMoFEM.cpp:1001
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > elasticFeLhs
SmartPetscObj< DM > dM
Coupled problem all fields.
Eshelbian plasticity interface.
boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > elasticBcLhs
SmartPetscObj< DM > dmElasticSchurBubble
Sub problem of dmElastic Schur.
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:304
boost::shared_ptr< TractionBcVec > bcSpatialTraction
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Index< 'j', 2 > j
Definition: ContactOps.hpp:27
structure to get information form mofem into DataForcesAndSourcesCore
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMMoFEM.cpp:360
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
BcDisp(std::string name, std::vector< double > &attr, Range &faces)
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMMoFEM.cpp:719
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:348
virtual MPI_Comm & get_comm() const =0
continuous field
Definition: definitions.h:177
std::string getName() const
virtual const Field * get_field_structure(const std::string &name)=0
get field structure
Calculate symmetric tensor field rates ant integratio pts.
auto createSmartDM
Creates smart DM object.
Definition: AuxPETSc.hpp:174
Calculate symmetric tensor field values at integration pts.
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:982
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
boost::shared_ptr< SPD > getData(const double *ptr=nullptr, const int nb_eval_points=0, const double eps=1e-12, VERBOSITY_LEVELS verb=QUIET)
Get the Data object.
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > elasticFeRhs
DEPRECATED MoFEMErrorCode loop_finite_elements(const Problem *problem_ptr, const std::string &fe_name, FEMethod &method, int lower_rank, int upper_rank, MoFEMTypes bh, int verb=DEFAULT_VERBOSITY)
int operator()(int p_row, int p_col, int p_data) const
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26
field with C-1 continuity
Definition: definitions.h:180