v0.16.0
Loading...
Searching...
No Matches
Classes | Functions | Variables
nonlinear_dynamics.cpp File Reference
#include <BasicFiniteElements.hpp>
#include <ElasticMaterials.hpp>
#include <SurfacePressureComplexForLazy.hpp>
#include <TimeForceScale.hpp>

Go to the source code of this file.

Classes

struct  MonitorPostProc
 
struct  MonitorRestart
 

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "...\n\n"
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)

Definition at line 244 of file nonlinear_dynamics.cpp.

244 {
245
246 const string default_options = "-ksp_type fgmres \n"
247 "-pc_type lu \n"
248 "-pc_factor_mat_solver_type mumps \n"
249 "-mat_mumps_icntl_20 0 \n"
250 "-ksp_atol 1e-10 \n"
251 "-ksp_rtol 1e-10 \n"
252 "-snes_monitor \n"
253 "-snes_type newtonls \n"
254 "-snes_linesearch_type basic \n"
255 "-snes_max_it 100 \n"
256 "-snes_atol 1e-7 \n"
257 "-snes_rtol 1e-7 \n"
258 "-ts_monitor \n"
259 "-ts_type alpha \n";
260
261 string param_file = "param_file.petsc";
262 if (!static_cast<bool>(ifstream(param_file))) {
263 std::ofstream file(param_file.c_str(), std::ios::ate);
264 if (file.is_open()) {
265 file << default_options;
266 file.close();
267 }
268 }
269
270 MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
271
272 // Add logging channel for example
273 auto core_log = logging::core::get();
274 core_log->add_sink(
276 LogManager::setLog("DYNAMIC");
277 MOFEM_LOG_TAG("DYNAMIC", "dynamic");
278
279 try {
280
281 moab::Core mb_instance;
282 moab::Interface &moab = mb_instance;
283
284 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
285 auto moab_comm_wrap =
286 boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
287 if (pcomm == NULL)
288 pcomm = new ParallelComm(&moab, moab_comm_wrap->get_comm());
289
291 char mesh_file_name[255];
292 PetscBool is_partitioned = PETSC_FALSE;
293 PetscBool linear = PETSC_TRUE;
294 PetscInt disp_order = 1;
295 PetscInt vel_order = 1;
296 PetscBool is_solve_at_time_zero = PETSC_FALSE;
297
298 auto read_command_line_parameters = [&]() {
300 PetscBool flg = PETSC_TRUE;
301 CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-my_file",
302 mesh_file_name, 255, &flg);
303 if (flg != PETSC_TRUE)
304 SETERRQ(PETSC_COMM_SELF, 1, "Error -my_file (mesh file needed)");
305
306 // use this if your mesh is partitioned and you run code on parts,
307 // you can solve very big problems
308 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, PETSC_NULLPTR, "-my_is_partitioned",
309 &is_partitioned, &flg);
310
311 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, PETSC_NULLPTR, "-is_linear", &linear,
312 PETSC_NULLPTR);
313
314 enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, LASBASETOP };
315 const char *list_bases[] = {"legendre", "lobatto", "bernstein_bezier"};
316 PetscInt choice_base_value = BERNSTEIN_BEZIER;
317 CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-base", list_bases,
318 LASBASETOP, &choice_base_value, PETSC_NULLPTR);
319 if (choice_base_value == LEGENDRE)
321 else if (choice_base_value == LOBATTO)
323 else if (choice_base_value == BERNSTEIN_BEZIER)
325
326 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, PETSC_NULLPTR, "-my_disp_order",
327 &disp_order, &flg);
328 if (flg != PETSC_TRUE)
329 disp_order = 1;
330
331 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, PETSC_NULLPTR, "-my_vel_order",
332 &vel_order, &flg);
333 if (flg != PETSC_TRUE)
334 vel_order = disp_order;
335
336 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, PETSC_NULLPTR,
337 "-my_solve_at_time_zero",
338 &is_solve_at_time_zero, &flg);
339
341 };
342
343 auto read_mesh = [&]() {
345 if (is_partitioned == PETSC_TRUE) {
346 // Read mesh to MOAB
347 const char *option;
348 option = "PARALLEL=BCAST_DELETE;"
349 "PARALLEL_RESOLVE_SHARED_ENTS;"
350 "PARTITION=PARALLEL_PARTITION;";
351 CHKERR moab.load_file(mesh_file_name, 0, option);
352 } else {
353 const char *option;
354 option = "";
355 CHKERR moab.load_file(mesh_file_name, 0, option);
356 }
358 };
359
360 CHKERR read_command_line_parameters();
362
363 MoFEM::Core core(moab);
364 MoFEM::Interface &m_field = core;
365
366 // ref meshset ref level 0
367 BitRefLevel bit_level0;
368 bit_level0.set(0);
369 EntityHandle meshset_level0;
370 CHKERR moab.create_meshset(MESHSET_SET, meshset_level0);
371 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
372 0, 3, bit_level0);
373 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByRefLevel(
374 bit_level0, BitRefLevel().set(), meshset_level0);
375
376 // Fields
377 CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
378 3, MB_TAG_SPARSE, MF_ZERO);
379 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
380 CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 2);
381 CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 2);
382 CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
383 CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
384
385 bool check_if_spatial_field_exist = m_field.check_field("DISPLACEMENT");
386 CHKERR m_field.add_field("DISPLACEMENT", H1, base, 3, MB_TAG_SPARSE,
387 MF_ZERO);
388 // add entities (by tets) to the field
389 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "DISPLACEMENT");
390
391 // set app. order
392 CHKERR m_field.set_field_order(0, MBTET, "DISPLACEMENT", disp_order);
393 CHKERR m_field.set_field_order(0, MBTRI, "DISPLACEMENT", disp_order);
394 CHKERR m_field.set_field_order(0, MBEDGE, "DISPLACEMENT", disp_order);
396 CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", disp_order);
397 else
398 CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
399
400 // Add nodal force element
401 CHKERR MetaNeumannForces::addNeumannBCElements(m_field, "DISPLACEMENT");
402 CHKERR MetaEdgeForces::addElement(m_field, "DISPLACEMENT");
403 CHKERR MetaNodalForces::addElement(m_field, "DISPLACEMENT");
404 // Add fluid pressure finite elements
405 FluidPressure fluid_pressure_fe(m_field);
406 fluid_pressure_fe.addNeumannFluidPressureBCElements("DISPLACEMENT");
408 fluid_pressure_fe.getLoopFe().getOpPtrVector(), {},
409 "MESH_NODE_POSITIONS");
410 fluid_pressure_fe.setNeumannFluidPressureFiniteElementOperators(
411 "DISPLACEMENT", PETSC_NULLPTR, false, true);
412
413 // Velocity
414 CHKERR m_field.add_field("VELOCITY", H1, base, 3, MB_TAG_SPARSE, MF_ZERO);
415 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "VELOCITY");
416
417 CHKERR m_field.set_field_order(0, MBTET, "VELOCITY", vel_order);
418 CHKERR m_field.set_field_order(0, MBTRI, "VELOCITY", vel_order);
419 CHKERR m_field.set_field_order(0, MBEDGE, "VELOCITY", vel_order);
421 CHKERR m_field.set_field_order(0, MBVERTEX, "VELOCITY", vel_order);
422 else
423 CHKERR m_field.set_field_order(0, MBVERTEX, "VELOCITY", 1);
424
425 CHKERR m_field.add_field("DOT_DISPLACEMENT", H1, base, 3, MB_TAG_SPARSE,
426 MF_ZERO);
427 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "DOT_DISPLACEMENT");
428 CHKERR m_field.set_field_order(0, MBTET, "DOT_DISPLACEMENT", disp_order);
429 CHKERR m_field.set_field_order(0, MBTRI, "DOT_DISPLACEMENT", disp_order);
430 CHKERR m_field.set_field_order(0, MBEDGE, "DOT_DISPLACEMENT", disp_order);
432 CHKERR m_field.set_field_order(0, MBVERTEX, "DOT_DISPLACEMENT",
433 disp_order);
434 else
435 CHKERR m_field.set_field_order(0, MBVERTEX, "DOT_DISPLACEMENT", 1);
436
437 CHKERR m_field.add_field("DOT_VELOCITY", H1, base, 3, MB_TAG_SPARSE,
438 MF_ZERO);
439 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "DOT_VELOCITY");
440 CHKERR m_field.set_field_order(0, MBTET, "DOT_VELOCITY", vel_order);
441 CHKERR m_field.set_field_order(0, MBTRI, "DOT_VELOCITY", vel_order);
442 CHKERR m_field.set_field_order(0, MBEDGE, "DOT_VELOCITY", vel_order);
444 CHKERR m_field.set_field_order(0, MBVERTEX, "DOT_VELOCITY", disp_order);
445 else
446 CHKERR m_field.set_field_order(0, MBVERTEX, "DOT_VELOCITY", 1);
447
448 // Set material model and mass element
449 NonlinearElasticElement elastic(m_field, 2);
450 ElasticMaterials elastic_materials(m_field);
451 CHKERR elastic_materials.setBlocks(elastic.setOfBlocks);
452 // NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI<adouble>
453 // st_venant_kirchhoff_material_adouble;
454 // NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI<double>
455 // st_venant_kirchhoff_material_double; CHKERR
456 // elastic.setBlocks(&st_venant_kirchhoff_material_double,&st_venant_kirchhoff_material_adouble);
457 CHKERR elastic.addElement("ELASTIC", "DISPLACEMENT");
458 CHKERR AddHOOps<3, 3, 3>::add(elastic.getLoopFeRhs().getOpPtrVector(), {H1},
459 "MESH_NODE_POSITIONS");
460 CHKERR AddHOOps<3, 3, 3>::add(elastic.getLoopFeLhs().getOpPtrVector(), {H1},
461 "MESH_NODE_POSITIONS");
463 elastic.getLoopFeEnergy().getOpPtrVector(), {H1},
464 "MESH_NODE_POSITIONS");
465 CHKERR elastic.setOperators("DISPLACEMENT", "MESH_NODE_POSITIONS", false,
466 true);
467
468 // set mass element
469 ConvectiveMassElement inertia(m_field, 1);
470 // CHKERR inertia.setBlocks();
471 CHKERR elastic_materials.setBlocks(inertia.setOfBlocks);
472 CHKERR inertia.addConvectiveMassElement("MASS_ELEMENT", "VELOCITY",
473 "DISPLACEMENT");
474 CHKERR inertia.addVelocityElement("VELOCITY_ELEMENT", "VELOCITY",
475 "DISPLACEMENT");
476
477 // Add possibility to load accelerogram
478 {
479 string name = "-my_accelerogram";
480 char time_file_name[255];
481 PetscBool flg;
482 CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, name.c_str(),
483 time_file_name, 255, &flg);
484 if (flg == PETSC_TRUE) {
485 inertia.methodsOp.push_back(new TimeAccelerogram(name));
486 }
487 }
488
489 // damper element
490 KelvinVoigtDamper damper(m_field);
491 CHKERR elastic_materials.setBlocks(damper.blockMaterialDataMap);
492 {
493 KelvinVoigtDamper::CommonData &common_data = damper.commonData;
494 common_data.spatialPositionName = "DISPLACEMENT";
495 common_data.spatialPositionNameDot = "DOT_DISPLACEMENT";
496 CHKERR m_field.add_finite_element("DAMPER", MF_ZERO);
498 "DISPLACEMENT");
500 "DISPLACEMENT");
502 "DISPLACEMENT");
503
504 if (m_field.check_field("MESH_NODE_POSITIONS")) {
506 "DAMPER", "MESH_NODE_POSITIONS");
507 }
508 std::map<int, KelvinVoigtDamper::BlockMaterialData>::iterator bit =
509 damper.blockMaterialDataMap.begin();
510 for (; bit != damper.blockMaterialDataMap.end(); bit++) {
511 bit->second.lInear = linear;
512 int id = bit->first;
513 KelvinVoigtDamper::BlockMaterialData &material_data = bit->second;
514 damper.constitutiveEquationMap.insert(
516 material_data));
517 CHKERR m_field.add_ents_to_finite_element_by_type(bit->second.tEts,
518 MBTET, "DAMPER");
519 }
520 CHKERR damper.setOperators(3);
521 }
522
523 MonitorPostProc post_proc(m_field, elastic.setOfBlocks,
524 elastic.getLoopFeEnergy(),
525 inertia.getLoopFeEnergy());
526
527 // elastic and mass element calculated in Kuu shell matrix problem. To
528 // calculate Mass element, velocity field is needed.
529 CHKERR m_field.modify_finite_element_add_field_data("ELASTIC", "VELOCITY");
531 "DOT_DISPLACEMENT");
533 "DOT_VELOCITY");
534
535 // build field
536 CHKERR m_field.build_fields();
537 // CHKERR m_field.list_dofs_by_field_name("DISPLACEMENT");
538
539 // 10 node tets
540 if (!check_if_spatial_field_exist) {
541 Projection10NodeCoordsOnField ent_method_material(m_field,
542 "MESH_NODE_POSITIONS");
543 CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
544 }
545
546 // build finite elements
548 // build adjacencies
549 CHKERR m_field.build_adjacencies(bit_level0);
550
551 // define problems
552 {
553 CHKERR m_field.add_problem("Kuu", MF_ZERO);
554 CHKERR m_field.modify_problem_add_finite_element("Kuu", "ELASTIC");
555 CHKERR m_field.modify_problem_add_finite_element("Kuu", "PRESSURE_FE");
556 CHKERR m_field.modify_problem_add_finite_element("Kuu", "FORCE_FE");
558 "FLUID_PRESSURE_FE");
559 CHKERR m_field.modify_problem_ref_level_add_bit("Kuu", bit_level0);
560
561 ProblemsManager *prb_mng_ptr;
562 CHKERR m_field.getInterface(prb_mng_ptr);
563 if (is_partitioned) {
564 CHKERR prb_mng_ptr->buildProblemOnDistributedMesh("Kuu", true);
565 CHKERR prb_mng_ptr->partitionFiniteElements("Kuu", true, 0,
566 pcomm->size());
567 } else {
568 CHKERR prb_mng_ptr->buildProblem("Kuu", true);
569 CHKERR prb_mng_ptr->partitionProblem("Kuu");
570 CHKERR prb_mng_ptr->partitionFiniteElements("Kuu");
571 }
572 CHKERR prb_mng_ptr->partitionGhostDofs("Kuu");
573 }
574
575 CHKERR m_field.add_problem("DYNAMICS", MF_ZERO);
576 // set finite elements for problems
577 CHKERR m_field.modify_problem_add_finite_element("DYNAMICS", "ELASTIC");
578 CHKERR m_field.modify_problem_add_finite_element("DYNAMICS", "DAMPER");
579 CHKERR m_field.modify_problem_add_finite_element("DYNAMICS", "PRESSURE_FE");
580 CHKERR m_field.modify_problem_add_finite_element("DYNAMICS", "FORCE_FE");
581 CHKERR m_field.modify_problem_add_finite_element("DYNAMICS",
582 "FLUID_PRESSURE_FE");
583 CHKERR m_field.modify_problem_add_finite_element("DYNAMICS",
584 "MASS_ELEMENT");
585 CHKERR m_field.modify_problem_add_finite_element("DYNAMICS",
586 "VELOCITY_ELEMENT");
587 // set refinement level for problem
588 CHKERR m_field.modify_problem_ref_level_add_bit("DYNAMICS", bit_level0);
589
590 ProblemsManager *prb_mng_ptr;
591 CHKERR m_field.getInterface(prb_mng_ptr);
592 if (is_partitioned) {
593 CHKERR prb_mng_ptr->buildProblemOnDistributedMesh("DYNAMICS", true);
594 CHKERR prb_mng_ptr->partitionFiniteElements("DYNAMICS", true, 0,
595 pcomm->size());
596 } else {
597 CHKERR prb_mng_ptr->buildProblem("DYNAMICS", true);
598 CHKERR prb_mng_ptr->partitionProblem("DYNAMICS");
599 CHKERR prb_mng_ptr->partitionFiniteElements("DYNAMICS");
600 }
601 CHKERR prb_mng_ptr->partitionGhostDofs("DYNAMICS");
602
603 Vec F;
604 CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("DYNAMICS", COL,
605 &F);
606 Vec D;
607 CHKERR VecDuplicate(F, &D);
608
609 // create tS
610 TS ts;
611 CHKERR TSCreate(PETSC_COMM_WORLD, &ts);
612 CHKERR TSSetType(ts, TSBEULER);
613
614 // shell matrix
618 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("Kuu",
619 &shellAij_ctx->K);
620 CHKERR MatDuplicate(shellAij_ctx->K, MAT_DO_NOT_COPY_VALUES,
621 &shellAij_ctx->M);
622 CHKERR shellAij_ctx->iNit();
623 CHKERR m_field.getInterface<VecManager>()->vecScatterCreate(
624 D, "DYNAMICS", COL, shellAij_ctx->u, "Kuu", COL,
625 &shellAij_ctx->scatterU);
626 CHKERR m_field.getInterface<VecManager>()->vecScatterCreate(
627 D, "DYNAMICS", "VELOCITY", COL, shellAij_ctx->v, "Kuu", "DISPLACEMENT",
628 COL, &shellAij_ctx->scatterV);
629 Mat shell_Aij;
630 const Problem *problem_ptr;
631 CHKERR m_field.get_problem("DYNAMICS", &problem_ptr);
632 CHKERR MatCreateShell(
633 PETSC_COMM_WORLD, problem_ptr->getNbLocalDofsRow(),
634 problem_ptr->getNbLocalDofsCol(), problem_ptr->getNbDofsRow(),
635 problem_ptr->getNbDofsRow(), (void *)shellAij_ctx, &shell_Aij);
636 CHKERR MatShellSetOperation(shell_Aij, MATOP_MULT,
637 (void (*)(void))ConvectiveMassElement::MultOpA);
638 CHKERR MatShellSetOperation(
639 shell_Aij, MATOP_ZERO_ENTRIES,
641 // blocked problem
642 ConvectiveMassElement::ShellMatrixElement shell_matrix_element(m_field);
643 DirichletDisplacementBc shell_dirichlet_bc(
644 m_field, "DISPLACEMENT", shellAij_ctx->barK, PETSC_NULLPTR, PETSC_NULLPTR);
645 DirichletDisplacementBc my_dirichlet_bc(m_field, "DISPLACEMENT", PETSC_NULLPTR,
646 D, F);
647 shell_matrix_element.problemName = "Kuu";
648 shell_matrix_element.shellMatCtx = shellAij_ctx;
649 shell_matrix_element.DirichletBcPtr = &shell_dirichlet_bc;
650 shell_matrix_element.loopK.push_back(
651 ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
652 "ELASTIC", &elastic.getLoopFeLhs()));
653 // damper
654 shell_matrix_element.loopK.push_back(
655 ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
656 "ELASTIC", &damper.feLhs));
657
658 CHKERR inertia.addHOOps();
659 CHKERR inertia.setShellMatrixMassOperators("VELOCITY", "DISPLACEMENT",
660 "MESH_NODE_POSITIONS", linear);
661 // element name "ELASTIC" is used, therefore M matrix is assembled as K
662 // matrix. This is added to M is shell matrix. M matrix is a derivative of
663 // inertia forces over spatial velocities
664 shell_matrix_element.loopM.push_back(
665 ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
666 "ELASTIC", &inertia.getLoopFeMassLhs()));
667 // this calculate derivatives of inertia forces over spatial positions and
668 // add this to shell K matrix
669 shell_matrix_element.loopAuxM.push_back(
670 ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
671 "ELASTIC", &inertia.getLoopFeMassAuxLhs()));
672
673 // Element to calculate shell matrix residual
674 ConvectiveMassElement::ShellResidualElement shell_matrix_residual(m_field);
675 shell_matrix_residual.shellMatCtx = shellAij_ctx;
676
677 // surface pressure
678 boost::ptr_map<std::string, NeumannForcesSurface> surface_forces;
679 {
680 string fe_name_str = "FORCE_FE";
681 surface_forces.insert(fe_name_str, new NeumannForcesSurface(m_field));
683 surface_forces.at(fe_name_str).getLoopFe().getOpPtrVector(), {},
684 "MESH_NODE_POSITIONS");
686 NODESET | FORCESET, it)) {
687 CHKERR surface_forces.at(fe_name_str)
688 .addForce("DISPLACEMENT", PETSC_NULLPTR, it->getMeshsetId(), true);
689 surface_forces.at(fe_name_str)
690 .methodsOp.push_back(new TimeForceScale());
691 }
692 }
693
694 boost::ptr_map<std::string, NeumannForcesSurface> surface_pressure;
695 {
696 string fe_name_str = "PRESSURE_FE";
697 surface_pressure.insert(fe_name_str, new NeumannForcesSurface(m_field));
699 surface_pressure.at(fe_name_str).getLoopFe().getOpPtrVector(), {},
700 "MESH_NODE_POSITIONS");
702 m_field, SIDESET | PRESSURESET, it)) {
703 CHKERR surface_pressure.at(fe_name_str)
704 .addPressure("DISPLACEMENT", PETSC_NULLPTR, it->getMeshsetId(), true);
705 surface_pressure.at(fe_name_str)
706 .methodsOp.push_back(new TimeForceScale());
707 }
708 }
709
710 // edge forces
711 boost::ptr_map<std::string, EdgeForce> edge_forces;
712 {
713 string fe_name_str = "FORCE_FE";
714 edge_forces.insert(fe_name_str, new EdgeForce(m_field));
716 NODESET | FORCESET, it)) {
717 CHKERR edge_forces.at(fe_name_str)
718 .addForce("DISPLACEMENT", PETSC_NULLPTR, it->getMeshsetId(), true);
719 edge_forces.at(fe_name_str).methodsOp.push_back(new TimeForceScale());
720 }
721 }
722
723 // nodal forces
724 boost::ptr_map<std::string, NodalForce> nodal_forces;
725 {
726 string fe_name_str = "FORCE_FE";
727 nodal_forces.insert(fe_name_str, new NodalForce(m_field));
729 NODESET | FORCESET, it)) {
730 CHKERR nodal_forces.at(fe_name_str)
731 .addForce("DISPLACEMENT", F, it->getMeshsetId(), true);
732 nodal_forces.at(fe_name_str).methodsOp.push_back(new TimeForceScale());
733 }
734 }
735
736 MonitorRestart monitor_restart(m_field, ts);
738 m_field, ts, "VELOCITY", "DISPLACEMENT");
739
740 // TS
741 TsCtx ts_ctx(m_field, "DYNAMICS");
742
743 // right hand side
744 // preprocess
745 ts_ctx.getPreProcessIFunction().push_back(&update_and_control);
746 ts_ctx.getPreProcessIFunction().push_back(&my_dirichlet_bc);
747
748 // fe looops
749 auto &loops_to_do_Rhs = ts_ctx.getLoopsIFunction();
750
751 auto add_static_rhs = [&](auto &loops_to_do_Rhs) {
753 loops_to_do_Rhs.push_back(
754 PairNameFEMethodPtr("ELASTIC", &elastic.getLoopFeRhs()));
755 for (auto fit = surface_forces.begin(); fit != surface_forces.end();
756 fit++) {
757 loops_to_do_Rhs.push_back(
758 PairNameFEMethodPtr(fit->first, &fit->second->getLoopFe()));
759 }
760 for (auto fit = surface_pressure.begin(); fit != surface_pressure.end();
761 fit++) {
762 loops_to_do_Rhs.push_back(
763 PairNameFEMethodPtr(fit->first, &fit->second->getLoopFe()));
764 }
765 for (auto fit = edge_forces.begin(); fit != edge_forces.end(); fit++) {
766 loops_to_do_Rhs.push_back(
767 PairNameFEMethodPtr(fit->first, &fit->second->getLoopFe()));
768 }
769 for (auto fit = nodal_forces.begin(); fit != nodal_forces.end(); fit++) {
770 loops_to_do_Rhs.push_back(
771 PairNameFEMethodPtr(fit->first, &fit->second->getLoopFe()));
772 }
773 loops_to_do_Rhs.push_back(PairNameFEMethodPtr(
774 "FLUID_PRESSURE_FE", &fluid_pressure_fe.getLoopFe()));
776 };
777
778 CHKERR add_static_rhs(loops_to_do_Rhs);
779
780 loops_to_do_Rhs.push_back(PairNameFEMethodPtr("DAMPER", &damper.feRhs));
781 loops_to_do_Rhs.push_back(
782 PairNameFEMethodPtr("MASS_ELEMENT", &inertia.getLoopFeMassRhs()));
783
784 // preporcess
785 // calculate residual for velocities
786 ts_ctx.getPreProcessIFunction().push_back(&shell_matrix_residual);
787 // postprocess
788 ts_ctx.getPostProcessIFunction().push_back(&my_dirichlet_bc);
789
790 // left hand side
791 // preprocess
792 ts_ctx.getPreProcessIJacobian().push_back(&update_and_control);
793 ts_ctx.getPreProcessIJacobian().push_back(&shell_matrix_element);
794 ts_ctx.getPostProcessIJacobian().push_back(&update_and_control);
795 // monitor
796 TsCtx::FEMethodsSequence &loopsMonitor =
798 loopsMonitor.push_back(
799 TsCtx::PairNameFEMethodPtr("MASS_ELEMENT", &post_proc));
800 loopsMonitor.push_back(
801 TsCtx::PairNameFEMethodPtr("MASS_ELEMENT", &monitor_restart));
802
803 CHKERR TSSetIFunction(ts, F, TsSetIFunction, &ts_ctx);
804 CHKERR TSSetIJacobian(ts, shell_Aij, shell_Aij, TsSetIJacobian, &ts_ctx);
805
806 CHKERR TSMonitorSet(ts, TsMonitorSet, &ts_ctx, PETSC_NULLPTR);
807
808 double ftime = 1;
809 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
810 CHKERR TSSetSolution(ts, D);
811 CHKERR TSSetFromOptions(ts);
812 // shell matrix pre-conditioner
813 SNES snes;
814 CHKERR TSGetSNES(ts, &snes);
815 // CHKERR SNESSetFromOptions(snes);
816 KSP ksp;
817 CHKERR SNESGetKSP(snes, &ksp);
818 CHKERR KSPSetFromOptions(ksp);
819 PC pc;
820 CHKERR KSPGetPC(ksp, &pc);
821 CHKERR PCSetType(pc, PCSHELL);
822 ConvectiveMassElement::PCShellCtx pc_shell_ctx(shell_Aij);
823 CHKERR PCShellSetContext(pc, (void *)&pc_shell_ctx);
826 CHKERR PCShellSetDestroy(pc, ConvectiveMassElement::PCShellDestroy);
827
828 CHKERR VecZeroEntries(D);
829 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
830 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
831 CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
832 "DYNAMICS", COL, D, INSERT_VALUES, SCATTER_REVERSE);
833
834 // Solve problem at time Zero
835 if (is_solve_at_time_zero) {
836
837 Mat Aij = shellAij_ctx->K;
838 Vec F;
839 CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("Kuu", COL, &F);
840 Vec D;
841 CHKERR VecDuplicate(F, &D);
842
843 // Set vector for Kuu problem from the mesh data
844 CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
845 "Kuu", COL, D, INSERT_VALUES, SCATTER_FORWARD);
846 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
847 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
848
849 SnesCtx snes_ctx(m_field, "Kuu");
850
851 SNES snes;
852 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
853 CHKERR SNESSetApplicationContext(snes, &snes_ctx);
854 CHKERR SNESSetFunction(snes, F, SnesRhs, &snes_ctx);
855 CHKERR SNESSetJacobian(snes, Aij, Aij, SnesMat, &snes_ctx);
856 CHKERR SNESSetFromOptions(snes);
857
858 DirichletDisplacementBc my_dirichlet_bc(m_field, "DISPLACEMENT",
859 PETSC_NULLPTR, D, F);
860
861 SnesCtx::FEMethodsSequence &loops_to_do_Rhs =
862 snes_ctx.getComputeRhs();
863 snes_ctx.getPreProcComputeRhs().push_back(&my_dirichlet_bc);
864 fluid_pressure_fe.getLoopFe().ts_t = 0;
865 CHKERR add_static_rhs(loops_to_do_Rhs);
866 snes_ctx.getPostProcComputeRhs().push_back(&my_dirichlet_bc);
867
868 SnesCtx::FEMethodsSequence &loops_to_do_Mat =
869 snes_ctx.getSetOperators();
870 snes_ctx.getPreProcSetOperators().push_back(&my_dirichlet_bc);
871 loops_to_do_Mat.push_back(
872 SnesCtx::PairNameFEMethodPtr("ELASTIC", &elastic.getLoopFeLhs()));
873 snes_ctx.getPostProcSetOperators().push_back(&my_dirichlet_bc);
874
875 CHKERR m_field.getInterface<FieldBlas>()->fieldScale(0, "VELOCITY");
876 CHKERR m_field.getInterface<FieldBlas>()->fieldScale(0,
877 "DOT_DISPLACEMENT");
878 CHKERR m_field.getInterface<FieldBlas>()->fieldScale(0, "DOT_VELOCITY");
879
880 CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
881 "Kuu", COL, D, INSERT_VALUES, SCATTER_FORWARD);
882
883 CHKERR SNESSolve(snes, PETSC_NULLPTR, D);
884 int its;
885 CHKERR SNESGetIterationNumber(snes, &its);
886 MOFEM_LOG_C("DYNAMIC", Sev::inform, "number of Newton iterations = %d\n",
887 its);
888
889 // Set data on the mesh
890 CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
891 "Kuu", COL, D, INSERT_VALUES, SCATTER_REVERSE);
892
893 CHKERR VecDestroy(&F);
894 CHKERR VecDestroy(&D);
895 CHKERR SNESDestroy(&snes);
896 }
897
898 if (is_solve_at_time_zero) {
899 CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
900 "DYNAMICS", COL, D, INSERT_VALUES, SCATTER_FORWARD);
901 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
902 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
903 CHKERR TSSetSolution(ts, D);
904 }
905
906#if PETSC_VERSION_GE(3, 7, 0)
907 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
908#endif
909 CHKERR TSSolve(ts, D);
910 CHKERR TSGetTime(ts, &ftime);
911
912 PetscInt steps, snesfails, rejects, nonlinits, linits;
913 CHKERR TSGetTimeStepNumber(ts, &steps);
914 CHKERR TSGetSNESFailures(ts, &snesfails);
915 CHKERR TSGetStepRejections(ts, &rejects);
916 CHKERR TSGetSNESIterations(ts, &nonlinits);
917 CHKERR TSGetKSPIterations(ts, &linits);
918 MOFEM_LOG_C("DYNAMIC", Sev::inform,
919 "steps %d (%d rejected, %D SNES fails), ftime %g, nonlinits "
920 "%d, linits %D\n",
921 steps, rejects, snesfails, ftime, nonlinits, linits);
922 CHKERR TSDestroy(&ts);
923
924 CHKERR VecDestroy(&F);
925 CHKERR VecDestroy(&D);
926 CHKERR MatDestroy(&shellAij_ctx->K);
927 CHKERR MatDestroy(&shellAij_ctx->M);
928 CHKERR VecScatterDestroy(&shellAij_ctx->scatterU);
929 CHKERR VecScatterDestroy(&shellAij_ctx->scatterV);
930 CHKERR MatDestroy(&shell_Aij);
931 delete shellAij_ctx;
932 }
934
936
937 return 0;
938}
#define MOFEM_LOG_C(channel, severity, format,...)
@ COL
#define CATCH_ERRORS
Catch errors.
@ MF_ZERO
FieldApproximationBase
approximation base
Definition definitions.h:58
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
@ AINSWORTH_LOBATTO_BASE
Definition definitions.h:62
@ NOBASE
Definition definitions.h:59
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition definitions.h:64
@ H1
continuous field
Definition definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ PRESSURESET
@ FORCESET
@ NODESET
@ SIDESET
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
@ F
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual bool check_field(const std::string &name) const =0
check if field is in database
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
MoFEMErrorCode buildProblemOnDistributedMesh(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures, assuming that mesh is distributed (collective)
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elements
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
auto bit
set bit
double D
MoFEM::TsCtx * ts_ctx
const FTensor::Tensor2< T, Dim, Dim > Vec
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
Definition TsCtx.cpp:169
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition TsCtx.cpp:263
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
Definition SnesCtx.cpp:491
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
Definition SnesCtx.cpp:227
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode TsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t, Vec F, void *ctx)
Set IFunction for TS solver.
Definition TsCtx.cpp:56
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
read_mesh(file, args)
static char help[]
structure grouping operators and data used for calculation of mass (convective) element \ nonlinear_e...
static MoFEMErrorCode MultOpA(Mat A, Vec x, Vec f)
Mult operator for shell matrix.
static MoFEMErrorCode ZeroEntriesOp(Mat A)
static MoFEMErrorCode PCShellDestroy(PC pc)
static MoFEMErrorCode PCShellApplyOp(PC pc, Vec f, Vec x)
apply pre-conditioner for shell matrix
static MoFEMErrorCode PCShellSetUpOp(PC pc)
Set Dirichlet boundary conditions on displacements.
Force on edges and lines.
Definition EdgeForce.hpp:13
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
Fluid pressure forces.
Common data for nonlinear_elastic_elem model.
Implementation of Kelvin Voigt Damper.
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition EdgeForce.hpp:62
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Add operators pushing bases from local to physical configuration.
Managing BitRefLevels.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
Core (interface) class.
Definition Core.hpp:83
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:68
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:123
Deprecated interface functions.
Basic algebra on fields.
Definition FieldBlas.hpp:21
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Matrix manager is used to build and partition problems.
keeps basic data about problem
DofIdx getNbLocalDofsRow() const
DofIdx getNbDofsRow() const
DofIdx getNbLocalDofsCol() const
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
Interface for nonlinear (SNES) solver.
Definition SnesCtx.hpp:15
MoFEM::FEMethodsSequence FEMethodsSequence
Definition SnesCtx.hpp:18
Interface for Time Stepping (TS) solver.
Definition TsCtx.hpp:17
BasicMethodsSequence & getPostProcessIJacobian()
Get the postProcess to do IJacobian object.
Definition TsCtx.hpp:132
MoFEM::FEMethodsSequence FEMethodsSequence
Definition TsCtx.hpp:26
FEMethodsSequence & getLoopsMonitor()
Get the loops to do Monitor object.
Definition TsCtx.hpp:102
BasicMethodsSequence & getPostProcessIFunction()
Get the postProcess to do IFunction object.
Definition TsCtx.hpp:116
FEMethodsSequence & getLoopsIFunction()
Get the loops to do IFunction object.
Definition TsCtx.hpp:63
BasicMethodsSequence & getPreProcessIFunction()
Get the preProcess to do IFunction object.
Definition TsCtx.hpp:109
BasicMethodsSequence & getPreProcessIJacobian()
Get the preProcess to do IJacobian object.
Definition TsCtx.hpp:125
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Finite element and operators to apply force/pressures applied to surfaces.
Force applied to nodes.
structure grouping operators and data used for calculation of nonlinear elastic element
Force scale operator for reading two columns.

Variable Documentation

◆ help

char help[] = "...\n\n"
static

Definition at line 28 of file nonlinear_dynamics.cpp.