v0.13.2
Loading...
Searching...
No Matches
contact.cpp
Go to the documentation of this file.
1/**
2 * \file contact.cpp
3 * \CONTACT contact.cpp
4 *
5 * CONTACT of contact problem
6 *
7 * @copyright Copyright (c) 2023
8 */
9
10#ifndef EXECUTABLE_DIMENSION
11#define EXECUTABLE_DIMENSION 3
12#endif
13
14#include <MoFEM.hpp>
15#include <MatrixFunction.hpp>
16
17#ifdef PYTHON_SFD
18#include <boost/python.hpp>
19#include <boost/python/def.hpp>
20namespace bp = boost::python;
21#endif
22
23using namespace MoFEM;
24
25constexpr AssemblyType A = AssemblyType::PETSC; //< selected assembly type
26constexpr IntegrationType I =
27 IntegrationType::GAUSS; //< selected integration type
28
29template <int DIM> struct ElementsAndOps {};
30
32 static constexpr FieldSpace CONTACT_SPACE = HCURL;
33};
34
36 static constexpr FieldSpace CONTACT_SPACE = HDIV;
37};
38
41
42constexpr int SPACE_DIM =
43 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
44
50
55
56//! [Operators used for contact]
62 GAUSS>::OpMixDivTimesU<3, SPACE_DIM, SPACE_DIM>;
73//! [Operators used for contact]
74
75//! [Only used for dynamics]
80//! [Only used for dynamics]
81
82//! [Essential boundary conditions]
87//! [Essential boundary conditions]
88
89// Only used with Hencky/nonlinear material
91 GAUSS>::OpGradTensorGrad<1, SPACE_DIM, SPACE_DIM, 1>;
94
95namespace ContactOps {
96
97struct DomainBCs {};
98struct BoundaryBCs {};
99
106
108 MoFEM::Interface &m_field,
109 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
110 std::string field_name, std::string block_name,
111 boost::shared_ptr<MatrixDouble> mat_D_Ptr, Sev sev);
112
113}; // namespace ContactOps
114
115constexpr bool is_quasi_static = true;
116
117int order = 2;
118constexpr int geom_order =
119 1; ///< Currently calculation of normals at integration points is missing on
120 ///< edges (i.e. 2d case). We have to restrict to linear geometry in 2d.
121double young_modulus = 100;
122double poisson_ratio = 0.25;
123double rho = 0.0;
124double cn = 0.1;
125double spring_stiffness = 0.1;
126
127double alpha_dumping = 0;
128
129#include <HenckyOps.hpp>
130using namespace HenckyOps;
131#include <ContactOps.hpp>
132#include <PostProcContact.hpp>
135
136using namespace ContactOps;
137
138struct Contact {
139
140 Contact(MoFEM::Interface &m_field) : mField(m_field) {}
141
143
144private:
146
154
155 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uXScatter;
156 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uYScatter;
157 std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter>> uZScatter;
158 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
159
160#ifdef PYTHON_SFD
161 boost::shared_ptr<SDFPython> sdfPythonPtr;
162#endif
163};
164
165//! [Run problem]
170 CHKERR bC();
171 CHKERR OPs();
172 CHKERR tsSolve();
176}
177//! [Run problem]
178
179//! [Set up problem]
183
184 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
185 // CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,
186 // PETSC_NULL);
187 MOFEM_LOG("CONTACT", Sev::inform) << "Order " << order;
188 MOFEM_LOG("CONTACT", Sev::inform) << "Geom order " << geom_order;
189
190 // Select base
191 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
192 const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
193 PetscInt choice_base_value = AINSWORTH;
194 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
195 LASBASETOPT, &choice_base_value, PETSC_NULL);
196
198 switch (choice_base_value) {
199 case AINSWORTH:
201 MOFEM_LOG("CONTACT", Sev::inform)
202 << "Set AINSWORTH_LEGENDRE_BASE for displacements";
203 break;
204 case DEMKOWICZ:
206 MOFEM_LOG("CONTACT", Sev::inform)
207 << "Set DEMKOWICZ_JACOBI_BASE for displacements";
208 break;
209 default:
210 base = LASTBASE;
211 break;
212 }
213
214 // Note: For tets we have only H1 Ainsworth base, for Hex we have only H1
215 // Demkowicz base. We need to implement Demkowicz H1 base on tet.
216 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
217 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
218 CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
219 SPACE_DIM);
220 CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
221 SPACE_DIM);
222 // CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
223
224
225 CHKERR simple->setFieldOrder("U", order);
226 // CHKERR simple->setFieldOrder("SIGMA", 0);
227 // CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
228
229 auto get_skin = [&]() {
230 Range body_ents;
231 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
232 Skinner skin(&mField.get_moab());
233 Range skin_ents;
234 CHKERR skin.find_skin(0, body_ents, false, skin_ents);
235 return skin_ents;
236 };
237
238 auto filter_blocks = [&](auto skin) {
239 Range contact_range;
240 for (auto m :
242
243 (boost::format("%s(.*)") % "CONTACT").str()
244
245 ))
246
247 ) {
248 MOFEM_LOG("CONTACT", Sev::inform)
249 << "Find contact block set: " << m->getName();
250 auto meshset = m->getMeshset();
251 CHKERR mField.get_moab().get_entities_by_dimension(meshset, SPACE_DIM - 1,
252 contact_range, true);
253
254 MOFEM_LOG("SYNC", Sev::inform)
255 << "Nb entities in contact surface: " << contact_range.size();
257 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
258 contact_range);
259 skin = intersect(skin, contact_range);
260 }
261 return skin;
262 };
263
264 auto filter_true_skin = [&](auto skin) {
265 Range boundary_ents;
266 ParallelComm *pcomm =
267 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
268 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
269 PSTATUS_NOT, -1, &boundary_ents);
270 return boundary_ents;
271 };
272
273 auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
274 CHKERR simple->setFieldOrder("SIGMA", 0);
275 CHKERR simple->setFieldOrder("SIGMA", order - 1, &boundary_ents);
276
277 CHKERR simple->setUp();
278
279 // auto project_ho_geometry = [&]() {
280 // Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
281 // return mField.loop_dofs("GEOMETRY", ent_method);
282 // };
283 // CHKERR project_ho_geometry();
284
286}
287//! [Set up problem]
288
289//! [Create common data]
292
293 auto get_options = [&]() {
295 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
296 &young_modulus, PETSC_NULL);
297 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
298 &poisson_ratio, PETSC_NULL);
299 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
300 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn", &cn, PETSC_NULL);
301 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-spring_stiffness",
302 &spring_stiffness, PETSC_NULL);
303 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-alpha_dumping",
304 &alpha_dumping, PETSC_NULL);
305
306 MOFEM_LOG("CONTACT", Sev::inform) << "Young modulus " << young_modulus;
307 MOFEM_LOG("CONTACT", Sev::inform) << "Poisson_ratio " << poisson_ratio;
308 MOFEM_LOG("CONTACT", Sev::inform) << "Density " << rho;
309 MOFEM_LOG("CONTACT", Sev::inform) << "cn " << cn;
310 MOFEM_LOG("CONTACT", Sev::inform)
311 << "spring_stiffness " << spring_stiffness;
312 MOFEM_LOG("CONTACT", Sev::inform) << "alpha_dumping " << alpha_dumping;
313
315 };
316
317 CHKERR get_options();
318
319#ifdef PYTHON_SFD
320 sdfPythonPtr = boost::make_shared<SDFPython>();
321 CHKERR sdfPythonPtr->sdfInit("sdf.py");
322 sdfPythonWeakPtr = sdfPythonPtr;
323#endif
324
326}
327//! [Create common data]
328
329//! [Boundary condition]
332 auto bc_mng = mField.getInterface<BcManager>();
334
335 for (auto f : {"U", "SIGMA"}) {
336 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
337 "REMOVE_X", f, 0, 0);
338 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
339 "REMOVE_Y", f, 1, 1);
340 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
341 "REMOVE_Z", f, 2, 2);
342 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
343 "REMOVE_ALL", f, 0, 3);
344 }
345
346 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_X",
347 "SIGMA", 0, 0, false, true);
348 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Y",
349 "SIGMA", 1, 1, false, true);
350 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Z",
351 "SIGMA", 2, 2, false, true);
352 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
353 "SIGMA", 0, 3, false, true);
354 CHKERR bc_mng->removeBlockDOFsOnEntities(
355 simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
356
357 // Note remove has to be always before push. Then node marking will be
358 // corrupted.
359 CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
360 simple->getProblemName(), "U");
361 boundaryMarker = bc_mng->getMergedBlocksMarker(vector<string>{"FIX_", "ROTATE_"});
362
364}
365//! [Boundary condition]
366
367//! [Push operators to pip]
371 auto *pip_mng = mField.getInterface<PipelineManager>();
372 auto bc_mng = mField.getInterface<BcManager>();
373 auto time_scale = boost::make_shared<TimeScale>();
374
375 auto add_domain_base_ops = [&](auto &pip) {
377 "GEOMETRY"*/);
378 };
379
380 auto common_data_ptr = boost::make_shared<ContactOps::CommonData>();
381 auto henky_common_data_ptr = boost::make_shared<HenckyOps::CommonData>();
382 henky_common_data_ptr->matGradPtr = common_data_ptr->mGradPtr();
383 henky_common_data_ptr->matDPtr = common_data_ptr->mDPtr();
384
385 auto add_domain_ops_lhs = [&](auto &pip) {
386 pip.push_back(new OpSetBc("U", true, boundaryMarker));
387
388 CHKERR addMatBlockOps(mField, pip, "U", "MAT_ELASTIC",
389 common_data_ptr->mDPtr(), Sev::verbose);
390
392 "U", common_data_ptr->mGradPtr()));
393 pip.push_back(
394 new OpCalculateEigenVals<SPACE_DIM>("U", henky_common_data_ptr));
395 pip.push_back(new OpCalculateLogC<SPACE_DIM>("U", henky_common_data_ptr));
396 pip.push_back(
397 new OpCalculateLogC_dC<SPACE_DIM>("U", henky_common_data_ptr));
398 pip.push_back(
399 new OpCalculateHenckyStress<SPACE_DIM>("U", henky_common_data_ptr));
400 pip.push_back(
401 new OpCalculatePiolaStress<SPACE_DIM>("U", henky_common_data_ptr));
402 pip.push_back(new OpHenckyTangent<SPACE_DIM>("U", henky_common_data_ptr));
403 pip.push_back(
404 new OpKPiola("U", "U", henky_common_data_ptr->getMatTangent()));
405
406 if (!is_quasi_static) {
407 auto get_inertia_and_mass_dumping = [this](const double, const double,
408 const double) {
409 auto *pip_mng = mField.getInterface<PipelineManager>();
410 auto &fe_domain_lhs = pip_mng->getDomainLhsFE();
411 return rho * fe_domain_lhs->ts_aa + alpha_dumping * fe_domain_lhs->ts_a;
412 };
413 pip.push_back(new OpMass("U", "U", get_inertia_and_mass_dumping));
414 } else if (alpha_dumping > 0) {
415 auto get_mass_dumping = [this](const double, const double,
416 const double) {
417 auto *pip_mng = mField.getInterface<PipelineManager>();
418 auto &fe_domain_lhs = pip_mng->getDomainLhsFE();
419 return alpha_dumping * fe_domain_lhs->ts_a;
420 };
421 pip.push_back(new OpMass("U", "U", get_mass_dumping));
422 }
423
424 auto unity = []() { return 1; };
425 pip.push_back(new OpMixDivULhs("SIGMA", "U", unity, true));
426 pip.push_back(new OpLambdaGraULhs("SIGMA", "U", unity, true));
427
428 pip.push_back(new OpUnSetBc("U"));
429 };
430
431 auto add_domain_ops_rhs = [&](auto &pip) {
432 pip.push_back(new OpSetBc("U", true, boundaryMarker));
433
435 pip, mField, "U", {time_scale}, Sev::inform);
436
437 CHKERR addMatBlockOps(mField, pip, "U", "MAT_ELASTIC",
438 common_data_ptr->mDPtr(), Sev::inform);
440 "U", common_data_ptr->mGradPtr()));
441
442 pip.push_back(
443 new OpCalculateEigenVals<SPACE_DIM>("U", henky_common_data_ptr));
444 pip.push_back(new OpCalculateLogC<SPACE_DIM>("U", henky_common_data_ptr));
445 pip.push_back(
446 new OpCalculateLogC_dC<SPACE_DIM>("U", henky_common_data_ptr));
447 pip.push_back(
448 new OpCalculateHenckyStress<SPACE_DIM>("U", henky_common_data_ptr));
449 pip.push_back(
450 new OpCalculatePiolaStress<SPACE_DIM>("U", henky_common_data_ptr));
451 pip.push_back(new OpInternalForcePiola(
452 "U", henky_common_data_ptr->getMatFirstPiolaStress()));
453
455 "U", common_data_ptr->contactDispPtr()));
456
458 "SIGMA", common_data_ptr->contactStressPtr()));
460 "SIGMA", common_data_ptr->contactStressDivergencePtr()));
461
462 pip.push_back(new OpMixDivURhs("SIGMA", common_data_ptr->contactDispPtr(),
463 [](double, double, double) { return 1; }));
464 pip.push_back(
465 new OpMixLambdaGradURhs("SIGMA", common_data_ptr->mGradPtr()));
466
467 pip.push_back(new OpMixUTimesDivLambdaRhs(
468 "U", common_data_ptr->contactStressDivergencePtr()));
469 pip.push_back(
470 new OpMixUTimesLambdaRhs("U", common_data_ptr->contactStressPtr()));
471
472 // only in case of dynamics
473 if (!is_quasi_static) {
474 auto mat_acceleration = boost::make_shared<MatrixDouble>();
476 "U", mat_acceleration));
477 pip.push_back(new OpInertiaForce(
478 "U", mat_acceleration, [](double, double, double) { return rho; }));
479 }
480 if (alpha_dumping > 0) {
481 auto mat_velocity = boost::make_shared<MatrixDouble>();
482 pip.push_back(
483 new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));
484 pip.push_back(
485 new OpInertiaForce("U", mat_velocity, [](double, double, double) {
486 return alpha_dumping;
487 }));
488 }
489 pip.push_back(new OpUnSetBc("U"));
490 };
491
492 auto add_boundary_base_ops = [&](auto &pip) {
494 "GEOMETRY"*/);
496 "U", common_data_ptr->contactDispPtr()));
498 "SIGMA", common_data_ptr->contactTractionPtr()));
499 };
500
501 auto add_boundary_ops_lhs = [&](auto &pip) {
505 mField, pip, simple->getProblemName(), "U");
506 pip.push_back(new OpSetBc("U", true, boundaryMarker));
508 pip, mField, "U", Sev::inform);
509 pip.push_back(new OpConstrainBoundaryLhs_dU("SIGMA", "U", common_data_ptr));
510 pip.push_back(new OpConstrainBoundaryLhs_dTraction("SIGMA", "SIGMA",
511 common_data_ptr));
512
513 if (spring_stiffness > 0)
514 pip.push_back(new OpSpringLhs(
515 "U", "U",
516
517 [this](double, double, double) { return spring_stiffness; }
518
519 ));
520
521 pip.push_back(new OpUnSetBc("U"));
523 };
524
525 auto add_boundary_ops_rhs = [&](auto &pip) {
527
528 auto u_mat_ptr = boost::make_shared<MatrixDouble>();
529 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_mat_ptr));
532 mField, pip, simple->getProblemName(), "U", u_mat_ptr,
533 {boost::make_shared<TimeScale>()}); // note displacements have no
534 // scaling
535
536 pip.push_back(new OpSetBc("U", true, boundaryMarker));
538 pip, mField, "U", {time_scale}, Sev::inform);
539 pip.push_back(new OpConstrainBoundaryRhs("SIGMA", common_data_ptr));
540 if (spring_stiffness > 0)
541 pip.push_back(new OpSpringRhs(
542 "U", common_data_ptr->contactDispPtr(),
543 [this](double, double, double) { return spring_stiffness; }));
544 pip.push_back(new OpUnSetBc("U"));
546 };
547
548 add_domain_base_ops(pip_mng->getOpDomainLhsPipeline());
549 add_domain_base_ops(pip_mng->getOpDomainRhsPipeline());
550 add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
551 add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
552
553 add_boundary_base_ops(pip_mng->getOpBoundaryLhsPipeline());
554 add_boundary_base_ops(pip_mng->getOpBoundaryRhsPipeline());
555 CHKERR add_boundary_ops_lhs(pip_mng->getOpBoundaryLhsPipeline());
556 CHKERR add_boundary_ops_rhs(pip_mng->getOpBoundaryRhsPipeline());
557
558 auto integration_rule_vol = [](int, int, int approx_order) {
559 return 2 * approx_order + geom_order - 1;
560 };
561 CHKERR pip_mng->setDomainRhsIntegrationRule(integration_rule_vol);
562 CHKERR pip_mng->setDomainLhsIntegrationRule(integration_rule_vol);
563 auto integration_rule_boundary = [](int, int, int approx_order) {
564 return 2 * approx_order + geom_order - 1;
565 };
566 CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_boundary);
567 CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_boundary);
568
570}
571//! [Push operators to pip]
572
573//! [Solve]
576
579 ISManager *is_manager = mField.getInterface<ISManager>();
580
581 auto set_section_monitor = [&](auto solver) {
583 SNES snes;
584 CHKERR TSGetSNES(solver, &snes);
585 PetscViewerAndFormat *vf;
586 CHKERR PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,
587 PETSC_VIEWER_DEFAULT, &vf);
588 CHKERR SNESMonitorSet(
589 snes,
590 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal, void *))SNESMonitorFields,
591 vf, (MoFEMErrorCode(*)(void **))PetscViewerAndFormatDestroy);
593 };
594
595 auto scatter_create = [&](auto D, auto coeff) {
597 CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
598 ROW, "U", coeff, coeff, is);
599 int loc_size;
600 CHKERR ISGetLocalSize(is, &loc_size);
601 Vec v;
602 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
603 VecScatter scatter;
604 CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
605 return std::make_tuple(SmartPetscObj<Vec>(v),
607 };
608
609 auto set_time_monitor = [&](auto dm, auto solver) {
611 boost::shared_ptr<Monitor> monitor_ptr(
613 boost::shared_ptr<ForcesAndSourcesCore> null;
614 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
615 monitor_ptr, null, null);
617 };
618
619 auto set_fieldsplit_preconditioner = [&](auto solver) {
621
622 SNES snes;
623 CHKERR TSGetSNES(solver, &snes);
624 KSP ksp;
625 CHKERR SNESGetKSP(snes, &ksp);
626 PC pc;
627 CHKERR KSPGetPC(ksp, &pc);
628 PetscBool is_pcfs = PETSC_FALSE;
629 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
630
631 // Setup fieldsplit (block) solver - optional: yes/no
632 if (is_pcfs == PETSC_TRUE) {
633 auto bc_mng = mField.getInterface<BcManager>();
634 auto name_prb = simple->getProblemName();
635 auto is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_X", "U", 0, 0);
636 is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_Y", "U", 1, 1, is_all_bc);
637 is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_Z", "U", 2, 2, is_all_bc);
638 is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_ALL", "U", 0, 2, is_all_bc);
639 int is_all_bc_size;
640 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
641 MOFEM_LOG("CONTACT", Sev::inform)
642 << "Field split block size " << is_all_bc_size;
643 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
644 is_all_bc); // boundary block
645 }
646
648 };
649
650 auto dm = simple->getDM();
651 auto D = createDMVector(dm);
652
654
655 uXScatter = scatter_create(D, 0);
656 uYScatter = scatter_create(D, 1);
657 if (SPACE_DIM == 3)
658 uZScatter = scatter_create(D, 2);
659
660 if (is_quasi_static) {
661 auto solver = pip_mng->createTSIM();
662 auto D = createDMVector(dm);
663 CHKERR set_section_monitor(solver);
664 CHKERR set_time_monitor(dm, solver);
665 CHKERR TSSetSolution(solver, D);
666 CHKERR TSSetFromOptions(solver);
667 CHKERR set_fieldsplit_preconditioner(solver);
668 CHKERR TSSetUp(solver);
669 CHKERR TSSolve(solver, NULL);
670 } else {
671 auto solver = pip_mng->createTSIM2();
672 auto dm = simple->getDM();
673 auto D = createDMVector(dm);
674 auto DD = vectorDuplicate(D);
675 CHKERR set_section_monitor(solver);
676 CHKERR set_time_monitor(dm, solver);
677 CHKERR TS2SetSolution(solver, D, DD);
678 CHKERR TSSetFromOptions(solver);
679 CHKERR set_fieldsplit_preconditioner(solver);
680 CHKERR TSSetUp(solver);
681 CHKERR TSSolve(solver, NULL);
682 }
683
685
686 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
687 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
688 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
689
691}
692//! [Solve]
693
694//! [Postprocess results]
696//! [Postprocess results]
697
698//! [Check]
700//! [Check]
701
702static char help[] = "...\n\n";
703
704int main(int argc, char *argv[]) {
705
706#ifdef PYTHON_SFD
707 Py_Initialize();
708#endif
709
710 // Initialisation of MoFEM/PETSc and MOAB data structures
711 const char param_file[] = "param_file.petsc";
713
714 // Add logging channel for CONTACT
715 auto core_log = logging::core::get();
716 core_log->add_sink(
718 LogManager::setLog("CONTACT");
719 MOFEM_LOG_TAG("CONTACT", "indent");
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 interface
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 //! [CONTACT]
745 Contact ex(m_field);
746 CHKERR ex.runProblem();
747 //! [CONTACT]
748 }
750
752
753#ifdef PYTHON_SFD
754 if (Py_FinalizeEx() < 0) {
755 exit(120);
756 }
757#endif
758}
759
761 MoFEM::Interface &m_field,
762 boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator> &pipeline,
763 std::string field_name, std::string block_name,
764 boost::shared_ptr<MatrixDouble> mat_D_Ptr, Sev sev) {
766
767 struct OpMatBlocks : public DomainEleOp {
768 OpMatBlocks(std::string field_name, boost::shared_ptr<MatrixDouble> m,
769 double bulk_modulus_K, double shear_modulus_G,
770 MoFEM::Interface &m_field, Sev sev,
771 std::vector<const CubitMeshSets *> meshset_vec_ptr)
773 bulkModulusKDefault(bulk_modulus_K),
774 shearModulusGDefault(shear_modulus_G) {
775 std::fill(&(doEntities[MBEDGE]), &(doEntities[MBMAXTYPE]), false);
776 CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
777 "Can not get data from block");
778 }
779
780 MoFEMErrorCode doWork(int side, EntityType type,
783
784 for (auto &b : blockData) {
785
786 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
787 CHKERR getMatDPtr(matDPtr, b.bulkModulusK, b.shearModulusG);
789 }
790 }
791
792 CHKERR getMatDPtr(matDPtr, bulkModulusKDefault, shearModulusGDefault);
794 }
795
796 private:
797 boost::shared_ptr<MatrixDouble> matDPtr;
798
799 struct BlockData {
800 double bulkModulusK;
801 double shearModulusG;
802 Range blockEnts;
803 };
804
805 double bulkModulusKDefault;
806 double shearModulusGDefault;
807 std::vector<BlockData> blockData;
808
810 extractBlockData(MoFEM::Interface &m_field,
811 std::vector<const CubitMeshSets *> meshset_vec_ptr,
812 Sev sev) {
814
815 for (auto m : meshset_vec_ptr) {
816 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
817 std::vector<double> block_data;
818 CHKERR m->getAttributes(block_data);
819 if (block_data.size() != 2) {
820 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
821 "Expected that block has two attribute");
822 }
823 auto get_block_ents = [&]() {
824 Range ents;
825 CHKERR
826 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
827 return ents;
828 };
829
830 double young_modulus = block_data[0];
831 double poisson_ratio = block_data[1];
832 double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
833 double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
834
835 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
836 << "E = " << young_modulus << " nu = " << poisson_ratio;
837
838 blockData.push_back(
839 {bulk_modulus_K, shear_modulus_G, get_block_ents()});
840 }
841 MOFEM_LOG_CHANNEL("WORLD");
843 }
844
845 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
846 double bulk_modulus_K, double shear_modulus_G) {
848 //! [Calculate elasticity tensor]
849 auto set_material_stiffness = [&]() {
855 double A = (SPACE_DIM == 2)
856 ? 2 * shear_modulus_G /
857 (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
858 : 1;
859 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
860 t_D(i, j, k, l) =
861 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
862 A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
863 t_kd(k, l);
864 };
865 //! [Calculate elasticity tensor]
866 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
867 mat_D_ptr->resize(size_symm * size_symm, 1);
868 set_material_stiffness();
870 }
871 };
872
873 double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
874 double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
875 pipeline.push_back(new OpMatBlocks(
876 field_name, mat_D_Ptr, bulk_modulus_K, shear_modulus_G, m_field, sev,
877
878 // Get blockset using regular expression
879 m_field.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
880
881 (boost::format("%s(.*)") % block_name).str()
882
883 ))
884
885 ));
886
888}
std::string param_file
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:362
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.
double young_modulus
Definition: contact.cpp:121
double alpha_dumping
Definition: contact.cpp:127
static char help[]
[Check]
Definition: contact.cpp:702
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< GAUSS >::OpMixDivTimesU< 3, SPACE_DIM, SPACE_DIM > OpMixDivURhs
Definition: contact.cpp:62
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< I >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpSpringRhs
Definition: contact.cpp:72
double spring_stiffness
Definition: contact.cpp:125
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< I >::OpMixDivTimesVec< SPACE_DIM > OpMixDivULhs
[Operators used for contact]
Definition: contact.cpp:58
constexpr int geom_order
Definition: contact.cpp:118
double rho
Definition: contact.cpp:123
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< I >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: contact.cpp:93
#define EXECUTABLE_DIMENSION
Definition: contact.cpp:11
constexpr int SPACE_DIM
Definition: contact.cpp:42
double poisson_ratio
Definition: contact.cpp:122
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< I >::OpMixTensorTimesGradU< SPACE_DIM > OpMixLambdaGradURhs
Definition: contact.cpp:64
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Operators used for contact]
Definition: contact.cpp:77
constexpr IntegrationType I
Definition: contact.cpp:26
constexpr AssemblyType A
Definition: contact.cpp:25
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< I >::OpMixTensorTimesGrad< SPACE_DIM > OpLambdaGraULhs
Definition: contact.cpp:60
int order
Definition: contact.cpp:117
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< I >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpMixUTimesLambdaRhs
Definition: contact.cpp:68
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Essential boundary conditions]
Definition: contact.cpp:91
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< I >::OpMixVecTimesDivLambda< SPACE_DIM > OpMixUTimesDivLambdaRhs
Definition: contact.cpp:66
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< I >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
Definition: contact.cpp:79
double cn
Definition: contact.cpp:124
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< I >::OpMass< 1, SPACE_DIM > OpSpringLhs
Definition: contact.cpp:70
constexpr bool is_quasi_static
Definition: contact.cpp:115
constexpr FieldSpace CONTACT_SPACE
Definition: contact.cpp:54
@ 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
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
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:79
constexpr double bulk_modulus_K
Definition: elastic.cpp:78
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'm', SPACE_DIM > m
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:509
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
IntegrationType
Form integrator integration types.
AssemblyType
[Storage and set boundary conditions]
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
FTensor::Index< 'i', SPACE_DIM > i
double D
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 >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< 3 > OpLambdaGraULhs
MoFEMErrorCode addMatBlockOps(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string field_name, std::string block_name, boost::shared_ptr< MatrixDouble > mat_D_Ptr, Sev sev)
Definition: contact.cpp:760
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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:1042
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
DomainEle::UserDataOperator DomainEleOp
constexpr AssemblyType A
double young_modulus
Definition: plastic.cpp:99
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: plastic.cpp:68
constexpr auto size_symm
Definition: plastic.cpp:33
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: plastic.cpp:66
constexpr auto field_name
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition: radiation.cpp:29
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: contact.cpp:158
MoFEMErrorCode tsSolve()
[Push operators to pip]
Definition: contact.cpp:574
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: contact.cpp:156
MoFEMErrorCode runProblem()
[Run problem]
Definition: contact.cpp:166
MoFEMErrorCode setupProblem()
[Run problem]
Definition: contact.cpp:180
MoFEMErrorCode OPs()
[Boundary condition]
Definition: contact.cpp:368
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: contact.cpp:290
MoFEMErrorCode postProcess()
[Solve]
Definition: contact.cpp:695
MoFEMErrorCode checkResults()
[Postprocess results]
Definition: contact.cpp:699
MoFEM::Interface & mField
Definition: contact.cpp:145
MoFEMErrorCode bC()
[Create common data]
Definition: contact.cpp:330
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: contact.cpp:155
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: contact.cpp:157
Contact(MoFEM::Interface &m_field)
Definition: contact.cpp:140
static SmartPetscObj< Vec > totalTraction
Definition: ContactOps.hpp:21
static auto createTotalTraction(MoFEM::Interface &m_field)
Definition: ContactOps.hpp:23
Add operators pushing bases from local to physical configuration.
Simple interface for fast problem set-up.
Definition: BcManager.hpp:24
Managing BitRefLevels.
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
virtual MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
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)
Essential boundary conditions.
Definition: Essential.hpp:101
@ OPROW
operator doWork function is executed on FE rows
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:298
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
Interface for managing meshsets containing materials and boundary conditions.
Assembly methods.
Definition: Natural.hpp:65
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.
Enforce essential constrains on lhs.
Definition: Essential.hpp:71
Enforce essential constrains on rhs.
Definition: Essential.hpp:55
Set indices on entities on finite element.
PipelineManager interface.
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.