v0.13.1
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
10
11#ifndef EXECUTABLE_DIMENSION
12#define EXECUTABLE_DIMENSION 3
13#endif
14
15#include <MoFEM.hpp>
16#include <MatrixFunction.hpp>
17
18using namespace MoFEM;
19
20template <int DIM> struct ElementsAndOps {};
21
22template <> struct ElementsAndOps<2> {
23 static constexpr FieldSpace CONTACT_SPACE = HCURL;
30};
31
32template <> struct ElementsAndOps<3> {
33 static constexpr FieldSpace CONTACT_SPACE = HDIV;
40};
41
44
45constexpr int SPACE_DIM =
46 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
47
48constexpr EntityType boundary_ent = SPACE_DIM == 3 ? MBTRI : MBEDGE;
55
63
64//! [Operators used for contact]
70 GAUSS>::OpMixDivTimesU<3, SPACE_DIM, SPACE_DIM>;
81//! [Operators used for contact]
82
83//! [Body force]
85 GAUSS>::OpSource<1, SPACE_DIM>;
86//! [Body force]
87
88//! [Only used with Hooke equation (linear material model)]
90 GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
93//! [Only used with Hooke equation (linear material model)]
94
95//! [Only used for dynamics]
100//! [Only used for dynamics]
101
102//! [Essential boundary conditions]
109//! [Essential boundary conditions]
110
111// Only used with Hencky/nonlinear material
113 GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
116
117constexpr bool is_quasi_static = true;
118constexpr bool is_large_strains = true;
119
120constexpr int order = 2;
121constexpr double young_modulus = 100;
122constexpr double poisson_ratio = 0.25;
123constexpr double rho = 0;
124constexpr double cn = 0.01;
125constexpr double spring_stiffness = 0.1;
126
127#include <ContactOps.hpp>
128#include <HenckyOps.hpp>
129#include <PostProcContact.hpp>
130using namespace ContactOps;
131using namespace HenckyOps;
132
133struct Example {
134
135 Example(MoFEM::Interface &m_field) : mField(m_field) {}
136
138
139private:
141
149
150 boost::shared_ptr<ContactOps::CommonData> commonDataPtr;
151 boost::shared_ptr<PostProcEle> postProcFe;
152 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
153 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
154 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
155
156 template <int DIM> Range getEntsOnMeshSkin();
157};
158
159//! [Run problem]
164 CHKERR bC();
165 CHKERR OPs();
166 CHKERR tsSolve();
170}
171//! [Run problem]
172
173//! [Set up problem]
177
178 // Select base
179 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
180 const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
181 PetscInt choice_base_value = AINSWORTH;
182 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
183 LASBASETOPT, &choice_base_value, PETSC_NULL);
184
186 switch (choice_base_value) {
187 case AINSWORTH:
189 MOFEM_LOG("EXAMPLE", Sev::inform)
190 << "Set AINSWORTH_LEGENDRE_BASE for displacements";
191 break;
192 case DEMKOWICZ:
194 MOFEM_LOG("EXAMPLE", Sev::inform)
195 << "Set DEMKOWICZ_JACOBI_BASE for displacents";
196 break;
197 default:
198 base = LASTBASE;
199 break;
200 }
201
202 // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
203 // Demkowicz base. We need to implement Demkowicz H1 base on tet.
204 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
205 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
206
207 CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
208 SPACE_DIM);
209 CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
210 SPACE_DIM);
211
212 CHKERR simple->setFieldOrder("U", order);
213 CHKERR simple->setFieldOrder("SIGMA", 0);
214
215 auto skin_edges = getEntsOnMeshSkin<SPACE_DIM>();
216
217 // filter not owned entities, those are not on boundary
218 Range boundary_ents;
219 ParallelComm *pcomm =
220 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
221 if (pcomm == NULL) {
222 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
223 "Communicator not created");
224 }
225
226 CHKERR pcomm->filter_pstatus(skin_edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
227 PSTATUS_NOT, -1, &boundary_ents);
228
229 CHKERR simple->setFieldOrder("SIGMA", order - 1, &boundary_ents);
230 // Range adj_edges;
231 // CHKERR mField.get_moab().get_adjacencies(boundary_ents, 1, false,
232 // adj_edges,
233 // moab::Interface::UNION);
234 // adj_edges.merge(boundary_ents);
235 // CHKERR simple->setFieldOrder("U", order, &adj_edges);
236
237 CHKERR simple->setUp();
239}
240//! [Set up problem]
241
242//! [Create common data]
245
246 auto set_matrial_stiffness = [&]() {
253 constexpr double bulk_modulus_K =
254 young_modulus / (3 * (1 - 2 * poisson_ratio));
255 constexpr double shear_modulus_G =
256 young_modulus / (2 * (1 + poisson_ratio));
257 constexpr double A =
258 (SPACE_DIM == 2) ? 2 * shear_modulus_G /
259 (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
260 : 1;
261 auto t_D =
262 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*commonDataPtr->mDPtr);
263 t_D(i, j, k, l) = 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
264 A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
265 t_kd(i, j) * t_kd(k, l);
267 };
268
269 commonDataPtr = boost::make_shared<ContactOps::CommonData>();
270
271 commonDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
272 commonDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
273 commonDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
274 commonDataPtr->contactStressPtr = boost::make_shared<MatrixDouble>();
275 commonDataPtr->contactStressDivergencePtr =
276 boost::make_shared<MatrixDouble>();
277 commonDataPtr->contactTractionPtr = boost::make_shared<MatrixDouble>();
278 commonDataPtr->contactDispPtr = boost::make_shared<MatrixDouble>();
279 commonDataPtr->curlContactStressPtr = boost::make_shared<MatrixDouble>();
280
281 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
282
283 commonDataPtr->mDPtr = boost::make_shared<MatrixDouble>();
284 commonDataPtr->mDPtr->resize(size_symm * size_symm, 1);
285
286 CHKERR set_matrial_stiffness();
288}
289//! [Create common data]
290
291//! [Boundary condition]
294 auto bc_mng = mField.getInterface<BcManager>();
296
297 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
298 "NO_CONTACT", "SIGMA", 0, 3);
299 CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(
300 simple->getProblemName(), "U");
301
303}
304//! [Boundary condition]
305
306//! [Push operators to pipeline]
310 auto bc_mng = mField.getInterface<BcManager>();
311
312 auto add_domain_base_ops = [&](auto &pipeline) {
313 auto det_ptr = boost::make_shared<VectorDouble>();
314 auto jac_ptr = boost::make_shared<MatrixDouble>();
315 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
316
317 pipeline.push_back(new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
318 pipeline.push_back(
319 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
320 pipeline.push_back(
321 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
322
323 if (SPACE_DIM == 2) {
324 pipeline.push_back(new OpMakeHdivFromHcurl());
325 pipeline.push_back(new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
326 pipeline.push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
327 pipeline.push_back(new OpSetHOWeightsOnFace());
328 } else {
329 pipeline.push_back(
330 new OpSetHOContravariantPiolaTransform(HDIV, det_ptr, jac_ptr));
331 pipeline.push_back(new OpSetHOInvJacVectorBase(HDIV, inv_jac_ptr));
332 pipeline.push_back(new OpSetHOWeights(det_ptr));
333 }
334 };
335
336 auto henky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
337 henky_common_data_ptr->matGradPtr = commonDataPtr->mGradPtr;
338 henky_common_data_ptr->matDPtr = commonDataPtr->mDPtr;
339
340 auto add_domain_ops_lhs = [&](auto &pipeline) {
341
342 if (is_large_strains) {
343 pipeline_mng->getOpDomainLhsPipeline().push_back(
345 "U", commonDataPtr->mGradPtr));
346 pipeline_mng->getOpDomainLhsPipeline().push_back(
347 new OpCalculateEigenVals<SPACE_DIM>("U", henky_common_data_ptr));
348 pipeline_mng->getOpDomainLhsPipeline().push_back(
349 new OpCalculateLogC<SPACE_DIM>("U", henky_common_data_ptr));
350 pipeline_mng->getOpDomainLhsPipeline().push_back(
351 new OpCalculateLogC_dC<SPACE_DIM>("U", henky_common_data_ptr));
352 pipeline_mng->getOpDomainLhsPipeline().push_back(
353 new OpCalculateHenckyStress<SPACE_DIM>("U", henky_common_data_ptr));
354 pipeline_mng->getOpDomainLhsPipeline().push_back(
355 new OpCalculatePiolaStress<SPACE_DIM>("U", henky_common_data_ptr));
356 pipeline_mng->getOpDomainLhsPipeline().push_back(
357 new OpHenckyTangent<SPACE_DIM>("U", henky_common_data_ptr));
358 pipeline_mng->getOpDomainLhsPipeline().push_back(
359 new OpKPiola("U", "U", henky_common_data_ptr->getMatTangent()));
360 } else {
361 pipeline.push_back(new OpKCauchy("U", "U", commonDataPtr->mDPtr));
362 }
363
364 if (!is_quasi_static) {
365 // Get pointer to U_tt shift in domain element
366 auto get_rho = [this](const double, const double, const double) {
367 auto *pipeline_mng = mField.getInterface<PipelineManager>();
368 auto &fe_domain_lhs = pipeline_mng->getDomainLhsFE();
369 return rho * fe_domain_lhs->ts_aa;
370 };
371 pipeline_mng->getOpDomainLhsPipeline().push_back(
372 new OpMass("U", "U", get_rho));
373 }
374
375 auto unity = []() { return 1; };
376 pipeline.push_back(new OpMixDivULhs("SIGMA", "U", unity, true));
377 pipeline.push_back(new OpLambdaGraULhs("SIGMA", "U", unity, true));
378
379 };
380
381 auto add_domain_ops_rhs = [&](auto &pipeline) {
382 auto get_body_force = [this](const double, const double, const double) {
383 auto *pipeline_mng = mField.getInterface<PipelineManager>();
386 auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
387 const auto time = fe_domain_rhs->ts_t;
388
389 // hardcoded gravity load
390 t_source(i) = 0;
391 t_source(1) = 1.0 * time;
392 return t_source;
393 };
394
395 pipeline.push_back(new OpBodyForce("U", get_body_force));
397 "U", commonDataPtr->mGradPtr));
398
399 if (is_large_strains) {
400 pipeline_mng->getOpDomainRhsPipeline().push_back(
401 new OpCalculateEigenVals<SPACE_DIM>("U", henky_common_data_ptr));
402 pipeline_mng->getOpDomainRhsPipeline().push_back(
403 new OpCalculateLogC<SPACE_DIM>("U", henky_common_data_ptr));
404 pipeline_mng->getOpDomainRhsPipeline().push_back(
405 new OpCalculateLogC_dC<SPACE_DIM>("U", henky_common_data_ptr));
406 pipeline_mng->getOpDomainRhsPipeline().push_back(
407 new OpCalculateHenckyStress<SPACE_DIM>("U", henky_common_data_ptr));
408 pipeline_mng->getOpDomainRhsPipeline().push_back(
409 new OpCalculatePiolaStress<SPACE_DIM>("U", henky_common_data_ptr));
410 pipeline_mng->getOpDomainRhsPipeline().push_back(new OpInternalForcePiola(
411 "U", henky_common_data_ptr->getMatFirstPiolaStress()));
412 } else {
413 pipeline.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
414 "U", commonDataPtr->mGradPtr, commonDataPtr->mStrainPtr));
416 "U", commonDataPtr->mStrainPtr, commonDataPtr->mStressPtr,
417 commonDataPtr->mDPtr));
418 pipeline.push_back(
419 new OpInternalForceCauchy("U", commonDataPtr->mStressPtr));
420 }
421
422 pipeline.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
423 "U", commonDataPtr->contactDispPtr));
424
426 "SIGMA", commonDataPtr->contactStressPtr));
427 pipeline.push_back(
429 "SIGMA", commonDataPtr->contactStressDivergencePtr));
430
431 pipeline.push_back(
432 new OpMixDivURhs("SIGMA", commonDataPtr->contactDispPtr,
433 [](double, double, double) { return 1; }));
434 pipeline.push_back(
435 new OpMixLambdaGradURhs("SIGMA", commonDataPtr->mGradPtr));
436
437 pipeline.push_back(new OpMixUTimesDivLambdaRhs(
438 "U", commonDataPtr->contactStressDivergencePtr));
439 pipeline.push_back(
440 new OpMixUTimesLambdaRhs("U", commonDataPtr->contactStressPtr));
441
442 // only in case of dynamics
443 if (!is_quasi_static) {
444 auto mat_acceleration = boost::make_shared<MatrixDouble>();
445 pipeline_mng->getOpDomainRhsPipeline().push_back(
447 mat_acceleration));
448 pipeline_mng->getOpDomainRhsPipeline().push_back(new OpInertiaForce(
449 "U", mat_acceleration, [](double, double, double) { return rho; }));
450 }
451
452 };
453
454 auto add_boundary_base_ops = [&](auto &pipeline) {
455 pipeline.push_back(new OpSetPiolaTransformOnBoundary(CONTACT_SPACE));
456 if (SPACE_DIM == 3)
457 pipeline.push_back(new OpSetHOWeightsOnFace());
458 pipeline.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
459 "U", commonDataPtr->contactDispPtr));
461 "SIGMA", commonDataPtr->contactTractionPtr));
462 };
463
464 auto add_boundary_ops_lhs = [&](auto &pipeline) {
466 auto &bc_map = bc_mng->getBcMapByBlockName();
467 for (auto bc : bc_map) {
468 if (bc_mng->checkBlock(bc, "FIX_")) {
469 MOFEM_LOG("EXAMPLE", Sev::inform)
470 << "Set boundary matrix for " << bc.first;
471 pipeline.push_back(
472 new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
473 pipeline.push_back(new OpBoundaryMass(
474 "U", "U", [](double, double, double) { return 1.; },
475 bc.second->getBcEntsPtr()));
476 }
477 }
478
479 pipeline.push_back(
480 new OpConstrainBoundaryLhs_dU("SIGMA", "U", commonDataPtr));
481 pipeline.push_back(
482 new OpConstrainBoundaryLhs_dTraction("SIGMA", "SIGMA", commonDataPtr));
483 pipeline.push_back(new OpSpringLhs(
484 "U", "U",
485
486 [this](double, double, double) { return spring_stiffness; }
487
488 ));
489
491 };
492
493 auto time_scaled = [&](double, double, double) {
494 auto *pipeline_mng = mField.getInterface<PipelineManager>();
495 auto &fe_domain_rhs = pipeline_mng->getDomainRhsFE();
496 return -fe_domain_rhs->ts_t;
497 };
498
499 auto add_boundary_ops_rhs = [&](auto &pipeline) {
501 for (auto &bc : mField.getInterface<BcManager>()->getBcMapByBlockName()) {
502 if (bc_mng->checkBlock(bc, "FIX_")) {
503 MOFEM_LOG("EXAMPLE", Sev::inform)
504 << "Set boundary residual for " << bc.first;
505 pipeline.push_back(
506 new OpSetBc("U", false, bc.second->getBcMarkersPtr()));
507 auto attr_vec = boost::make_shared<MatrixDouble>(SPACE_DIM, 1);
508 attr_vec->clear();
509 if (bc.second->bcAttributes.size() != SPACE_DIM)
510 SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
511 "Wrong size of boundary attributes vector. Set right block "
512 "size attributes. Size of attributes %d",
513 bc.second->bcAttributes.size());
514 std::copy(&bc.second->bcAttributes[0],
515 &bc.second->bcAttributes[SPACE_DIM],
516 attr_vec->data().begin());
517
518 pipeline.push_back(new OpBoundaryVec("U", attr_vec, time_scaled,
519 bc.second->getBcEntsPtr()));
520 pipeline.push_back(new OpBoundaryInternal(
521 "U", commonDataPtr->contactDispPtr,
522 [](double, double, double) { return 1.; },
523 bc.second->getBcEntsPtr()));
524
525 }
526 }
527
528 pipeline.push_back(new OpConstrainBoundaryRhs("SIGMA", commonDataPtr));
529 pipeline.push_back(new OpSpringRhs(
530 "U", commonDataPtr->contactDispPtr,
531 [this](double, double, double) { return spring_stiffness; }));
533 };
534
535 add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
536 add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
537 add_domain_ops_lhs(pipeline_mng->getOpDomainLhsPipeline());
538 add_domain_ops_rhs(pipeline_mng->getOpDomainRhsPipeline());
539
540 add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
541 add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
542 CHKERR add_boundary_ops_lhs(pipeline_mng->getOpBoundaryLhsPipeline());
543 CHKERR add_boundary_ops_rhs(pipeline_mng->getOpBoundaryRhsPipeline());
544
545 auto integration_rule_vol = [](int, int, int approx_order) {
546 return 3 * order;
547 };
548 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule_vol);
549 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule_vol);
550 auto integration_rule_boundary = [](int, int, int approx_order) {
551 return 3 * order;
552 };
553 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule_boundary);
554 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule_boundary);
555
556 // Enforce non-homegonus boundary conditions
557 auto get_bc_hook = [&]() {
559 mField, pipeline_mng->getDomainRhsFE(),
560 {boost::make_shared<TimeScale>()});
561 return hook;
562 };
563
564 pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook();
565 pipeline_mng->getDomainLhsFE()->preProcessHook = get_bc_hook();
566
568}
569//! [Push operators to pipeline]
570
571//! [Solve]
574
577 ISManager *is_manager = mField.getInterface<ISManager>();
578
579 auto set_section_monitor = [&](auto solver) {
581 SNES snes;
582 CHKERR TSGetSNES(solver, &snes);
583 PetscViewerAndFormat *vf;
584 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
585 PETSC_VIEWER_DEFAULT, &vf);
586 CHKERR SNESMonitorSet(
587 snes,
588 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
589 vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
591 };
592
593 auto scatter_create = [&](auto D, auto coeff) {
595 CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
596 ROW, "U", coeff, coeff, is);
597 int loc_size;
598 CHKERR ISGetLocalSize(is, &loc_size);
599 Vec v;
600 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
601 VecScatter scatter;
602 CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
603 return std::make_tuple(SmartPetscObj<Vec>(v),
605 };
606
607 auto set_time_monitor = [&](auto dm, auto solver) {
609 boost::shared_ptr<Monitor> monitor_ptr(
611 boost::shared_ptr<ForcesAndSourcesCore> null;
612 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
613 monitor_ptr, null, null);
615 };
616
617 auto dm = simple->getDM();
618 auto D = smartCreateDMVector(dm);
619 uXScatter = scatter_create(D, 0);
620 uYScatter = scatter_create(D, 1);
621 if (SPACE_DIM == 3)
622 uZScatter = scatter_create(D, 2);
623
624 if (is_quasi_static) {
625 auto solver = pipeline_mng->createTSIM();
626 auto D = smartCreateDMVector(dm);
627 CHKERR set_section_monitor(solver);
628 CHKERR set_time_monitor(dm, solver);
629 CHKERR TSSetSolution(solver, D);
630 CHKERR TSSetFromOptions(solver);
631 CHKERR TSSetUp(solver);
632 CHKERR TSSolve(solver, NULL);
633 } else {
634 auto solver = pipeline_mng->createTSIM2();
635 auto dm = simple->getDM();
636 auto D = smartCreateDMVector(dm);
637 auto DD = smartVectorDuplicate(D);
638 CHKERR set_section_monitor(solver);
639 CHKERR set_time_monitor(dm, solver);
640 CHKERR TS2SetSolution(solver, D, DD);
641 CHKERR TSSetFromOptions(solver);
642 CHKERR TSSetUp(solver);
643 CHKERR TSSolve(solver, NULL);
644 }
645
646 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
647 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
648 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
649
651}
652//! [Solve]
653
654//! [Postprocess results]
656//! [Postprocess results]
657
658//! [Check]
660//! [Check]
661
662template <int DIM> Range Example::getEntsOnMeshSkin() {
663 Range body_ents;
664 CHKERR mField.get_moab().get_entities_by_dimension(0, DIM, body_ents);
665 Skinner skin(&mField.get_moab());
666 Range skin_ents;
667 CHKERR skin.find_skin(0, body_ents, false, skin_ents);
668
669 return skin_ents;
670};
671
672static char help[] = "...\n\n";
673
674int main(int argc, char *argv[]) {
675
676 // Initialisation of MoFEM/PETSc and MOAB data structures
677 const char param_file[] = "param_file.petsc";
679
680 // Add logging channel for example
681 auto core_log = logging::core::get();
682 core_log->add_sink(
683 LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
684 LogManager::setLog("EXAMPLE");
685 MOFEM_LOG_TAG("EXAMPLE", "example");
686
687 try {
688
689 //! [Register MoFEM discrete manager in PETSc]
690 DMType dm_name = "DMMOFEM";
691 CHKERR DMRegister_MoFEM(dm_name);
692 //! [Register MoFEM discrete manager in PETSc
693
694 //! [Create MoAB]
695 moab::Core mb_instance; ///< mesh database
696 moab::Interface &moab = mb_instance; ///< mesh database interface
697 //! [Create MoAB]
698
699 //! [Create MoFEM]
700 MoFEM::Core core(moab); ///< finite element database
701 MoFEM::Interface &m_field = core; ///< finite element database insterface
702 //! [Create MoFEM]
703
704 //! [Load mesh]
705 Simple *simple = m_field.getInterface<Simple>();
706 CHKERR simple->getOptions();
707 CHKERR simple->loadFile("");
708 //! [Load mesh]
709
710 //! [Example]
711 Example ex(m_field);
712 CHKERR ex.runProblem();
713 //! [Example]
714 }
716
718}
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEM::FaceElementForcesAndSourcesCore FaceEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
constexpr double spring_stiffness
Definition: contact.cpp:125
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixVecTimesDivLambda< SPACE_DIM > OpMixUTimesDivLambdaRhs
Definition: contact.cpp:74
int main(int argc, char *argv[])
Definition: contact.cpp:674
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: contact.cpp:115
EntitiesFieldData::EntData EntData
Definition: contact.cpp:49
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Essential boundary conditions]
Definition: contact.cpp:113
constexpr double rho
Definition: contact.cpp:123
static char help[]
Definition: contact.cpp:672
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: contact.cpp:97
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesVec< SPACE_DIM > OpMixDivULhs
[Operators used for contact]
Definition: contact.cpp:66
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: contact.cpp:92
#define EXECUTABLE_DIMENSION
Definition: contact.cpp:12
constexpr int SPACE_DIM
Definition: contact.cpp:45
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpMixUTimesLambdaRhs
Definition: contact.cpp:76
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpSpringRhs
Definition: contact.cpp:80
ElementsAndOps< SPACE_DIM >::OpSetPiolaTransformOnBoundary OpSetPiolaTransformOnBoundary
Definition: contact.cpp:61
constexpr double poisson_ratio
Definition: contact.cpp:122
constexpr bool is_large_strains
Definition: contact.cpp:118
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: contact.cpp:106
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< SPACE_DIM > OpLambdaGraULhs
Definition: contact.cpp:68
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixDivTimesU< 3, SPACE_DIM, SPACE_DIM > OpMixDivURhs
Definition: contact.cpp:70
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Body force]
Definition: contact.cpp:90
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: contact.cpp:108
constexpr double cn
Definition: contact.cpp:124
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpMixTensorTimesGradU< SPACE_DIM > OpMixLambdaGradURhs
Definition: contact.cpp:72
constexpr int order
Definition: contact.cpp:120
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used for dynamics]
Definition: contact.cpp:104
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpSpringLhs
Definition: contact.cpp:78
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
[Operators used for contact]
Definition: contact.cpp:85
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
Definition: contact.cpp:99
constexpr bool is_quasi_static
Definition: contact.cpp:117
constexpr double young_modulus
Definition: contact.cpp:121
constexpr FieldSpace CONTACT_SPACE
Definition: contact.cpp:62
constexpr EntityType boundary_ent
Definition: contact.cpp:48
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ LASTBASE
Definition: definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
FieldSpace
approximation spaces
Definition: definitions.h:82
@ H1
continuous field
Definition: definitions.h:85
@ HCURL
field with continuous tangents
Definition: definitions.h:86
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr double shear_modulus_G
Definition: elastic.cpp:63
constexpr double bulk_modulus_K
Definition: elastic.cpp:62
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
constexpr auto t_kd
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:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
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.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
FTensor::Index< 'i', SPACE_DIM > i
const 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< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
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 >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
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:56
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
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:1003
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1086
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:1955
const double D
diffusivity
double A
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: plastic.cpp:63
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: plastic.cpp:61
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: plastic.cpp:49
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: plastic.cpp:70
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Only used with Hooke equation (linear material model)]
Definition: plastic.cpp:47
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: plastic.cpp:72
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:68
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
[Example]
Definition: plastic.cpp:120
Range getEntsOnMeshSkin()
[Check]
Definition: contact.cpp:662
boost::shared_ptr< PostProcEle > postProcFe
Definition: plastic.cpp:137
boost::shared_ptr< ContactOps::CommonData > commonDataPtr
Definition: contact.cpp:150
MoFEMErrorCode tsSolve()
FieldApproximationBase base
Definition: plot_base.cpp:68
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:143
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: contact.cpp:659
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
Definition: contact.cpp:135
MoFEMErrorCode OPs()
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
Definition: plastic.cpp:127
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:144
MoFEMErrorCode postProcess()
[Solve]
Definition: contact.cpp:655
MoFEMErrorCode setupProblem()
MoFEMErrorCode bC()
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:142
Simple interface for fast problem set-up.
Definition: BcManager.hpp:23
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:165
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Definition of the displacement bc data structure.
Definition: BCData.hpp:72
Data on single entity (This is passed as argument to DataOperator::doWork)
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
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 geometry.
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()
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryLhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Postprocess on face.
Post processing.