v0.9.0
EshelbianPlasticity.cpp
Go to the documentation of this file.
1 /**
2  * \file EshelbianPlasticity.cpp
3  * \example EshelbianPlasticity.cpp
4  *
5  * \brief Eshelbian plasticity implementation
6  */
7 
8 #include <MoFEM.hpp>
9 using namespace MoFEM;
10 
11 #include <BasicFiniteElements.hpp>
13 
14 #include <EshelbianPlasticity.hpp>
15 #include <boost/math/constants/constants.hpp>
16 
17 namespace EshelbianPlasticity {
18 
19 MoFEMErrorCode EshelbianCore::query_interface(const MOFEMuuid &uuid,
20  UnknownInterface **iface) const {
22  *iface = NULL;
24  *iface = const_cast<EshelbianCore *>(this);
26  }
27  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
29 }
30 
31 MoFEMErrorCode OpJacobian::doWork(int side, EntityType type, EntData &data) {
33  if (type != MBVERTEX)
35 
36  if (evalRhs)
37  CHKERR evaluateRhs(data);
38 
39  if (evalLhs)
40  CHKERR evaluateLhs(data);
41 
43 }
44 
45 EshelbianCore::EshelbianCore(MoFEM::Interface &m_field)
46  : mField(m_field), 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  // CHKERR add_field_to_fe(elementVolumeName, materialDisp);
247  // CHKERR add_field_to_fe(elementVolumeName, tauField);
248  // CHKERR add_field_to_fe(elementVolumeName, lambdaField);
249 
251  materialGradient + "0");
252  }
253 
254  // build finite elements data structures
256 
258 }
259 
263 
264  auto bc_elements_add_to_range = [&](const std::string disp_block_set_name,
265  Range &r) {
268  if (it->getName().compare(0, disp_block_set_name.length(),
269  disp_block_set_name) == 0) {
270  Range faces;
271  CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
272  true);
273  r.merge(faces);
274  }
275  }
277  };
278 
279  // set finite element fields
280  auto add_field_to_fe = [this](const std::string fe,
281  const std::string field_name) {
287  };
288 
289  Range natural_bc_elements;
290  CHKERR bc_elements_add_to_range("SPATIAL_DISP_BC", natural_bc_elements);
291  CHKERR bc_elements_add_to_range("SPATIAL_ROTATION_BC", natural_bc_elements);
292  Range essentail_bc_elements;
293  CHKERR bc_elements_add_to_range("SPATIAL_TRACTION_BC", essentail_bc_elements);
294 
296  CHKERR mField.add_ents_to_finite_element_by_type(natural_bc_elements, MBTRI,
298  CHKERR add_field_to_fe(naturalBcElement, piolaStress);
299  CHKERR add_field_to_fe(naturalBcElement, eshelbyStress);
301 
303  CHKERR mField.add_ents_to_finite_element_by_type(essentail_bc_elements, MBTRI,
305  CHKERR add_field_to_fe(essentialBcElement, piolaStress);
306  CHKERR add_field_to_fe(essentialBcElement, eshelbyStress);
308 
310 }
311 
314 
315  // find adjacencies between finite elements and dofs
317 
318  // Create coupled problem
319  dM = createSmartDM(mField.get_comm(), "DMMOFEM");
320  CHKERR DMMoFEMCreateMoFEM(dM, &mField, "ESHELBY_PLASTICITY", bit,
321  BitRefLevel().set());
322  CHKERR DMMoFEMSetDestroyProblem(dM, PETSC_TRUE);
323  CHKERR DMMoFEMSetIsPartitioned(dM, PETSC_TRUE);
327  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_TRUE;
328  CHKERR DMSetUp(dM);
329  mField.getInterface<ProblemsManager>()->buildProblemFromFields = PETSC_FALSE;
330 
331  auto remove_dofs_on_essential_spatial_stress_boundary =
332  [&](const std::string prb_name) {
334  for (int d : {0, 1, 2})
336  prb_name, piolaStress, (*bcSpatialFreeTraction)[d], d, d, NOISY,
337  true);
339  };
340  CHKERR remove_dofs_on_essential_spatial_stress_boundary("ESHELBY_PLASTICITY");
341 
342  // Create elastic sub-problem
343  dmElastic = createSmartDM(mField.get_comm(), "DMMOFEM");
344  CHKERR DMMoFEMCreateSubDM(dmElastic, dM, "ELASTIC_PROBLEM");
356  CHKERR DMSetUp(dmElastic);
357 
358  // Create elastic streach-problem
361  "ELASTIC_PROBLEM_STREACH_SCHUR");
372  CHKERR DMSetUp(dmElasticSchurStreach);
373 
374  // Create elastic bubble-problem
377  "ELASTIC_PROBLEM_BUBBLE_SCHUR");
387  CHKERR DMSetUp(dmElasticSchurBubble);
388 
389  // Create elastic omega-problem
392  "ELASTIC_PROBLEM_OMEGA_SCHUR");
401  CHKERR DMSetUp(dmElasticSchurOmega);
402 
403  // Create elastic tet_stress-problem
406  "ELASTIC_PROBLEM_SPATIAL_DISP_SCHUR");
410  elementVolumeName.c_str());
413  essentialBcElement.c_str());
417 
418  {
419  PetscSection section;
420  CHKERR mField.getInterface<ISManager>()->sectionCreate("ELASTIC_PROBLEM",
421  &section);
422  CHKERR DMSetDefaultSection(dmElastic, section);
423  CHKERR DMSetDefaultGlobalSection(dmElastic, section);
424  CHKERR PetscSectionDestroy(&section);
425  }
426 
428 }
429 
430 BcDisp::BcDisp(std::string name, std::vector<double> &attr, Range &faces)
431  : blockName(name), faces(faces) {
432  vals.resize(3, false);
433  flags.resize(3, false);
434  for (int ii = 0; ii != 3; ++ii) {
435  vals[ii] = attr[ii];
436  flags[ii] = static_cast<int>(attr[ii + 3]);
437  }
438 }
439 
440 BcRot::BcRot(std::string name, std::vector<double> &attr, Range &faces)
441  : blockName(name), faces(faces) {
442  vals.resize(3, false);
443  for (int ii = 0; ii != 3; ++ii) {
444  vals[ii] = attr[ii];
445  }
446  theta = attr[3];
447 }
448 
449 TractionBc::TractionBc(std::string name, std::vector<double> &attr,
450  Range &faces)
451  : blockName(name), faces(faces) {
452  vals.resize(3, false);
453  flags.resize(3, false);
454  for (int ii = 0; ii != 3; ++ii) {
455  vals[ii] = attr[ii];
456  flags[ii] = static_cast<int>(attr[ii + 3]);
457  }
458 }
459 
462  boost::shared_ptr<TractionFreeBc> &bc_ptr,
463  const std::string disp_block_set_name,
464  const std::string rot_block_set_name) {
466  Range tets;
467  CHKERR mField.get_moab().get_entities_by_type(meshset, MBTET, tets);
468  Range tets_skin_part;
469  Skinner skin(&mField.get_moab());
470  CHKERR skin.find_skin(0, tets, false, tets_skin_part);
471  ParallelComm *pcomm =
472  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
473  Range tets_skin;
474  CHKERR pcomm->filter_pstatus(tets_skin_part,
475  PSTATUS_SHARED | PSTATUS_MULTISHARED,
476  PSTATUS_NOT, -1, &tets_skin);
477 
478  bc_ptr->resize(3);
479  for (int dd = 0; dd != 3; ++dd)
480  (*bc_ptr)[dd] = tets_skin;
481 
483  if (it->getName().compare(0, disp_block_set_name.length(),
484  disp_block_set_name) == 0) {
485  std::vector<double> block_attributes;
486  CHKERR it->getAttributes(block_attributes);
487  if (block_attributes.size() != 6) {
488  SETERRQ2(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
489  "In block %s six attributes are required for given BC "
490  "blockset (3 values + "
491  "3 flags) != %d",
492  it->getName().c_str(), block_attributes.size());
493  }
494  Range faces;
495  CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
496  true);
497  if (block_attributes[3] != 0)
498  (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
499  if (block_attributes[4] != 0)
500  (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
501  if (block_attributes[5] != 0)
502  (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
503  }
504  if (it->getName().compare(0, rot_block_set_name.length(),
505  rot_block_set_name) == 0) {
506  Range faces;
507  CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 2, faces,
508  true);
509  (*bc_ptr)[0] = subtract((*bc_ptr)[0], faces);
510  (*bc_ptr)[1] = subtract((*bc_ptr)[1], faces);
511  (*bc_ptr)[2] = subtract((*bc_ptr)[2], faces);
512  }
513  }
514 
515  // for (int dd = 0; dd != 3; ++dd) {
516  // EntityHandle meshset;
517  // CHKERR mField.get_moab().create_meshset(MESHSET_SET, meshset);
518  // CHKERR mField.get_moab().add_entities(meshset, (*bc_ptr)[dd]);
519  // std::string file_name = disp_block_set_name +
520  // "_traction_free_bc_" + boost::lexical_cast<std::string>(dd) + ".vtk";
521  // CHKERR mField.get_moab().write_file(file_name.c_str(), " VTK ", "",
522  // &meshset, 1);
523  // CHKERR mField.get_moab().delete_entities(&meshset, 1);
524  // }
525 
527 }
528 
529 /**
530  * @brief Set integration rule on element
531  * \param order on row
532  * \param order on column
533  * \param order on data
534  *
535  * Use maximal oder on data in order to determine integration rule
536  *
537  */
538 struct VolRule {
539  int operator()(int p_row, int p_col, int p_data) const {
540  return 2 * (p_data + 1);
541  }
542 };
543 
544 struct FaceRule {
545  int operator()(int p_row, int p_col, int p_data) const {
546  return 2 * (p_data);
547  }
548 };
549 
551 
554 
556  BaseFunctionUnknownInterface **iface) const {
558  *iface = NULL;
559  if (uuid == IDD_TET_BASE_FUNCTION) {
560  *iface = const_cast<CGGUserPolynomialBase *>(this);
562  } else {
563  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "wrong interference");
564  }
565  CHKERR BaseFunction::query_interface(uuid, iface);
567  }
568 
570  boost::shared_ptr<BaseFunctionCtx> ctx_ptr) {
572 
574  CHKERR ctx_ptr->query_interface(IDD_TET_BASE_FUNCTION, &iface);
575  cTx = reinterpret_cast<EntPolynomialBaseCtx *>(iface);
576 
577  int nb_gauss_pts = pts.size2();
578  if (!nb_gauss_pts) {
580  }
581 
582  if (pts.size1() < 3) {
583  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
584  "Wrong dimension of pts, should be at least 3 rows with "
585  "coordinates");
586  }
587 
588  switch (cTx->sPace) {
589  case HDIV:
591  break;
592  default:
593  SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Not yet implemented");
594  }
595 
597  }
598 
599 private:
601 
603 
606 
607  // This should be used only in case USER_BASE is selected
608  if (cTx->bAse != USER_BASE) {
609  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
610  "Wrong base, should be USER_BASE");
611  }
612  // get access to data structures on element
614  // Get approximation order on element. Note that bubble functions are only
615  // on tetrahedron.
616  const int order = data.dataOnEntities[MBTET][0].getDataOrder();
617  /// number of integration points
618  const int nb_gauss_pts = pts.size2();
619  // number of base functions
620 
621  // calculate shape functions, i.e. barycentric coordinates
622  shapeFun.resize(nb_gauss_pts, 4, false);
623  CHKERR ShapeMBTET(&*shapeFun.data().begin(), &pts(0, 0), &pts(1, 0),
624  &pts(2, 0), nb_gauss_pts);
625  // direvatives of shape functions
626  double diff_shape_fun[12];
627  CHKERR ShapeDiffMBTET(diff_shape_fun);
628 
629  const int nb_base_functions = NBVOLUMETET_CCG_BUBBLE(order);
630  // get base functions and set size
631  MatrixDouble &phi = data.dataOnEntities[MBTET][0].getN(USER_BASE);
632  phi.resize(nb_gauss_pts, 9 * nb_base_functions, false);
633  // finally calculate base functions
635  &phi(0, 0), &phi(0, 1), &phi(0, 2),
636 
637  &phi(0, 3), &phi(0, 4), &phi(0, 5),
638 
639  &phi(0, 6), &phi(0, 7), &phi(0, 8));
640  CHKERR CGG_BubbleBase_MBTET(order, &shapeFun(0, 0), diff_shape_fun, t_phi,
641  nb_gauss_pts);
642 
644  }
645 };
646 
648  const int tag, const bool do_rhs, const bool do_lhs,
649  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe) {
651  fe = boost::make_shared<EpElement<VolumeElementForcesAndSourcesCore>>(mField);
652 
653  fe->getUserPolynomialBase() =
654  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
655 
656  // set integration rule
657  fe->getRuleHook = VolRule();
658 
659  if (!dataAtPts) {
660  dataAtPts =
661  boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
662  dataAtPts->approxPAtPts = boost::make_shared<MatrixDouble>();
663  dataAtPts->divPAtPts = boost::make_shared<MatrixDouble>();
664  dataAtPts->approxSigmaAtPts = boost::make_shared<MatrixDouble>();
665  dataAtPts->divSigmaAtPts = boost::make_shared<MatrixDouble>();
666  dataAtPts->wAtPts = boost::make_shared<MatrixDouble>();
667  dataAtPts->wDotAtPts = boost::make_shared<MatrixDouble>();
668  dataAtPts->rotAxisAtPts = boost::make_shared<MatrixDouble>();
669  dataAtPts->logStreachTensorAtPts = boost::make_shared<MatrixDouble>();
670  dataAtPts->streachTensorAtPts = boost::make_shared<MatrixDouble>();
671  dataAtPts->diffStreachTensorAtPts = boost::make_shared<MatrixDouble>();
672  dataAtPts->logStreachDotTensorAtPts = boost::make_shared<MatrixDouble>();
673  dataAtPts->rotMatAtPts = boost::make_shared<MatrixDouble>();
674  dataAtPts->diffRotMatAtPts = boost::make_shared<MatrixDouble>();
675  dataAtPts->hAtPts = boost::make_shared<MatrixDouble>();
676  dataAtPts->GAtPts = boost::make_shared<MatrixDouble>();
677  dataAtPts->G0AtPts = boost::make_shared<MatrixDouble>();
678  dataAtPts->physicsPtr = physicalEquations;
679  }
680 
681  // calculate fields values
682  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
683  piolaStress, dataAtPts->approxPAtPts));
684  fe->getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
685  bubbleField, dataAtPts->approxPAtPts, MBMAXTYPE));
686  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
687  piolaStress, dataAtPts->divPAtPts));
688  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
689  eshelbyStress, dataAtPts->approxSigmaAtPts));
690  fe->getOpPtrVector().push_back(new OpCalculateHVecTensorDivergence<3, 3>(
691  eshelbyStress, dataAtPts->divSigmaAtPts));
692  fe->getOpPtrVector().push_back(new OpCalculateTensor2SymmetricFieldValues<3>(
693  streachTensor, dataAtPts->logStreachTensorAtPts, MBTET));
694  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
695  rotAxis, dataAtPts->rotAxisAtPts, MBTET));
696  fe->getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
697  materialGradient, dataAtPts->GAtPts, MBTET));
698  fe->getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
699  materialGradient + "0", dataAtPts->G0AtPts, MBTET));
700  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
701  spatialDisp, dataAtPts->wAtPts, MBTET));
702 
703  // velocities
704  fe->getOpPtrVector().push_back(new OpCalculateVectorFieldValuesDot<3>(
705  spatialDisp, dataAtPts->wDotAtPts, MBTET));
706  fe->getOpPtrVector().push_back(
708  streachTensor, dataAtPts->logStreachDotTensorAtPts, MBTET));
709 
710  // calculate other derived quantities
711  fe->getOpPtrVector().push_back(
713 
714  // evaluate integration points
715  fe->getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
716  spatialDisp, tag, do_rhs, do_lhs, dataAtPts, physicalEquations));
717 
719 }
720 
722  const int tag, const bool add_elastic, const bool add_material,
723  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe_rhs,
724  boost::shared_ptr<EpElement<VolumeElementForcesAndSourcesCore>> &fe_lhs) {
726 
727  // Right hand side
728  CHKERR setBaseVolumeElementOps(tag, true, false, fe_rhs);
729 
730  // elastic
731  if (add_elastic) {
732  fe_rhs->getOpPtrVector().push_back(
734  fe_rhs->getOpPtrVector().push_back(
736  fe_rhs->getOpPtrVector().push_back(
738  fe_rhs->getOpPtrVector().push_back(
740  fe_rhs->getOpPtrVector().push_back(
742  fe_rhs->getOpPtrVector().push_back(
744  }
745 
746  // Left hand side
747  CHKERR setBaseVolumeElementOps(tag, true, true, fe_lhs);
748 
749  // elastic
750  if (add_elastic) {
751 
752  // Schur
753  fe_lhs->getOpPtrVector().push_back(
755 
756  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_du(
758  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dBubble(
760 
761  fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dP(
763  fe_lhs->getOpPtrVector().push_back(new OpSpatialEquilibrium_dw_dw(
765  fe_lhs->getOpPtrVector().push_back(
767  fe_lhs->getOpPtrVector().push_back(new OpSpatialConsistency_dBubble_domega(
768  bubbleField, rotAxis, dataAtPts, false));
769 
770  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_dP(
772  fe_lhs->getOpPtrVector().push_back(new OpSpatialPhysical_du_domega(
773  streachTensor, rotAxis, dataAtPts, false));
774 
775  fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dP(
776  rotAxis, piolaStress, dataAtPts, false));
777  fe_lhs->getOpPtrVector().push_back(new OpSpatialRotation_domega_dBubble(
778  rotAxis, bubbleField, dataAtPts, false));
779  fe_lhs->getOpPtrVector().push_back(
781  // fe_lhs->getOpPtrVector().push_back(
782  // new OpSpatialRotation_domega_du(rotAxis, streachTensor, dataAtPts));
783 
784  // Schur
785  dataAtPts->ooMatPtr = boost::make_shared<MatrixDouble>();
786  fe_lhs->getOpPtrVector().push_back(
788  if (alpha_w < std::numeric_limits<double>::epsilon()) {
789  dataAtPts->wwMatPtr = boost::make_shared<MatrixDouble>();
790  fe_lhs->getOpPtrVector().push_back(
792  } else {
793  dataAtPts->wwMatPtr.reset();
794  }
795  fe_lhs->getOpPtrVector().push_back(
797 
798  if (add_material) {
799  }
800  }
801 
803 }
804 
806  const bool add_elastic, const bool add_material,
807  boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> &fe_rhs,
808  boost::shared_ptr<EpElement<FaceElementForcesAndSourcesCore>> &fe_lhs) {
810 
811  fe_rhs =
812  boost::make_shared<EpElement<FaceElementForcesAndSourcesCore>>(mField);
813  fe_lhs =
814  boost::make_shared<EpElement<FaceElementForcesAndSourcesCore>>(mField);
815 
816  // set integration rule
817  fe_rhs->getRuleHook = FaceRule();
818  fe_lhs->getRuleHook = FaceRule();
819 
820  if (add_elastic) {
821  fe_rhs->getOpPtrVector().push_back(
823  fe_rhs->getOpPtrVector().push_back(
825 
826  }
827 
829 }
830 
834  elasticFeLhs);
837 }
838 
841  boost::shared_ptr<FEMethod> null;
842  boost::shared_ptr<EpElement<FeTractionBc>> spatial_traction_bc(
844 
845  CHKERR DMMoFEMTSSetIFunction(dm, DM_NO_ELEMENT, null, spatial_traction_bc,
846  null);
850  spatial_traction_bc);
851 
852  schurAssembly = boost::make_shared<EpFEMethod>();
858 }
859 
862  boost::shared_ptr<TsCtx> ts_ctx;
864  CHKERR TSSetIFunction(ts, f, PETSC_NULL, PETSC_NULL);
865  CHKERR TSSetIJacobian(ts, m, m, PETSC_NULL, PETSC_NULL);
866  CHKERR TSMonitorSet(ts, TsMonitorSet, ts_ctx.get(), PETSC_NULL);
867 
868  auto add_schur_streach_op = [this](auto &list, Mat S, AO aoS) {
870  for (auto &fe : list)
871  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
872  fe_cast->addStreachSchurMatrix(S, aoS);
873  else
874  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
876  };
877 
878  auto add_schur_streach_pre = [this](auto &list, Mat S, AO aoS) {
880  for (auto &fe : list)
881  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
882  fe_cast->addStreachSchurMatrix(S, aoS);
883  else
884  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
886  };
887 
888  auto add_schur_bubble_op = [this](auto &list, Mat S, AO aoS) {
890  for (auto &fe : list)
891  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
892  fe_cast->addBubbleSchurMatrix(S, aoS);
893  else
894  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
896  };
897 
898  auto add_schur_bubble_pre = [this](auto &list, Mat S, AO aoS) {
900  for (auto &fe : list)
901  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
902  fe_cast->addBubbleSchurMatrix(S, aoS);
903  else
904  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
906  };
907 
908  auto add_schur_spatial_disp_op = [this](auto &list, Mat S, AO aoS) {
910  for (auto &fe : list)
911  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
912  fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
913  else
914  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
916  };
917 
918  auto add_schur_spatial_disp_pre = [this](auto &list, Mat S, AO aoS) {
920  for (auto &fe : list)
921  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
922  fe_cast->addSpatialDispStressSchurMatrix(S, aoS);
923  else
924  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
926  };
927 
928  auto add_schur_omega_op = [this](auto &list, Mat S, AO aoS) {
930  for (auto &fe : list)
931  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.second))
932  fe_cast->addOmegaSchurMatrix(S, aoS);
933  else
934  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
936  };
937 
938  auto add_schur_omega_pre = [this](auto &list, Mat S, AO ao) {
940  for (auto &fe : list)
941  if (auto fe_cast = dynamic_cast<EpElementBase *>(fe.getSharedPtr().get()))
942  fe_cast->addOmegaSchurMatrix(S, ao);
943  else
944  SETERRQ(mField.get_comm(), MOFEM_DATA_INCONSISTENCY, "No EpElement");
946  };
947 
948  const MoFEM::Problem *schur_streach_prb_ptr;
949  CHKERR DMMoFEMGetProblemPtr(dmElasticSchurStreach, &schur_streach_prb_ptr);
950  if (auto sub_data = schur_streach_prb_ptr->subProblemData) {
951  Mat Suu;
952  CHKERR DMCreateMatrix(dmElasticSchurStreach, &Suu);
953  AO aoSuu;
954  CHKERR sub_data->getRowMap(&aoSuu);
955 
956  CHKERR add_schur_streach_op(ts_ctx->loops_to_do_IJacobian, Suu, aoSuu);
957  CHKERR add_schur_streach_pre(ts_ctx->preProcess_IJacobian, Suu, aoSuu);
958  CHKERR add_schur_streach_pre(ts_ctx->postProcess_IJacobian, Suu, aoSuu);
959 
960  CHKERR MatDestroy(&Suu);
961  CHKERR AODestroy(&aoSuu);
962 
963  const MoFEM::Problem *schur_bubble_prb_ptr;
964  CHKERR DMMoFEMGetProblemPtr(dmElasticSchurBubble, &schur_bubble_prb_ptr);
965  if (auto bubble_data = schur_bubble_prb_ptr->subProblemData) {
966  Mat SBubble;
967  CHKERR DMCreateMatrix(dmElasticSchurBubble, &SBubble);
968  AO aoSBubble;
969  CHKERR bubble_data->getRowMap(&aoSBubble);
970 
971  CHKERR add_schur_bubble_op(ts_ctx->loops_to_do_IJacobian, SBubble,
972  aoSBubble);
973  CHKERR add_schur_bubble_pre(ts_ctx->preProcess_IJacobian, SBubble,
974  aoSBubble);
975  CHKERR add_schur_bubble_pre(ts_ctx->postProcess_IJacobian, SBubble,
976  aoSBubble);
977 
978  CHKERR MatDestroy(&SBubble);
979  CHKERR AODestroy(&aoSBubble);
980 
981  const MoFEM::Problem *schur_omega_prb_ptr;
982  CHKERR DMMoFEMGetProblemPtr(dmElasticSchurOmega, &schur_omega_prb_ptr);
983  if (auto tet_stress_data = schur_omega_prb_ptr->subProblemData) {
984  Mat SOmega;
985  CHKERR DMCreateMatrix(dmElasticSchurOmega, &SOmega);
986  AO aoSOmega;
987  CHKERR tet_stress_data->getRowMap(&aoSOmega);
988 
989  CHKERR add_schur_omega_op(ts_ctx->loops_to_do_IJacobian, SOmega,
990  aoSOmega);
991  CHKERR add_schur_omega_pre(ts_ctx->preProcess_IJacobian, SOmega,
992  aoSOmega);
993  CHKERR add_schur_omega_pre(ts_ctx->postProcess_IJacobian, SOmega,
994  aoSOmega);
995 
996  const MoFEM::Problem *schur_spatial_disp_prb_ptr;
998  &schur_spatial_disp_prb_ptr);
999  if (auto spatial_disp_data =
1000  schur_spatial_disp_prb_ptr->subProblemData) {
1001 
1002  Mat Sw;
1003  CHKERR DMCreateMatrix(dmElasticSchurSpatialDisp, &Sw);
1004  AO aoSw;
1005  CHKERR spatial_disp_data->getRowMap(&aoSw);
1006 
1007  CHKERR add_schur_spatial_disp_op(ts_ctx->loops_to_do_IJacobian, Sw,
1008  aoSw);
1009  CHKERR add_schur_spatial_disp_pre(ts_ctx->preProcess_IJacobian, Sw,
1010  aoSw);
1011  CHKERR add_schur_spatial_disp_pre(ts_ctx->postProcess_IJacobian, Sw,
1012  aoSw);
1013 
1014  CHKERR MatDestroy(&Sw);
1015  CHKERR AODestroy(&aoSw);
1016  } else
1018  "Problem does not have sub-problem data");
1019 
1020  CHKERR MatDestroy(&SOmega);
1021  CHKERR AODestroy(&aoSOmega);
1022 
1023  } else
1025  "Problem does not have sub-problem data");
1026 
1027  } else
1029  "Problem does not have sub-problem data");
1030 
1031  } else
1033  "Problem does not have sub-problem data");
1034 
1035  struct Monitor : public FEMethod {
1036 
1037  using Ele = ForcesAndSourcesCore;
1040  using SetPtsData = FieldEvaluatorInterface::SetPtsData;
1041 
1042  EshelbianCore &eP;
1043  boost::shared_ptr<SetPtsData> dataFieldEval;
1044  boost::shared_ptr<VolEle> volPostProcEle;
1045  boost::shared_ptr<double> gEnergy;
1046 
1047  Monitor(EshelbianCore &ep)
1048  : eP(ep),
1049  dataFieldEval(ep.mField.getInterface<FieldEvaluatorInterface>()
1050  ->getData<VolEle>()),
1051  volPostProcEle(new VolEle(ep.mField)), gEnergy(new double) {
1052  ierr = ep.mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
1053  dataFieldEval, "EP");
1054  CHKERRABORT(PETSC_COMM_WORLD, ierr);
1055 
1056  auto no_rule = [](int, int, int) { return -1; };
1057 
1058  auto set_element_for_field_eval = [&]() {
1059  boost::shared_ptr<Ele> vol_ele(dataFieldEval->feMethodPtr.lock());
1060  vol_ele->getRuleHook = no_rule;
1061  vol_ele->getUserPolynomialBase() =
1062  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1063  vol_ele->getOpPtrVector().push_back(
1065  eP.dataAtPts->approxPAtPts));
1066  vol_ele->getOpPtrVector().push_back(
1068  eP.bubbleField, eP.dataAtPts->approxPAtPts, MBMAXTYPE));
1069  vol_ele->getOpPtrVector().push_back(
1071  eP.streachTensor, eP.dataAtPts->logStreachTensorAtPts, MBTET));
1072  vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1073  eP.rotAxis, eP.dataAtPts->rotAxisAtPts, MBTET));
1074  vol_ele->getOpPtrVector().push_back(
1076  eP.materialGradient, eP.dataAtPts->GAtPts, MBTET));
1077  vol_ele->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1078  eP.spatialDisp, eP.dataAtPts->wAtPts, MBTET));
1079  vol_ele->getOpPtrVector().push_back(
1081  eP.dataAtPts));
1082  };
1083 
1084 
1085  auto set_element_for_post_process = [&]() {
1086  volPostProcEle->getRuleHook = VolRule();
1087  volPostProcEle->getUserPolynomialBase() =
1088  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1089  volPostProcEle->getOpPtrVector().push_back(
1091  eP.dataAtPts->approxPAtPts));
1092  volPostProcEle->getOpPtrVector().push_back(
1094  eP.bubbleField, eP.dataAtPts->approxPAtPts, MBMAXTYPE));
1095  volPostProcEle->getOpPtrVector().push_back(
1097  eP.streachTensor, eP.dataAtPts->logStreachTensorAtPts, MBTET));
1098  volPostProcEle->getOpPtrVector().push_back(
1100  eP.rotAxis, eP.dataAtPts->rotAxisAtPts, MBTET));
1101  volPostProcEle->getOpPtrVector().push_back(
1103  eP.materialGradient, eP.dataAtPts->GAtPts, MBTET));
1104  volPostProcEle->getOpPtrVector().push_back(
1106  eP.dataAtPts->wAtPts, MBTET));
1107  volPostProcEle->getOpPtrVector().push_back(
1109  eP.dataAtPts));
1110  volPostProcEle->getOpPtrVector().push_back(
1111  new OpCalculateStrainEnergy(eP.spatialDisp, eP.dataAtPts, gEnergy));
1112  };
1113 
1114  set_element_for_field_eval();
1115  set_element_for_post_process();
1116  }
1117 
1118  MoFEMErrorCode preProcess() { return 0; }
1119 
1120  MoFEMErrorCode operator()() { return 0; }
1121 
1122  MoFEMErrorCode postProcess() {
1124 
1125  auto get_str_time = [](auto ts_t) {
1126  std::ostringstream ss;
1127  ss << boost::str(boost::format("%d") %
1128  static_cast<int>(std::ceil(ts_t * 1e6)));
1129  std::string s = ss.str();
1130  return s;
1131  };
1132 
1133  PetscViewer viewer;
1134  CHKERR PetscViewerBinaryOpen(
1135  PETSC_COMM_WORLD, ("restart_" + get_str_time(ts_t) + ".dat").c_str(),
1136  FILE_MODE_WRITE, &viewer);
1137  CHKERR VecView(ts_u, viewer);
1138  CHKERR PetscViewerDestroy(&viewer);
1139 
1140  CHKERR eP.postProcessResults(1, "out_sol_elastic_" + get_str_time(ts_t) +
1141  ".h5m");
1142 
1143  // Loop boundary elements with traction boundary conditions
1144  *gEnergy = 0;
1145  CHKERR eP.mField.loop_finite_elements(problemPtr->getName(), "EP",
1146  *volPostProcEle,nullptr);
1147 
1148  double body_energy;
1149  MPI_Allreduce(gEnergy.get(), &body_energy, 1, MPI_DOUBLE, MPI_SUM,
1150  eP.mField.get_comm());
1151  CHKERR PetscPrintf(eP.mField.get_comm(),
1152  "Step %d time %3.4g strain energy %3.6e\n", ts_step,
1153  ts_t, body_energy);
1154 
1155  auto post_proc_at_points = [&](std::array<double, 3> point,
1156  std::string str) {
1158 
1159  dataFieldEval->setEvalPoints(point.data(), point.size() / 3);
1160 
1161  struct OpPrint : public VolOp {
1162 
1163  EshelbianCore &eP;
1164  std::array<double, 3> point;
1165  std::string str;
1166 
1167  OpPrint(EshelbianCore &ep, std::array<double, 3> &point,
1168  std::string &str)
1169  : VolOp(ep.spatialDisp, VolOp::OPROW), eP(ep), point(point),
1170  str(str) {}
1171 
1172  MoFEMErrorCode doWork(int side, EntityType type,
1175  if (type == MBTET) {
1176  if (getGaussPts().size2()) {
1177 
1178  auto t_h = getFTensor2FromMat<3, 3>(*(eP.dataAtPts->hAtPts));
1179  auto t_approx_P =
1180  getFTensor2FromMat<3, 3>(*(eP.dataAtPts->approxPAtPts));
1181 
1185  const double jac = dEterminant(t_h);
1187  t_cauchy(i, j) = (1. / jac) * (t_approx_P(i, k) * t_h(j, k));
1188 
1189  auto add = [&]() {
1190  std::ostringstream s;
1191  s << str << " elem " << getFEEntityHandle() << " ";
1192  return s.str();
1193  };
1194 
1195  std::ostringstream print;
1196  print << add() << "comm rank " << eP.mField.get_comm_rank()
1197  << std::endl;
1198  print << add() << "point " << getVectorAdaptor(point.data(), 3)
1199  << std::endl;
1200  print << add() << "coords at gauss pts " << getCoordsAtGaussPts()
1201  << std::endl;
1202  print << add() << "w " << (*eP.dataAtPts->wAtPts) << std::endl;
1203  print << add() << "Piola " << (*eP.dataAtPts->approxPAtPts)
1204  << std::endl;
1205  print << add() << "Cauchy " << t_cauchy << std::endl;
1206  print << std::endl;
1207  CHKERR PetscSynchronizedPrintf(eP.mField.get_comm(), "%s",
1208  print.str().c_str());
1209 
1210  }
1211  }
1213  }
1214  };
1215 
1216  if (auto fe_ptr = dataFieldEval->feMethodPtr.lock()) {
1217 
1218  fe_ptr->getOpPtrVector().push_back(new OpPrint(eP, point, str));
1220  ->evalFEAtThePoint3D(point.data(), 1e-12, problemPtr->getName(),
1221  "EP", dataFieldEval,
1222  eP.mField.get_comm_rank(),
1224  fe_ptr->getOpPtrVector().pop_back();
1225  }
1226 
1228  };
1229 
1230  // Points for Cook beam
1231  std::array<double, 3> pointA = {48., 60., 5.};
1232  CHKERR post_proc_at_points(pointA, "Point A");
1233  PetscSynchronizedFlush(eP.mField.get_comm(), PETSC_STDOUT);
1234 
1235  std::array<double, 3> pointB = {48. / 2., 44. + (60. - 44.) / 2., 0.};
1236  CHKERR post_proc_at_points(pointB, "Point B");
1237  PetscSynchronizedFlush(eP.mField.get_comm(), PETSC_STDOUT);
1238 
1239  std::array<double, 3> pointC = {48. / 2., (44. - 0.) / 2., 0.};
1240  CHKERR post_proc_at_points(pointC, "Point C");
1241  PetscSynchronizedFlush(eP.mField.get_comm(), PETSC_STDOUT);
1242 
1244  }
1245  };
1246 
1247  boost::shared_ptr<FEMethod> monitor_ptr(new Monitor(*this));
1248  ts_ctx->get_loops_to_do_Monitor().push_back(
1250 
1251  CHKERR TSAppendOptionsPrefix(ts, "elastic_");
1252  CHKERR TSSetFromOptions(ts);
1254 }
1255 
1258 
1259  CHKERR TSSetDM(ts, dmElastic);
1260 
1261  SNES snes;
1262  CHKERR TSGetSNES(ts, &snes);
1263 
1264  PetscViewerAndFormat *vf;
1265  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
1266  PETSC_VIEWER_DEFAULT, &vf);
1267  CHKERR SNESMonitorSet(
1268  snes,
1269  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
1270  vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
1271 
1272  PetscSection section;
1273  CHKERR DMGetDefaultSection(dmElastic, &section);
1274  int num_fields;
1275  CHKERR PetscSectionGetNumFields(section, &num_fields);
1276  for (int ff = 0; ff != num_fields; ff++) {
1277  const char *field_name;
1278  CHKERR PetscSectionGetFieldName(section, ff, &field_name);
1279  CHKERR PetscPrintf(PETSC_COMM_WORLD, "Field %d name %s\n", ff, field_name);
1280  }
1281 
1282  CHKERR DMoFEMMeshToLocalVector(dmElastic, x, INSERT_VALUES, SCATTER_FORWARD);
1283 
1284  // PetscRandom rctx;
1285  // PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
1286  // VecSetRandom(x, rctx);
1287  // VecScale(x, -1e-1);
1288  // PetscRandomDestroy(&rctx);
1289  // CHKERR VecView(x, PETSC_VIEWER_STDOUT_WORLD);
1290 
1291  CHKERR VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD);
1292  CHKERR VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD);
1293 
1294  // CHKERR VecView(x, PETSC_VIEWER_STDOUT_WORLD);
1295 
1296  // Adding field split solver
1297 
1298  KSP ksp;
1299  CHKERR SNESGetKSP(snes, &ksp);
1300  PC pc;
1301  CHKERR KSPGetPC(ksp, &pc);
1302  PetscBool is_uu_field_split;
1303  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_uu_field_split);
1304  if (is_uu_field_split) {
1305 
1306  const MoFEM::Problem *schur_uu_ptr;
1308  if (auto uu_data = schur_uu_ptr->subProblemData) {
1309 
1310  const MoFEM::Problem *prb_ptr;
1312  map<std::string, IS> is_map;
1313  for (int ff = 0; ff != num_fields; ff++) {
1314  const char *field_name;
1315  CHKERR PetscSectionGetFieldName(section, ff, &field_name);
1316  CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1317  prb_ptr->getName(), ROW, field_name, 0, MAX_DOFS_ON_ENTITY,
1318  &is_map[field_name]);
1319  }
1320  // CHKERR mField.getInterface<ISManager>()
1321  // ->isCreateProblemFieldAndEntityType(
1322  // prb_ptr->getName(), ROW, piolaStress, MBTET, MBTET, 0,
1323  // MAX_DOFS_ON_ENTITY, &is_map["T_STRESS_ON_TETS"]);
1324  // CHKERR mField.getInterface<ISManager>()
1325  // ->isCreateProblemFieldAndEntityType(
1326  // prb_ptr->getName(), ROW, piolaStress, MBTRI, MBTRI, 0,
1327  // MAX_DOFS_ON_ENTITY, &is_map["T_STRESS_ON_TRIS"]);
1328 
1329  CHKERR uu_data->getRowIs(&is_map["E_IS_SUU"]);
1330 
1331  CHKERR PCFieldSplitSetIS(pc, NULL, is_map[streachTensor]);
1332  CHKERR PCFieldSplitSetIS(pc, NULL, is_map["E_IS_SUU"]);
1333 
1334  CHKERR PCFieldSplitSetSchurPre(pc, PC_FIELDSPLIT_SCHUR_PRE_USER,
1335  schurAssembly->Suu);
1336 
1337  CHKERR PCSetUp(pc);
1338  PetscInt n;
1339  KSP *uu_ksp;
1340  CHKERR PCFieldSplitGetSubKSP(pc, &n, &uu_ksp);
1341  PC bubble_pc;
1342  CHKERR KSPGetPC(uu_ksp[1], &bubble_pc);
1343  PetscBool is_bubble_field_split;
1344  PetscObjectTypeCompare((PetscObject)bubble_pc, PCFIELDSPLIT,
1345  &is_bubble_field_split);
1346  if (is_bubble_field_split) {
1347 
1348  const MoFEM::Problem *schur_bubble_ptr;
1349  CHKERR DMMoFEMGetProblemPtr(dmElasticSchurBubble, &schur_bubble_ptr);
1350  if (auto bubble_data = schur_bubble_ptr->subProblemData) {
1351 
1352  CHKERR bubble_data->getRowIs(&is_map["E_IS_BUBBLE"]);
1353 
1354  AO uu_ao;
1355  CHKERR uu_data->getRowMap(&uu_ao);
1356 
1357  CHKERR AOApplicationToPetscIS(uu_ao, is_map[bubbleField]);
1358  CHKERR PCFieldSplitSetIS(bubble_pc, NULL, is_map[bubbleField]);
1359  CHKERR PCFieldSplitSetIS(bubble_pc, NULL, is_map["E_IS_BUBBLE"]);
1360  CHKERR PCFieldSplitSetSchurPre(
1361  bubble_pc, PC_FIELDSPLIT_SCHUR_PRE_USER, schurAssembly->SBubble);
1362 
1363  CHKERR PCSetUp(bubble_pc);
1364  PetscInt bubble_n;
1365  KSP *bubble_ksp;
1366  CHKERR PCFieldSplitGetSubKSP(bubble_pc, &bubble_n, &bubble_ksp);
1367  PC omega_pc;
1368  CHKERR KSPGetPC(bubble_ksp[1], &omega_pc);
1369  PetscBool is_omega_field_split;
1370  PetscObjectTypeCompare((PetscObject)omega_pc, PCFIELDSPLIT,
1371  &is_omega_field_split);
1372 
1373  if (is_omega_field_split) {
1374 
1375  const MoFEM::Problem *schur_omega_ptr;
1376  CHKERR DMMoFEMGetProblemPtr(dmElasticSchurOmega, &schur_omega_ptr);
1377  if (auto omega_data = schur_omega_ptr->subProblemData) {
1378 
1379  AO bubble_ao;
1380  CHKERR bubble_data->getRowMap(&bubble_ao);
1381 
1382  CHKERR AOApplicationToPetscIS(uu_ao, is_map[rotAxis]);
1383  CHKERR AOApplicationToPetscIS(bubble_ao, is_map[rotAxis]);
1384  CHKERR omega_data->getRowIs(&is_map["E_IS_OMEGA"]);
1385 
1386  CHKERR PCFieldSplitSetIS(omega_pc, NULL, is_map[rotAxis]);
1387  CHKERR PCFieldSplitSetIS(omega_pc, NULL, is_map["E_IS_OMEGA"]);
1388  CHKERR PCFieldSplitSetSchurPre(omega_pc,
1389  PC_FIELDSPLIT_SCHUR_PRE_USER,
1390  schurAssembly->SOmega);
1391 
1392  CHKERR PCSetUp(omega_pc);
1393  PetscInt omega_n;
1394  KSP *omega_ksp;
1395  CHKERR PCFieldSplitGetSubKSP(omega_pc, &omega_n, &omega_ksp);
1396  PC w_pc;
1397  CHKERR KSPGetPC(omega_ksp[1], &w_pc);
1398  PetscBool is_w_field_split;
1399  PetscObjectTypeCompare((PetscObject)w_pc, PCFIELDSPLIT,
1400  &is_w_field_split);
1401  if (is_w_field_split) {
1402 
1403  const MoFEM::Problem *schur_w_ptr;
1405  &schur_w_ptr);
1406  if (auto w_data = schur_w_ptr->subProblemData) {
1407 
1408  AO omega_ao;
1409  CHKERR omega_data->getRowMap(&omega_ao);
1410 
1411  CHKERR AOApplicationToPetscIS(uu_ao, is_map[spatialDisp]);
1412  CHKERR AOApplicationToPetscIS(bubble_ao, is_map[spatialDisp]);
1413  CHKERR AOApplicationToPetscIS(omega_ao, is_map[spatialDisp]);
1414  CHKERR w_data->getRowIs(&is_map["E_IS_W"]);
1415 
1416  CHKERR PCFieldSplitSetIS(w_pc, NULL, is_map[spatialDisp]);
1417  CHKERR PCFieldSplitSetIS(w_pc, NULL, is_map["E_IS_W"]);
1418  CHKERR PCFieldSplitSetSchurPre(
1419  w_pc, PC_FIELDSPLIT_SCHUR_PRE_USER, schurAssembly->Sw);
1420 
1421  CHKERR AODestroy(&omega_ao);
1422  }
1423  }
1424 
1425  CHKERR PetscFree(omega_ksp);
1426  }
1427  }
1428  CHKERR PetscFree(bubble_ksp);
1429  CHKERR AODestroy(&uu_ao);
1430  }
1431  }
1432  CHKERR PetscFree(uu_ksp);
1433 
1434  for (auto &m : is_map)
1435  CHKERR ISDestroy(&m.second);
1436  }
1437  }
1438 
1439  CHKERR TSSolve(ts, x);
1441 }
1442 
1444  const std::string file) {
1446 
1447  if (!dataAtPts) {
1448  dataAtPts =
1449  boost::shared_ptr<DataAtIntegrationPts>(new DataAtIntegrationPts());
1450  dataAtPts->approxPAtPts = boost::make_shared<MatrixDouble>();
1451  dataAtPts->approxSigmaAtPts = boost::make_shared<MatrixDouble>();
1452  dataAtPts->divSigmaAtPts = boost::make_shared<MatrixDouble>();
1453  dataAtPts->wAtPts = boost::make_shared<MatrixDouble>();
1454  dataAtPts->rotAxisAtPts = boost::make_shared<MatrixDouble>();
1455  dataAtPts->logStreachTensorAtPts = boost::make_shared<MatrixDouble>();
1456  dataAtPts->streachTensorAtPts = boost::make_shared<MatrixDouble>();
1457  dataAtPts->diffStreachTensorAtPts = boost::make_shared<MatrixDouble>();
1458  dataAtPts->rotMatAtPts = boost::make_shared<MatrixDouble>();
1459  dataAtPts->diffRotMatAtPts = boost::make_shared<MatrixDouble>();
1460  dataAtPts->hAtPts = boost::make_shared<MatrixDouble>();
1461  dataAtPts->GAtPts = boost::make_shared<MatrixDouble>();
1462  dataAtPts->G0AtPts = boost::make_shared<MatrixDouble>();
1463  }
1464 
1466  post_proc.getUserPolynomialBase() =
1467  boost::shared_ptr<BaseFunction>(new CGGUserPolynomialBase());
1468 
1470  post_proc.getOpPtrVector().push_back(new OpCalculateHVecTensorField<3, 3>(
1471  piolaStress, dataAtPts->approxPAtPts));
1472  post_proc.getOpPtrVector().push_back(new OpCalculateHTensorTensorField<3, 3>(
1473  bubbleField, dataAtPts->approxPAtPts, MBMAXTYPE));
1474  post_proc.getOpPtrVector().push_back(
1476  streachTensor, dataAtPts->logStreachTensorAtPts, MBTET));
1477  post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1478  rotAxis, dataAtPts->rotAxisAtPts, MBTET));
1479  post_proc.getOpPtrVector().push_back(new OpCalculateTensor2FieldValues<3, 3>(
1480  materialGradient, dataAtPts->GAtPts, MBTET));
1481  post_proc.getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
1482  spatialDisp, dataAtPts->wAtPts, MBTET));
1483 
1484  // evaluate derived quantities
1485  post_proc.getOpPtrVector().push_back(
1487 
1488  // evaluate integration points
1489  post_proc.getOpPtrVector().push_back(physicalEquations->returnOpJacobian(
1490  spatialDisp, tag, true, false, dataAtPts, physicalEquations));
1491  post_proc.getOpPtrVector().push_back(new OpPostProcDataStructure(
1492  post_proc.postProcMesh, post_proc.mapGaussPts, spatialDisp, dataAtPts));
1493 
1494  CHKERR DMoFEMLoopFiniteElements(dM, elementVolumeName.c_str(), &post_proc);
1495  CHKERR post_proc.writeFile(file.c_str());
1497 }
1498 
1499 } // 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 mofem_fo...
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
#define DM_NO_ELEMENT
Definition: DMMoFEM.hpp:22
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
MoFEM interface unique ID.
Problem manager is used to build and partition problems.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEMErrorCode setBaseVolumeElementOps(const int tag, const bool do_rhs, const bool do_lhs, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore >> &fe)
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
field with continuous normal traction
Definition: definitions.h:177
user implemented approximation base
Definition: definitions.h:158
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
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:505
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:434
boost::shared_ptr< SubProblemData > subProblemData
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problemIt if true is assumed that matrix has the same indexing on rows and columns....
Definition: DMMMoFEM.cpp:367
MoFEMErrorCode getValueHdivForCGGBubble(MatrixDouble &pts)
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
MoFEMErrorCode addDMs(const BitRefLevel bit=BitRefLevel().set(0))
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:414
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
auto & getUserPolynomialBase()
Get the User Polynomial Base object.
Class used to pass element data to calculate base functions on tet,triangle,edge.
const FieldApproximationBase bAse
base class for all interface classes
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:481
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:548
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:190
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:512
keeps basic data about problemThis is low level structure with information about problem,...
MoFEMErrorCode removeDofsOnEntities(const std::string problem_name, const std::string field_name, const Range ents, const int lo_coeff=0, const int hi_coeff=MAX_DOFS_ON_ENTITY, int verb=VERBOSE, const bool debug=false)
Remove DOFs from problemRemove DOFs from problem which are on entities on the given range and given f...
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition: DMMMoFEM.cpp:761
Get values at integration pts for tensor filed rank 1, i.e. vector field.
static const MOFEMuuid IDD_TET_BASE_FUNCTION
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:496
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:166
MoFEMErrorCode writeFile(const std::string file_name)
wrote results in (MOAB) format, use "file_name.h5m"
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
Definition: DMMMoFEM.cpp:187
FaceElementForcesAndSourcesCoreBase VolEle
int operator()(int p_row, int p_col, int p_data) const
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
virtual int get_comm_rank() const =0
PetscErrorCode ShapeDiffMBTET(double *diffN)
calculate derivatives of shape functions
Definition: fem_tools.c:331
boost::shared_ptr< BcDispVec > bcSpatialDispVecPtr
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
Base class if inherited used to calculate base functions.
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:431
DataForcesAndSourcesCore & dAta
MoFEMErrorCode setGenericVolumeElementOps(const int tag, const bool add_elastic, const bool add_material, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore >> &fe_rhs, boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore >> &fe_lhs)
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
SmartPetscObj< DM > dmElasticSchurSpatialDisp
Sub problem of dmElastic Schur.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
boost::shared_ptr< TractionFreeBc > bcSpatialFreeTraction
TractionBc(std::string name, std::vector< double > &attr, Range &faces)
MoFEMErrorCode setUpTSElastic(TS ts, Mat m, Vec f)
MoFEMErrorCode CGG_BubbleBase_MBTET(const int p, const double *N, const double *diffN, FTensor::Tensor2< FTensor::PackPtr< double *, 9 >, 3, 3 > &phi, const int gdim)
Calculate CGGT tonsorial bubble base.
MoFEMErrorCode getValue(MatrixDouble &pts, boost::shared_ptr< BaseFunctionCtx > ctx_ptr)
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:150
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:109
FieldOrderTable & getFieldOrderTable()
Get the Field Order Table.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
boost::shared_ptr< EpElement< FaceElementForcesAndSourcesCore > > elasticBcRhs
#define CHKERR
Inline error check.
Definition: definitions.h:600
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
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:289
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
#define NBVOLUMETET_CCG_BUBBLE(P)
Bubble function for CGG H div space.
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Set integration rule on element.
PetscErrorCode DMMoFEMGetTsCtx(DM dm, MoFEM::TsCtx **ts_ctx)
get MoFEM::TsCtx data structure
Definition: DMMMoFEM.cpp:990
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > elasticFeLhs
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:302
boost::shared_ptr< TractionBcVec > bcSpatialTraction
Data on single entity (This is passed as argument to DataOperator::doWork)
structure to get information form mofem into DataForcesAndSourcesCore
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition: DMMMoFEM.cpp:349
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string &name_row)=0
set field row which finite element use
BcDisp(std::string name, std::vector< double > &attr, Range &faces)
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMMoFEM.cpp:708
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:411
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:337
virtual MPI_Comm & get_comm() const =0
continuous field
Definition: definitions.h:175
std::string getName() const
virtual const Field * get_field_structure(const std::string &name)=0
get field structure
auto createSmartDM
Creates smart DM object.
Definition: AuxPETSc.hpp:174
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:971
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
Get values at integration pts for tensor filed rank 2, i.e. matrix field.
boost::shared_ptr< SPD > getData(const double *ptr=nullptr, const int nb_eval_points=0, const double eps=1e-12, VERBOSITY_LEVELS verb=QUIET)
Get the Data object.
boost::shared_ptr< EpElement< VolumeElementForcesAndSourcesCore > > elasticFeRhs
DEPRECATED MoFEMErrorCode loop_finite_elements(const Problem *problem_ptr, const std::string &fe_name, FEMethod &method, int lower_rank, int upper_rank, MoFEMTypes bh, int verb=DEFAULT_VERBOSITY)
int operator()(int p_row, int p_col, int p_data) const
Face finite elementUser is implementing own operator at Gauss point level, by own object derived from...
field with C-1 continuity
Definition: definitions.h:178