v0.13.0
contact.cpp
Go to the documentation of this file.
1 /**
2  * \file contact.cpp
3  * \example contact.cpp
4  *
5  * Example of contact problem
6  *
7  */
8 
9 /* This file is part of MoFEM.
10  * MoFEM is free software: you can redistribute it and/or modify it under
11  * the terms of the GNU Lesser General Public License as published by the
12  * Free Software Foundation, either version 3 of the License, or (at your
13  * option) any later version.
14  *
15  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22 
23 #ifndef EXECUTABLE_DIMENSION
24 #define EXECUTABLE_DIMENSION 3
25 #endif
26 
27 #include <MoFEM.hpp>
28 #include <MatrixFunction.hpp>
29 
30 using namespace MoFEM;
31 
32 template <int DIM> struct ElementsAndOps {};
33 
34 template <> struct ElementsAndOps<2> {
35  static constexpr FieldSpace CONTACT_SPACE = HCURL;
43 };
44 
45 template <> struct ElementsAndOps<3> {
46  static constexpr FieldSpace CONTACT_SPACE = HDIV;
54 };
55 
58 
59 constexpr int SPACE_DIM =
60  EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
61 
62 constexpr EntityType boundary_ent = SPACE_DIM == 3 ? MBTRI : MBEDGE;
69 
77 
78 //! [Operators used for contact]
84  GAUSS>::OpMixDivTimesU<3, SPACE_DIM, SPACE_DIM>;
95 //! [Operators used for contact]
96 
97 //! [Body force]
99  GAUSS>::OpSource<1, SPACE_DIM>;
100 //! [Body force]
101 
102 //! [Only used with Hooke equation (linear material model)]
104  GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
107 //! [Only used with Hooke equation (linear material model)]
108 
109 //! [Only used for dynamics]
114 //! [Only used for dynamics]
115 
116 //! [Essential boundary conditions]
123 //! [Essential boundary conditions]
124 
125 // Only used with Hencky/nonlinear material
127  GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
130 
131 constexpr bool is_quasi_static = true;
132 constexpr bool is_large_strains = true;
133 
134 constexpr int order = 2;
135 constexpr double young_modulus = 100;
136 constexpr double poisson_ratio = 0.25;
137 constexpr double rho = 0;
138 constexpr double cn = 0.01;
139 constexpr double spring_stiffness = 0.1;
140 
141 #include <OpPostProcElastic.hpp>
142 #include <ContactOps.hpp>
143 #include <HenckyOps.hpp>
144 #include <PostProcContact.hpp>
145 using namespace ContactOps;
146 using namespace HenckyOps;
147 
148 struct Example {
149 
150  Example(MoFEM::Interface &m_field) : mField(m_field) {}
151 
153 
154 private:
155  MoFEM::Interface &mField;
156 
162  MoFEMErrorCode postProcess();
164 
165  boost::shared_ptr<ContactOps::CommonData> commonDataPtr;
166  boost::shared_ptr<PostProcEle> postProcFe;
167  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
168  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
169  std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
170 
171  boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
172 
173  template <int DIM> Range getEntsOnMeshSkin();
174 };
175 
176 //! [Run problem]
179  CHKERR setupProblem();
180  CHKERR createCommonData();
181  CHKERR bC();
182  CHKERR OPs();
183  CHKERR tsSolve();
184  CHKERR postProcess();
185  CHKERR checkResults();
187 }
188 //! [Run problem]
189 
190 //! [Set up problem]
193  Simple *simple = mField.getInterface<Simple>();
194 
195  // Select base
196  enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
197  const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
198  PetscInt choice_base_value = AINSWORTH;
199  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
200  LASBASETOPT, &choice_base_value, PETSC_NULL);
201 
203  switch (choice_base_value) {
204  case AINSWORTH:
206  MOFEM_LOG("EXAMPLE", Sev::inform)
207  << "Set AINSWORTH_LEGENDRE_BASE for displacements";
208  break;
209  case DEMKOWICZ:
210  base = DEMKOWICZ_JACOBI_BASE;
211  MOFEM_LOG("EXAMPLE", Sev::inform)
212  << "Set DEMKOWICZ_JACOBI_BASE for displacents";
213  break;
214  default:
215  base = LASTBASE;
216  break;
217  }
218 
219  // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
220  // Demkowicz base. We need to implement Demkowicz H1 base on tet.
221  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
222  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
223 
224  CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
225  SPACE_DIM);
226  CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
227  SPACE_DIM);
228 
229  CHKERR simple->setFieldOrder("U", order);
230  CHKERR simple->setFieldOrder("SIGMA", 0);
231 
232  auto skin_edges = getEntsOnMeshSkin<SPACE_DIM>();
233 
234  // filter not owned entities, those are not on boundary
235  Range boundary_ents;
236  ParallelComm *pcomm =
237  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
238  if (pcomm == NULL) {
239  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
240  "Communicator not created");
241  }
242 
243  CHKERR pcomm->filter_pstatus(skin_edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
244  PSTATUS_NOT, -1, &boundary_ents);
245 
246  CHKERR simple->setFieldOrder("SIGMA", order - 1, &boundary_ents);
247  // Range adj_edges;
248  // CHKERR mField.get_moab().get_adjacencies(boundary_ents, 1, false,
249  // adj_edges,
250  // moab::Interface::UNION);
251  // adj_edges.merge(boundary_ents);
252  // CHKERR simple->setFieldOrder("U", order, &adj_edges);
253 
254  CHKERR simple->setUp();
256 }
257 //! [Set up problem]
258 
259 //! [Create common data]
262 
263  auto set_matrial_stiffness = [&]() {
270  constexpr double bulk_modulus_K =
271  young_modulus / (3 * (1 - 2 * poisson_ratio));
272  constexpr double shear_modulus_G =
273  young_modulus / (2 * (1 + poisson_ratio));
274  constexpr double A =
275  (SPACE_DIM == 2) ? 2 * shear_modulus_G /
276  (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
277  : 1;
278  auto t_D =
279  getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*commonDataPtr->mDPtr);
280  t_D(i, j, k, l) = 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
281  A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
282  t_kd(i, j) * t_kd(k, l);
284  };
285 
286  commonDataPtr = boost::make_shared<ContactOps::CommonData>();
287 
288  commonDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
289  commonDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
290  commonDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
291  commonDataPtr->contactStressPtr = boost::make_shared<MatrixDouble>();
292  commonDataPtr->contactStressDivergencePtr =
293  boost::make_shared<MatrixDouble>();
294  commonDataPtr->contactTractionPtr = boost::make_shared<MatrixDouble>();
295  commonDataPtr->contactDispPtr = boost::make_shared<MatrixDouble>();
296  commonDataPtr->curlContactStressPtr = boost::make_shared<MatrixDouble>();
297 
298  constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
299 
300  commonDataPtr->mDPtr = boost::make_shared<MatrixDouble>();
301  commonDataPtr->mDPtr->resize(size_symm * size_symm, 1);
302 
303  CHKERR set_matrial_stiffness();
305 }
306 //! [Create common data]
307 
308 //! [Boundary condition]
311  auto bc_mng = mField.getInterface<BcManager>();
312  auto simple = mField.getInterface<Simple>();
313 
314  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
315  "NO_CONTACT", "SIGMA", 0, 3);
316 
317  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
318  "U", 0, 0);
319  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
320  "U", 1, 1);
321  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
322  "U", 2, 2);
323  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
324  "REMOVE_ALL", "U", 0, 3);
325 
326  CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_X", "U",
327  0, 0);
328  CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_Y", "U",
329  1, 1);
330  CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_Z", "U",
331  2, 2);
332  CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
333  "U", 0, 3);
334 
335  boundaryMarker = bc_mng->getMergedBlocksMarker(vector<string>{"FIX_"});
336 
338 }
339 //! [Boundary condition]
340 
341 //! [Push operators to pipeline]
344  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
345  auto bc_mng = mField.getInterface<BcManager>();
346 
347  auto add_domain_base_ops = [&](auto &pipeline) {
348  auto det_ptr = boost::make_shared<VectorDouble>();
349  auto jac_ptr = boost::make_shared<MatrixDouble>();
350  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
351  if (SPACE_DIM == 2) {
352  pipeline.push_back(new OpCalculateHOJacForFace(jac_ptr));
353  pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
354  pipeline.push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
355  pipeline.push_back(new OpMakeHdivFromHcurl());
356  pipeline.push_back(new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
357  pipeline.push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
358  pipeline.push_back(new OpSetHOWeightsOnFace());
359  } else {
360  pipeline.push_back(new OpCalculateHOJacVolume(jac_ptr));
361  pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
362  pipeline.push_back(
363  new OpSetHOContravariantPiolaTransform(HDIV, det_ptr, jac_ptr));
364  pipeline.push_back(new OpSetHOInvJacVectorBase(HDIV, inv_jac_ptr));
365  pipeline.push_back(new OpSetHOInvJacToScalarBases(H1, inv_jac_ptr));
366  pipeline.push_back(new OpSetHOWeights(det_ptr));
367  }
368  };
369 
370  auto henky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
371  henky_common_data_ptr->matGradPtr = commonDataPtr->mGradPtr;
372  henky_common_data_ptr->matDPtr = commonDataPtr->mDPtr;
373 
374  auto add_domain_ops_lhs = [&](auto &pipeline) {
375  pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
376 
377  if (is_large_strains) {
378  pipeline_mng->getOpDomainLhsPipeline().push_back(
380  "U", commonDataPtr->mGradPtr));
381  pipeline_mng->getOpDomainLhsPipeline().push_back(
382  new OpCalculateEigenVals<SPACE_DIM>("U", henky_common_data_ptr));
383  pipeline_mng->getOpDomainLhsPipeline().push_back(
384  new OpCalculateLogC<SPACE_DIM>("U", henky_common_data_ptr));
385  pipeline_mng->getOpDomainLhsPipeline().push_back(
386  new OpCalculateLogC_dC<SPACE_DIM>("U", henky_common_data_ptr));
387  pipeline_mng->getOpDomainLhsPipeline().push_back(
388  new OpCalculateHenckyStress<SPACE_DIM>("U", henky_common_data_ptr));
389  pipeline_mng->getOpDomainLhsPipeline().push_back(
390  new OpCalculatePiolaStress<SPACE_DIM>("U", henky_common_data_ptr));
391  pipeline_mng->getOpDomainLhsPipeline().push_back(
392  new OpHenckyTangent<SPACE_DIM>("U", henky_common_data_ptr));
393  pipeline_mng->getOpDomainLhsPipeline().push_back(
394  new OpKPiola("U", "U", henky_common_data_ptr->getMatTangent()));
395  } else {
396  pipeline.push_back(new OpKCauchy("U", "U", commonDataPtr->mDPtr));
397  }
398 
399  if (!is_quasi_static) {
400  // Get pointer to U_tt shift in domain element
401  auto get_rho = [this](const double, const double, const double) {
402  auto *pipeline_mng = mField.getInterface<PipelineManager>();
403  auto &fe_domain_lhs = pipeline_mng->getDomainLhsFE();
404  return rho * fe_domain_lhs->ts_aa;
405  };
406  pipeline_mng->getOpDomainLhsPipeline().push_back(
407  new OpMass("U", "U", get_rho));
408  }
409 
410  auto unity = []() { return 1; };
411  pipeline.push_back(new OpMixDivULhs("SIGMA", "U", unity, true));
412  pipeline.push_back(new OpLambdaGraULhs("SIGMA", "U", unity, true));
413 
414  pipeline.push_back(new OpUnSetBc("U"));
415  };
416 
417  auto add_domain_ops_rhs = [&](auto &pipeline) {
418  auto get_body_force = [this](const double, const double, const double) {
419  auto *pipeline_mng = mField.getInterface<PipelineManager>();
422  auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
423  const auto time = fe_domain_rhs->ts_t;
424 
425  // hardcoded gravity load
426  t_source(i) = 0;
427  t_source(1) = 1.0 * time;
428  return t_source;
429  };
430 
431  pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
432 
433  pipeline.push_back(new OpBodyForce("U", get_body_force));
435  "U", commonDataPtr->mGradPtr));
436 
437  if (is_large_strains) {
438  pipeline_mng->getOpDomainRhsPipeline().push_back(
439  new OpCalculateEigenVals<SPACE_DIM>("U", henky_common_data_ptr));
440  pipeline_mng->getOpDomainRhsPipeline().push_back(
441  new OpCalculateLogC<SPACE_DIM>("U", henky_common_data_ptr));
442  pipeline_mng->getOpDomainRhsPipeline().push_back(
443  new OpCalculateLogC_dC<SPACE_DIM>("U", henky_common_data_ptr));
444  pipeline_mng->getOpDomainRhsPipeline().push_back(
445  new OpCalculateHenckyStress<SPACE_DIM>("U", henky_common_data_ptr));
446  pipeline_mng->getOpDomainRhsPipeline().push_back(
447  new OpCalculatePiolaStress<SPACE_DIM>("U", henky_common_data_ptr));
448  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpInternalForcePiola(
449  "U", henky_common_data_ptr->getMatFirstPiolaStress()));
450  } else {
451  pipeline.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
452  "U", commonDataPtr->mGradPtr, commonDataPtr->mStrainPtr));
454  "U", commonDataPtr->mStrainPtr, commonDataPtr->mStressPtr,
455  commonDataPtr->mDPtr));
456  pipeline.push_back(
457  new OpInternalForceCauchy("U", commonDataPtr->mStressPtr));
458  }
459 
460  pipeline.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
461  "U", commonDataPtr->contactDispPtr));
462 
464  "SIGMA", commonDataPtr->contactStressPtr));
465  pipeline.push_back(
467  "SIGMA", commonDataPtr->contactStressDivergencePtr));
468 
469  pipeline.push_back(
470  new OpMixDivURhs("SIGMA", commonDataPtr->contactDispPtr,
471  [](double, double, double) { return 1; }));
472  pipeline.push_back(
473  new OpMixLambdaGradURhs("SIGMA", commonDataPtr->mGradPtr));
474 
475  pipeline.push_back(new OpMixUTimesDivLambdaRhs(
476  "U", commonDataPtr->contactStressDivergencePtr));
477  pipeline.push_back(
478  new OpMixUTimesLambdaRhs("U", commonDataPtr->contactStressPtr));
479 
480  // only in case of dynamics
481  if (!is_quasi_static) {
482  auto mat_acceleration = boost::make_shared<MatrixDouble>();
483  pipeline_mng->getOpDomainRhsPipeline().push_back(
485  mat_acceleration));
486  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpInertiaForce(
487  "U", mat_acceleration, [](double, double, double) { return rho; }));
488  }
489 
490  pipeline.push_back(new OpUnSetBc("U"));
491  };
492 
493  auto add_boundary_base_ops = [&](auto &pipeline) {
494  pipeline.push_back(new OpSetPiolaTransformOnBoundary(CONTACT_SPACE));
495  if (SPACE_DIM == 3)
496  pipeline.push_back(new OpSetHOWeightsOnFace());
497  pipeline.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
498  "U", commonDataPtr->contactDispPtr));
500  "SIGMA", commonDataPtr->contactTractionPtr));
501  };
502 
503  auto add_boundary_ops_lhs = [&](auto &pipeline) {
505  auto &bc_map = bc_mng->getBcMapByBlockName();
506  for (auto bc : bc_map) {
507  if (bc_mng->checkBlock(bc, "FIX_")) {
508  MOFEM_LOG("EXAMPLE", Sev::inform)
509  << "Set boundary matrix for " << bc.first;
510  pipeline.push_back(
511  new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
512  pipeline.push_back(new OpBoundaryMass(
513  "U", "U", [](double, double, double) { return 1.; },
514  bc.second->getBcEntsPtr()));
515  pipeline.push_back(new OpUnSetBc("U"));
516  }
517  }
518 
519  pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
520  pipeline.push_back(
521  new OpConstrainBoundaryLhs_dU("SIGMA", "U", commonDataPtr));
522  pipeline.push_back(
523  new OpConstrainBoundaryLhs_dTraction("SIGMA", "SIGMA", commonDataPtr));
524  pipeline.push_back(new OpSpringLhs(
525  "U", "U",
526 
527  [this](double, double, double) { return spring_stiffness; }
528 
529  ));
530  if (boundaryMarker)
531  pipeline.push_back(new OpUnSetBc("U"));
533  };
534 
535  auto time_scaled = [&](double, double, double) {
536  auto *pipeline_mng = mField.getInterface<PipelineManager>();
537  auto &fe_domain_rhs = pipeline_mng->getDomainRhsFE();
538  return -fe_domain_rhs->ts_t;
539  };
540 
541  auto add_boundary_ops_rhs = [&](auto &pipeline) {
543  for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName()) {
544  if (bc_mng->checkBlock(bc, "FIX_")) {
545  MOFEM_LOG("EXAMPLE", Sev::inform)
546  << "Set boundary residual for " << bc.first;
547  pipeline.push_back(
548  new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
549  auto attr_vec = boost::make_shared<MatrixDouble>(SPACE_DIM, 1);
550  attr_vec->clear();
551  if (bc.second->bcAttributes.size() != SPACE_DIM)
552  SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
553  "Wrong size of boundary attributes vector. Set right block "
554  "size attributes. Size of attributes %d",
555  bc.second->bcAttributes.size());
556  std::copy(&bc.second->bcAttributes[0],
557  &bc.second->bcAttributes[SPACE_DIM],
558  attr_vec->data().begin());
559 
560  pipeline.push_back(new OpBoundaryVec("U", attr_vec, time_scaled,
561  bc.second->getBcEntsPtr()));
562  pipeline.push_back(new OpBoundaryInternal(
563  "U", commonDataPtr->contactDispPtr,
564  [](double, double, double) { return 1.; },
565  bc.second->getBcEntsPtr()));
566 
567  pipeline.push_back(new OpUnSetBc("U"));
568  }
569  }
570 
571  pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
572  pipeline.push_back(new OpConstrainBoundaryRhs("SIGMA", commonDataPtr));
573  pipeline.push_back(new OpSpringRhs(
574  "U", commonDataPtr->contactDispPtr,
575  [this](double, double, double) { return spring_stiffness; }));
576  pipeline.push_back(new OpUnSetBc("U"));
578  };
579 
580  add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
581  add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
582  add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline());
583  add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline());
584 
585  add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
586  add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
587  CHKERR add_boundary_ops_lhs(pipeline_mng->getOpBoundaryLhsPipeline());
588  CHKERR add_boundary_ops_rhs(pipeline_mng->getOpBoundaryRhsPipeline());
589 
590  auto integration_rule_vol = [](int, int, int approx_order) {
591  return 3 * order;
592  };
593  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule_vol);
594  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule_vol);
595  auto integration_rule_boundary = [](int, int, int approx_order) {
596  return 3 * order;
597  };
598  CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule_boundary);
599  CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule_boundary);
600 
602 }
603 //! [Push operators to pipeline]
604 
605 //! [Solve]
608 
609  Simple *simple = mField.getInterface<Simple>();
610  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
611  ISManager *is_manager = mField.getInterface<ISManager>();
612 
613  auto set_section_monitor = [&](auto solver) {
615  SNES snes;
616  CHKERR TSGetSNES(solver, &snes);
617  PetscViewerAndFormat *vf;
618  CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
619  PETSC_VIEWER_DEFAULT, &vf);
620  CHKERR SNESMonitorSet(
621  snes,
622  (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
623  vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
625  };
626 
627  auto scatter_create = [&](auto D, auto coeff) {
629  CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
630  ROW, "U", coeff, coeff, is);
631  int loc_size;
632  CHKERR ISGetLocalSize(is, &loc_size);
633  Vec v;
634  CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
635  VecScatter scatter;
636  CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
637  return std::make_tuple(SmartPetscObj<Vec>(v),
638  SmartPetscObj<VecScatter>(scatter));
639  };
640 
641  auto set_time_monitor = [&](auto dm, auto solver) {
643  boost::shared_ptr<Monitor> monitor_ptr(
644  new Monitor(dm, commonDataPtr, uXScatter, uYScatter, uZScatter));
645  boost::shared_ptr<ForcesAndSourcesCore> null;
646  CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
647  monitor_ptr, null, null);
649  };
650 
651  auto dm = simple->getDM();
652  auto D = smartCreateDMVector(dm);
653  uXScatter = scatter_create(D, 0);
654  uYScatter = scatter_create(D, 1);
655  if (SPACE_DIM == 3)
656  uZScatter = scatter_create(D, 2);
657 
658  if (is_quasi_static) {
659  auto solver = pipeline_mng->createTSIM();
660  auto D = smartCreateDMVector(dm);
661  CHKERR set_section_monitor(solver);
662  CHKERR set_time_monitor(dm, solver);
663  CHKERR TSSetSolution(solver, D);
664  CHKERR TSSetFromOptions(solver);
665  CHKERR TSSetUp(solver);
666  CHKERR TSSolve(solver, NULL);
667  } else {
668  auto solver = pipeline_mng->createTSIM2();
669  auto dm = simple->getDM();
670  auto D = smartCreateDMVector(dm);
671  auto DD = smartVectorDuplicate(D);
672  CHKERR set_section_monitor(solver);
673  CHKERR set_time_monitor(dm, solver);
674  CHKERR TS2SetSolution(solver, D, DD);
675  CHKERR TSSetFromOptions(solver);
676  CHKERR TSSetUp(solver);
677  CHKERR TSSolve(solver, NULL);
678  }
679 
680  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
681  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
682  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
683 
685 }
686 //! [Solve]
687 
688 //! [Postprocess results]
690 //! [Postprocess results]
691 
692 //! [Check]
694 //! [Check]
695 
696 template <int DIM> Range Example::getEntsOnMeshSkin() {
697  Range body_ents;
698  CHKERR mField.get_moab().get_entities_by_dimension(0, DIM, body_ents);
699  Skinner skin(&mField.get_moab());
700  Range skin_ents;
701  CHKERR skin.find_skin(0, body_ents, false, skin_ents);
702 
703  return skin_ents;
704 };
705 
706 static char help[] = "...\n\n";
707 
708 int main(int argc, char *argv[]) {
709 
710  // Initialisation of MoFEM/PETSc and MOAB data structures
711  const char param_file[] = "param_file.petsc";
712  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
713 
714  // Add logging channel for example
715  auto core_log = logging::core::get();
716  core_log->add_sink(
718  LogManager::setLog("EXAMPLE");
719  MOFEM_LOG_TAG("EXAMPLE", "example");
720 
721  try {
722 
723  //! [Register MoFEM discrete manager in PETSc]
724  DMType dm_name = "DMMOFEM";
725  CHKERR DMRegister_MoFEM(dm_name);
726  //! [Register MoFEM discrete manager in PETSc
727 
728  //! [Create MoAB]
729  moab::Core mb_instance; ///< mesh database
730  moab::Interface &moab = mb_instance; ///< mesh database interface
731  //! [Create MoAB]
732 
733  //! [Create MoFEM]
734  MoFEM::Core core(moab); ///< finite element database
735  MoFEM::Interface &m_field = core; ///< finite element database insterface
736  //! [Create MoFEM]
737 
738  //! [Load mesh]
739  Simple *simple = m_field.getInterface<Simple>();
740  CHKERR simple->getOptions();
741  CHKERR simple->loadFile("");
742  //! [Load mesh]
743 
744  //! [Example]
745  Example ex(m_field);
746  CHKERR ex.runProblem();
747  //! [Example]
748  }
749  CATCH_ERRORS;
750 
752 }
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class symmetric.
constexpr double spring_stiffness
Definition: contact.cpp:139
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixVecTimesDivLambda< SPACE_DIM > OpMixUTimesDivLambdaRhs
Definition: contact.cpp:88
int main(int argc, char *argv[])
Definition: contact.cpp:708
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: contact.cpp:129
EntitiesFieldData::EntData EntData
Definition: contact.cpp:63
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Essential boundary conditions]
Definition: contact.cpp:127
constexpr double rho
Definition: contact.cpp:137
static char help[]
Definition: contact.cpp:706
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: contact.cpp:111
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesVec< SPACE_DIM > OpMixDivULhs
[Operators used for contact]
Definition: contact.cpp:80
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: contact.cpp:106
#define EXECUTABLE_DIMENSION
Definition: contact.cpp:24
constexpr int SPACE_DIM
Definition: contact.cpp:59
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpMixUTimesLambdaRhs
Definition: contact.cpp:90
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpSpringRhs
Definition: contact.cpp:94
ElementsAndOps< SPACE_DIM >::OpSetPiolaTransformOnBoundary OpSetPiolaTransformOnBoundary
Definition: contact.cpp:75
constexpr double poisson_ratio
Definition: contact.cpp:136
constexpr bool is_large_strains
Definition: contact.cpp:132
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: contact.cpp:120
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< SPACE_DIM > OpLambdaGraULhs
Definition: contact.cpp:82
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixDivTimesU< 3, SPACE_DIM, SPACE_DIM > OpMixDivURhs
Definition: contact.cpp:84
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Body force]
Definition: contact.cpp:104
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: contact.cpp:122
constexpr double cn
Definition: contact.cpp:138
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixTensorTimesGradU< SPACE_DIM > OpMixLambdaGradURhs
Definition: contact.cpp:86
constexpr int order
Definition: contact.cpp:134
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used for dynamics]
Definition: contact.cpp:118
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpSpringLhs
Definition: contact.cpp:92
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
[Operators used for contact]
Definition: contact.cpp:99
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
Definition: contact.cpp:113
constexpr bool is_quasi_static
Definition: contact.cpp:131
constexpr double young_modulus
Definition: contact.cpp:135
constexpr FieldSpace CONTACT_SPACE
Definition: contact.cpp:76
constexpr EntityType boundary_ent
Definition: contact.cpp:62
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
constexpr double shear_modulus_G
constexpr double bulk_modulus_K
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
@ ROW
Definition: definitions.h:136
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
FieldApproximationBase
approximation base
Definition: definitions.h:71
@ LASTBASE
Definition: definitions.h:82
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
FieldSpace
approximation spaces
Definition: definitions.h:95
@ H1
continuous field
Definition: definitions.h:98
@ HCURL
field with continuous tangents
Definition: definitions.h:99
@ HDIV
field with continuous normal traction
Definition: definitions.h:100
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
constexpr auto t_kd
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
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:481
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
FTensor::Index< 'i', SPACE_DIM > i
double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpSpringLhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 3, 3 > OpMixUTimesLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixDivTimesVec< 3 > OpMixDivULhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixVecTimesDivLambda< 3 > OpMixUTimesDivLambdaRhs
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpSpringRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 3, 3 > OpMixDivURhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< 3 > OpLambdaGraULhs
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, 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 Monitor To TS solver.
Definition: DMMMoFEM.cpp:994
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
CoreTmp< 0 > Core
Definition: Core.hpp:1096
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
const double D
diffusivity
double A
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: plastic.cpp:88
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: plastic.cpp:86
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: plastic.cpp:74
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: plastic.cpp:95
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Body force]
Definition: plastic.cpp:72
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: plastic.cpp:97
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:93
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:41
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
[Postprocessing]
PipelineManager::FaceEle DomainEle
[Example]
Range getEntsOnMeshSkin()
[Check]
Definition: contact.cpp:696
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
Definition: contact.cpp:165
MoFEMErrorCode tsSolve()
MoFEMErrorCode checkResults()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
Definition: contact.cpp:150
MoFEMErrorCode OPs()
MoFEMErrorCode runProblem()
MoFEMErrorCode postProcess()
[Solve]
Definition: contact.cpp:689
MoFEMErrorCode setupProblem()
MoFEMErrorCode bC()
Simple interface for fast problem set-up.
Definition: BcManager.hpp:27
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:131
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:33
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
Calculate jacobian on Hex or other volume which is not simplex.
Calculate divergence of tonsorial field using vectorial base.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Calculate trace of vector (Hdiv/Hcurl) space.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Approximate field valuse for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
transform Hdiv base fluxes from reference element to physical triangle
Make Hdiv space from Hcurl space in 2d.
Set indices on entities on finite element.
Apply contravariant (Piola) transfer to Hdiv space for HO geometr.
Set inverse jacobian to base functions.
transform local reference derivatives of shape function to global derivatives if higher order geometr...
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryLhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Simple interface for fast problem set-up.
Definition: Simple.hpp:33
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
Postprocess on face.
Post processing.