v0.13.2
Loading...
Searching...
No Matches
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_rhs = [&]() {
559 mField, pipeline_mng->getDomainRhsFE(),
560 {boost::make_shared<TimeScale>()});
561 return hook;
562 };
563
564 auto get_bc_hook_lhs = [&]() {
566 mField, pipeline_mng->getDomainLhsFE(),
567 {boost::make_shared<TimeScale>()});
568 return hook;
569 };
570
571 pipeline_mng->getDomainRhsFE()->preProcessHook = get_bc_hook_rhs();
572 pipeline_mng->getDomainLhsFE()->preProcessHook = get_bc_hook_lhs();
573
575}
576//! [Push operators to pipeline]
577
578//! [Solve]
581
584 ISManager *is_manager = mField.getInterface<ISManager>();
585
586 auto set_section_monitor = [&](auto solver) {
588 SNES snes;
589 CHKERR TSGetSNES(solver, &snes);
590 PetscViewerAndFormat *vf;
591 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
592 PETSC_VIEWER_DEFAULT, &vf);
593 CHKERR SNESMonitorSet(
594 snes,
595 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
596 vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
598 };
599
600 auto scatter_create = [&](auto D, auto coeff) {
602 CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
603 ROW, "U", coeff, coeff, is);
604 int loc_size;
605 CHKERR ISGetLocalSize(is, &loc_size);
606 Vec v;
607 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
608 VecScatter scatter;
609 CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
610 return std::make_tuple(SmartPetscObj<Vec>(v),
612 };
613
614 auto set_time_monitor = [&](auto dm, auto solver) {
616 boost::shared_ptr<Monitor> monitor_ptr(
618 boost::shared_ptr<ForcesAndSourcesCore> null;
619 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
620 monitor_ptr, null, null);
622 };
623
624 auto dm = simple->getDM();
625 auto D = smartCreateDMVector(dm);
626 uXScatter = scatter_create(D, 0);
627 uYScatter = scatter_create(D, 1);
628 if (SPACE_DIM == 3)
629 uZScatter = scatter_create(D, 2);
630
631 if (is_quasi_static) {
632 auto solver = pipeline_mng->createTSIM();
633 auto D = smartCreateDMVector(dm);
634 CHKERR set_section_monitor(solver);
635 CHKERR set_time_monitor(dm, solver);
636 CHKERR TSSetSolution(solver, D);
637 CHKERR TSSetFromOptions(solver);
638 CHKERR TSSetUp(solver);
639 CHKERR TSSolve(solver, NULL);
640 } else {
641 auto solver = pipeline_mng->createTSIM2();
642 auto dm = simple->getDM();
643 auto D = smartCreateDMVector(dm);
644 auto DD = smartVectorDuplicate(D);
645 CHKERR set_section_monitor(solver);
646 CHKERR set_time_monitor(dm, solver);
647 CHKERR TS2SetSolution(solver, D, DD);
648 CHKERR TSSetFromOptions(solver);
649 CHKERR TSSetUp(solver);
650 CHKERR TSSolve(solver, NULL);
651 }
652
653 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
654 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
655 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
656
658}
659//! [Solve]
660
661//! [Postprocess results]
663//! [Postprocess results]
664
665//! [Check]
667//! [Check]
668
669template <int DIM> Range Example::getEntsOnMeshSkin() {
670 Range body_ents;
671 CHKERR mField.get_moab().get_entities_by_dimension(0, DIM, body_ents);
672 Skinner skin(&mField.get_moab());
673 Range skin_ents;
674 CHKERR skin.find_skin(0, body_ents, false, skin_ents);
675
676 return skin_ents;
677};
678
679static char help[] = "...\n\n";
680
681int main(int argc, char *argv[]) {
682
683 // Initialisation of MoFEM/PETSc and MOAB data structures
684 const char param_file[] = "param_file.petsc";
686
687 // Add logging channel for example
688 auto core_log = logging::core::get();
689 core_log->add_sink(
691 LogManager::setLog("EXAMPLE");
692 MOFEM_LOG_TAG("EXAMPLE", "example");
693
694 try {
695
696 //! [Register MoFEM discrete manager in PETSc]
697 DMType dm_name = "DMMOFEM";
698 CHKERR DMRegister_MoFEM(dm_name);
699 //! [Register MoFEM discrete manager in PETSc
700
701 //! [Create MoAB]
702 moab::Core mb_instance; ///< mesh database
703 moab::Interface &moab = mb_instance; ///< mesh database interface
704 //! [Create MoAB]
705
706 //! [Create MoFEM]
707 MoFEM::Core core(moab); ///< finite element database
708 MoFEM::Interface &m_field = core; ///< finite element database insterface
709 //! [Create MoFEM]
710
711 //! [Load mesh]
712 Simple *simple = m_field.getInterface<Simple>();
713 CHKERR simple->getOptions();
714 CHKERR simple->loadFile("");
715 //! [Load mesh]
716
717 //! [Example]
718 Example ex(m_field);
719 CHKERR ex.runProblem();
720 //! [Example]
721 }
723
725}
std::string param_file
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
int main()
Definition: adol-c_atom.cpp:46
constexpr int SPACE_DIM
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
Kronecker Delta class symmetric.
constexpr double spring_stiffness
Definition: contact.cpp:125
constexpr double rho
Definition: contact.cpp:123
static char help[]
Definition: contact.cpp:679
#define EXECUTABLE_DIMENSION
Definition: contact.cpp:12
constexpr int SPACE_DIM
Definition: contact.cpp:45
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
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
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
@ 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: DMMoFEM.cpp:491
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:978
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.
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:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
double D
const double v
phase velocity of light in medium (cm/ns)
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: Common.hpp:10
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: DMMoFEM.cpp:1024
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
constexpr AssemblyType A
double young_modulus
Definition: plastic.cpp:96
double rho
Definition: plastic.cpp:98
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: plastic.cpp:66
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Only used with Hooke equation (linear material model)]
Definition: plastic.cpp:50
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< G >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:71
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: plastic.cpp:52
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< G >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBoundaryVec
Definition: plastic.cpp:73
constexpr auto size_symm
Definition: plastic.cpp:33
PetscBool is_large_strains
Definition: plastic.cpp:92
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: plastic.cpp:64
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< G >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpBoundaryInternal
Definition: plastic.cpp:75
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
[Example]
Definition: plastic.cpp:138
Range getEntsOnMeshSkin()
[Check]
Definition: contact.cpp:669
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:158
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: contact.cpp:666
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
Definition: contact.cpp:135
MoFEMErrorCode OPs()
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
Definition: plastic.cpp:145
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:159
MoFEMErrorCode postProcess()
[Solve]
Definition: contact.cpp:662
MoFEMErrorCode setupProblem()
boost::shared_ptr< PostProcEle > postProcFe
Definition: contact.cpp:151
MoFEMErrorCode bC()
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:157
Simple interface for fast problem set-up.
Definition: BcManager.hpp:24
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:200
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)
Class (Function) to enforce essential constrains.
Definition: Essential.hpp:39
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
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 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 values 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()
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:27
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Postprocess on face.
Post processing.