v0.14.0
Loading...
Searching...
No Matches
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 226 of file nonlinear_dynamics.cpp.

226 {
227
228 const string default_options = "-ksp_type fgmres \n"
229 "-pc_type lu \n"
230 "-pc_factor_mat_solver_type mumps \n"
231 "-mat_mumps_icntl_20 0 \n"
232 "-ksp_atol 1e-10 \n"
233 "-ksp_rtol 1e-10 \n"
234 "-snes_monitor \n"
235 "-snes_type newtonls \n"
236 "-snes_linesearch_type basic \n"
237 "-snes_max_it 100 \n"
238 "-snes_atol 1e-7 \n"
239 "-snes_rtol 1e-7 \n"
240 "-ts_monitor \n"
241 "-ts_type alpha \n";
242
243 string param_file = "param_file.petsc";
244 if (!static_cast<bool>(ifstream(param_file))) {
245 std::ofstream file(param_file.c_str(), std::ios::ate);
246 if (file.is_open()) {
247 file << default_options;
248 file.close();
249 }
250 }
251
252 MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
253
254 // Add logging channel for example
255 auto core_log = logging::core::get();
256 core_log->add_sink(
258 LogManager::setLog("DYNAMIC");
259 MOFEM_LOG_TAG("DYNAMIC", "dynamic");
260
261 try {
262
263 moab::Core mb_instance;
264 moab::Interface &moab = mb_instance;
265
266 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
267 auto moab_comm_wrap =
268 boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
269 if (pcomm == NULL)
270 pcomm = new ParallelComm(&moab, moab_comm_wrap->get_comm());
271
273 char mesh_file_name[255];
274 PetscBool is_partitioned = PETSC_FALSE;
275 PetscBool linear = PETSC_TRUE;
276 PetscInt disp_order = 1;
277 PetscInt vel_order = 1;
278 PetscBool is_solve_at_time_zero = PETSC_FALSE;
279
280 auto read_command_line_parameters = [&]() {
282 PetscBool flg = PETSC_TRUE;
283 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
284 mesh_file_name, 255, &flg);
285 if (flg != PETSC_TRUE)
286 SETERRQ(PETSC_COMM_SELF, 1, "Error -my_file (mesh file needed)");
287
288 // use this if your mesh is partitioned and you run code on parts,
289 // you can solve very big problems
290 CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL, "-my_is_partitioned",
291 &is_partitioned, &flg);
292
293 CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL, "-is_linear", &linear,
294 PETSC_NULL);
295
296 enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, LASBASETOP };
297 const char *list_bases[] = {"legendre", "lobatto", "bernstein_bezier"};
298 PetscInt choice_base_value = BERNSTEIN_BEZIER;
299 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
300 LASBASETOP, &choice_base_value, PETSC_NULL);
301 if (choice_base_value == LEGENDRE)
303 else if (choice_base_value == LOBATTO)
305 else if (choice_base_value == BERNSTEIN_BEZIER)
307
308 CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_disp_order",
309 &disp_order, &flg);
310 if (flg != PETSC_TRUE)
311 disp_order = 1;
312
313 CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_vel_order",
314 &vel_order, &flg);
315 if (flg != PETSC_TRUE)
316 vel_order = disp_order;
317
318 CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL,
319 "-my_solve_at_time_zero",
320 &is_solve_at_time_zero, &flg);
321
323 };
324
325 auto read_mesh = [&]() {
327 if (is_partitioned == PETSC_TRUE) {
328 // Read mesh to MOAB
329 const char *option;
330 option = "PARALLEL=BCAST_DELETE;"
331 "PARALLEL_RESOLVE_SHARED_ENTS;"
332 "PARTITION=PARALLEL_PARTITION;";
333 CHKERR moab.load_file(mesh_file_name, 0, option);
334 } else {
335 const char *option;
336 option = "";
337 CHKERR moab.load_file(mesh_file_name, 0, option);
338 }
340 };
341
342 CHKERR read_command_line_parameters();
343 CHKERR read_mesh();
344
345 MoFEM::Core core(moab);
346 MoFEM::Interface &m_field = core;
347
348 // ref meshset ref level 0
349 BitRefLevel bit_level0;
350 bit_level0.set(0);
351 EntityHandle meshset_level0;
352 CHKERR moab.create_meshset(MESHSET_SET, meshset_level0);
353 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
354 0, 3, bit_level0);
355 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByRefLevel(
356 bit_level0, BitRefLevel().set(), meshset_level0);
357
358 // Fields
359 CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
360 3, MB_TAG_SPARSE, MF_ZERO);
361 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
362 CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 2);
363 CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 2);
364 CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
365 CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
366
367 bool check_if_spatial_field_exist = m_field.check_field("DISPLACEMENT");
368 CHKERR m_field.add_field("DISPLACEMENT", H1, base, 3, MB_TAG_SPARSE,
369 MF_ZERO);
370 // add entities (by tets) to the field
371 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "DISPLACEMENT");
372
373 // set app. order
374 CHKERR m_field.set_field_order(0, MBTET, "DISPLACEMENT", disp_order);
375 CHKERR m_field.set_field_order(0, MBTRI, "DISPLACEMENT", disp_order);
376 CHKERR m_field.set_field_order(0, MBEDGE, "DISPLACEMENT", disp_order);
378 CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", disp_order);
379 else
380 CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
381
382 // Add nodal force element
383 CHKERR MetaNeumannForces::addNeumannBCElements(m_field, "DISPLACEMENT");
384 CHKERR MetaEdgeForces::addElement(m_field, "DISPLACEMENT");
385 CHKERR MetaNodalForces::addElement(m_field, "DISPLACEMENT");
386 // Add fluid pressure finite elements
387 FluidPressure fluid_pressure_fe(m_field);
388 fluid_pressure_fe.addNeumannFluidPressureBCElements("DISPLACEMENT");
389 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fluid_pressure_fe.getLoopFe(),
390 false, false);
391 fluid_pressure_fe.setNeumannFluidPressureFiniteElementOperators(
392 "DISPLACEMENT", PETSC_NULL, false, true);
393
394 // Velocity
395 CHKERR m_field.add_field("VELOCITY", H1, base, 3, MB_TAG_SPARSE, MF_ZERO);
396 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "VELOCITY");
397
398 CHKERR m_field.set_field_order(0, MBTET, "VELOCITY", vel_order);
399 CHKERR m_field.set_field_order(0, MBTRI, "VELOCITY", vel_order);
400 CHKERR m_field.set_field_order(0, MBEDGE, "VELOCITY", vel_order);
402 CHKERR m_field.set_field_order(0, MBVERTEX, "VELOCITY", vel_order);
403 else
404 CHKERR m_field.set_field_order(0, MBVERTEX, "VELOCITY", 1);
405
406 CHKERR m_field.add_field("DOT_DISPLACEMENT", H1, base, 3, MB_TAG_SPARSE,
407 MF_ZERO);
408 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "DOT_DISPLACEMENT");
409 CHKERR m_field.set_field_order(0, MBTET, "DOT_DISPLACEMENT", disp_order);
410 CHKERR m_field.set_field_order(0, MBTRI, "DOT_DISPLACEMENT", disp_order);
411 CHKERR m_field.set_field_order(0, MBEDGE, "DOT_DISPLACEMENT", disp_order);
413 CHKERR m_field.set_field_order(0, MBVERTEX, "DOT_DISPLACEMENT",
414 disp_order);
415 else
416 CHKERR m_field.set_field_order(0, MBVERTEX, "DOT_DISPLACEMENT", 1);
417
418 CHKERR m_field.add_field("DOT_VELOCITY", H1, base, 3, MB_TAG_SPARSE,
419 MF_ZERO);
420 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "DOT_VELOCITY");
421 CHKERR m_field.set_field_order(0, MBTET, "DOT_VELOCITY", vel_order);
422 CHKERR m_field.set_field_order(0, MBTRI, "DOT_VELOCITY", vel_order);
423 CHKERR m_field.set_field_order(0, MBEDGE, "DOT_VELOCITY", vel_order);
425 CHKERR m_field.set_field_order(0, MBVERTEX, "DOT_VELOCITY", disp_order);
426 else
427 CHKERR m_field.set_field_order(0, MBVERTEX, "DOT_VELOCITY", 1);
428
429 // Set material model and mass element
430 NonlinearElasticElement elastic(m_field, 2);
431 ElasticMaterials elastic_materials(m_field);
432 CHKERR elastic_materials.setBlocks(elastic.setOfBlocks);
433 // NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI<adouble>
434 // st_venant_kirchhoff_material_adouble;
435 // NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI<double>
436 // st_venant_kirchhoff_material_double; CHKERR
437 // elastic.setBlocks(&st_venant_kirchhoff_material_double,&st_venant_kirchhoff_material_adouble);
438 CHKERR elastic.addElement("ELASTIC", "DISPLACEMENT");
439 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", elastic.getLoopFeRhs(), true, false,
440 false, false);
441 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", elastic.getLoopFeLhs(), true, false,
442 false, false);
443 CHKERR addHOOpsVol("MESH_NODE_POSITIONS", elastic.getLoopFeEnergy(), true,
444 false, false, false);
445 CHKERR elastic.setOperators("DISPLACEMENT", "MESH_NODE_POSITIONS", false,
446 true);
447
448 // set mass element
449 ConvectiveMassElement inertia(m_field, 1);
450 // CHKERR inertia.setBlocks();
451 CHKERR elastic_materials.setBlocks(inertia.setOfBlocks);
452 CHKERR inertia.addConvectiveMassElement("MASS_ELEMENT", "VELOCITY",
453 "DISPLACEMENT");
454 CHKERR inertia.addVelocityElement("VELOCITY_ELEMENT", "VELOCITY",
455 "DISPLACEMENT");
456
457 // Add possibility to load accelerogram
458 {
459 string name = "-my_accelerogram";
460 char time_file_name[255];
461 PetscBool flg;
462 CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, name.c_str(),
463 time_file_name, 255, &flg);
464 if (flg == PETSC_TRUE) {
465 inertia.methodsOp.push_back(new TimeAccelerogram(name));
466 }
467 }
468
469 // damper element
470 KelvinVoigtDamper damper(m_field);
471 CHKERR elastic_materials.setBlocks(damper.blockMaterialDataMap);
472 {
473 KelvinVoigtDamper::CommonData &common_data = damper.commonData;
474 common_data.spatialPositionName = "DISPLACEMENT";
475 common_data.spatialPositionNameDot = "DOT_DISPLACEMENT";
476 CHKERR m_field.add_finite_element("DAMPER", MF_ZERO);
478 "DISPLACEMENT");
480 "DISPLACEMENT");
482 "DISPLACEMENT");
483
484 if (m_field.check_field("MESH_NODE_POSITIONS")) {
486 "DAMPER", "MESH_NODE_POSITIONS");
487 }
488 std::map<int, KelvinVoigtDamper::BlockMaterialData>::iterator bit =
489 damper.blockMaterialDataMap.begin();
490 for (; bit != damper.blockMaterialDataMap.end(); bit++) {
491 bit->second.lInear = linear;
492 int id = bit->first;
493 KelvinVoigtDamper::BlockMaterialData &material_data = bit->second;
494 damper.constitutiveEquationMap.insert(
496 material_data));
497 CHKERR m_field.add_ents_to_finite_element_by_type(bit->second.tEts,
498 MBTET, "DAMPER");
499 }
500 CHKERR damper.setOperators(3);
501 }
502
503 MonitorPostProc post_proc(m_field, elastic.setOfBlocks,
504 elastic.getLoopFeEnergy(),
505 inertia.getLoopFeEnergy());
506
507 // elastic and mass element calculated in Kuu shell matrix problem. To
508 // calculate Mass element, velocity field is needed.
509 CHKERR m_field.modify_finite_element_add_field_data("ELASTIC", "VELOCITY");
511 "DOT_DISPLACEMENT");
513 "DOT_VELOCITY");
514
515 // build field
516 CHKERR m_field.build_fields();
517 // CHKERR m_field.list_dofs_by_field_name("DISPLACEMENT");
518
519 // 10 node tets
520 if (!check_if_spatial_field_exist) {
521 Projection10NodeCoordsOnField ent_method_material(m_field,
522 "MESH_NODE_POSITIONS");
523 CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
524 }
525
526 // build finite elements
528 // build adjacencies
529 CHKERR m_field.build_adjacencies(bit_level0);
530
531 // define problems
532 {
533 CHKERR m_field.add_problem("Kuu", MF_ZERO);
534 CHKERR m_field.modify_problem_add_finite_element("Kuu", "ELASTIC");
535 CHKERR m_field.modify_problem_add_finite_element("Kuu", "PRESSURE_FE");
536 CHKERR m_field.modify_problem_add_finite_element("Kuu", "FORCE_FE");
538 "FLUID_PRESSURE_FE");
539 CHKERR m_field.modify_problem_ref_level_add_bit("Kuu", bit_level0);
540
541 ProblemsManager *prb_mng_ptr;
542 CHKERR m_field.getInterface(prb_mng_ptr);
543 if (is_partitioned) {
544 CHKERR prb_mng_ptr->buildProblemOnDistributedMesh("Kuu", true);
545 CHKERR prb_mng_ptr->partitionFiniteElements("Kuu", true, 0,
546 pcomm->size());
547 } else {
548 CHKERR prb_mng_ptr->buildProblem("Kuu", true);
549 CHKERR prb_mng_ptr->partitionProblem("Kuu");
550 CHKERR prb_mng_ptr->partitionFiniteElements("Kuu");
551 }
552 CHKERR prb_mng_ptr->partitionGhostDofs("Kuu");
553 }
554
555 CHKERR m_field.add_problem("DYNAMICS", MF_ZERO);
556 // set finite elements for problems
557 CHKERR m_field.modify_problem_add_finite_element("DYNAMICS", "ELASTIC");
558 CHKERR m_field.modify_problem_add_finite_element("DYNAMICS", "DAMPER");
559 CHKERR m_field.modify_problem_add_finite_element("DYNAMICS", "PRESSURE_FE");
560 CHKERR m_field.modify_problem_add_finite_element("DYNAMICS", "FORCE_FE");
561 CHKERR m_field.modify_problem_add_finite_element("DYNAMICS",
562 "FLUID_PRESSURE_FE");
563 CHKERR m_field.modify_problem_add_finite_element("DYNAMICS",
564 "MASS_ELEMENT");
565 CHKERR m_field.modify_problem_add_finite_element("DYNAMICS",
566 "VELOCITY_ELEMENT");
567 // set refinement level for problem
568 CHKERR m_field.modify_problem_ref_level_add_bit("DYNAMICS", bit_level0);
569
570 ProblemsManager *prb_mng_ptr;
571 CHKERR m_field.getInterface(prb_mng_ptr);
572 if (is_partitioned) {
573 CHKERR prb_mng_ptr->buildProblemOnDistributedMesh("DYNAMICS", true);
574 CHKERR prb_mng_ptr->partitionFiniteElements("DYNAMICS", true, 0,
575 pcomm->size());
576 } else {
577 CHKERR prb_mng_ptr->buildProblem("DYNAMICS", true);
578 CHKERR prb_mng_ptr->partitionProblem("DYNAMICS");
579 CHKERR prb_mng_ptr->partitionFiniteElements("DYNAMICS");
580 }
581 CHKERR prb_mng_ptr->partitionGhostDofs("DYNAMICS");
582
583 Vec F;
584 CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("DYNAMICS", COL,
585 &F);
586 Vec D;
587 CHKERR VecDuplicate(F, &D);
588
589 // create tS
590 TS ts;
591 CHKERR TSCreate(PETSC_COMM_WORLD, &ts);
592 CHKERR TSSetType(ts, TSBEULER);
593
594 // shell matrix
598 ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("Kuu",
599 &shellAij_ctx->K);
600 CHKERR MatDuplicate(shellAij_ctx->K, MAT_DO_NOT_COPY_VALUES,
601 &shellAij_ctx->M);
602 CHKERR shellAij_ctx->iNit();
603 CHKERR m_field.getInterface<VecManager>()->vecScatterCreate(
604 D, "DYNAMICS", COL, shellAij_ctx->u, "Kuu", COL,
605 &shellAij_ctx->scatterU);
606 CHKERR m_field.getInterface<VecManager>()->vecScatterCreate(
607 D, "DYNAMICS", "VELOCITY", COL, shellAij_ctx->v, "Kuu", "DISPLACEMENT",
608 COL, &shellAij_ctx->scatterV);
609 Mat shell_Aij;
610 const Problem *problem_ptr;
611 CHKERR m_field.get_problem("DYNAMICS", &problem_ptr);
612 CHKERR MatCreateShell(
613 PETSC_COMM_WORLD, problem_ptr->getNbLocalDofsRow(),
614 problem_ptr->getNbLocalDofsCol(), problem_ptr->getNbDofsRow(),
615 problem_ptr->getNbDofsRow(), (void *)shellAij_ctx, &shell_Aij);
616 CHKERR MatShellSetOperation(shell_Aij, MATOP_MULT,
617 (void (*)(void))ConvectiveMassElement::MultOpA);
618 CHKERR MatShellSetOperation(
619 shell_Aij, MATOP_ZERO_ENTRIES,
621 // blocked problem
622 ConvectiveMassElement::ShellMatrixElement shell_matrix_element(m_field);
623 DirichletDisplacementBc shell_dirichlet_bc(
624 m_field, "DISPLACEMENT", shellAij_ctx->barK, PETSC_NULL, PETSC_NULL);
625 DirichletDisplacementBc my_dirichlet_bc(m_field, "DISPLACEMENT", PETSC_NULL,
626 D, F);
627 shell_matrix_element.problemName = "Kuu";
628 shell_matrix_element.shellMatCtx = shellAij_ctx;
629 shell_matrix_element.DirichletBcPtr = &shell_dirichlet_bc;
630 shell_matrix_element.loopK.push_back(
631 ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
632 "ELASTIC", &elastic.getLoopFeLhs()));
633 // damper
634 shell_matrix_element.loopK.push_back(
635 ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
636 "ELASTIC", &damper.feLhs));
637
638 CHKERR inertia.addHOOpsVol();
639 CHKERR inertia.setShellMatrixMassOperators("VELOCITY", "DISPLACEMENT",
640 "MESH_NODE_POSITIONS", linear);
641 // element name "ELASTIC" is used, therefore M matrix is assembled as K
642 // matrix. This is added to M is shell matrix. M matrix is a derivative of
643 // inertia forces over spatial velocities
644 shell_matrix_element.loopM.push_back(
645 ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
646 "ELASTIC", &inertia.getLoopFeMassLhs()));
647 // this calculate derivatives of inertia forces over spatial positions and
648 // add this to shell K matrix
649 shell_matrix_element.loopAuxM.push_back(
650 ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
651 "ELASTIC", &inertia.getLoopFeMassAuxLhs()));
652
653 // Element to calculate shell matrix residual
654 ConvectiveMassElement::ShellResidualElement shell_matrix_residual(m_field);
655 shell_matrix_residual.shellMatCtx = shellAij_ctx;
656
657 // surface pressure
658 boost::ptr_map<std::string, NeumannForcesSurface> surface_forces;
659 {
660 string fe_name_str = "FORCE_FE";
661 surface_forces.insert(fe_name_str, new NeumannForcesSurface(m_field));
662 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS",
663 surface_forces.at(fe_name_str).getLoopFe(), false,
664 false);
666 NODESET | FORCESET, it)) {
667 CHKERR surface_forces.at(fe_name_str)
668 .addForce("DISPLACEMENT", PETSC_NULL, it->getMeshsetId(), true);
669 surface_forces.at(fe_name_str)
670 .methodsOp.push_back(new TimeForceScale());
671 }
672 }
673
674 boost::ptr_map<std::string, NeumannForcesSurface> surface_pressure;
675 {
676 string fe_name_str = "PRESSURE_FE";
677 surface_pressure.insert(fe_name_str, new NeumannForcesSurface(m_field));
678 CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS",
679 surface_pressure.at(fe_name_str).getLoopFe(), false,
680 false);
682 m_field, SIDESET | PRESSURESET, it)) {
683 CHKERR surface_pressure.at(fe_name_str)
684 .addPressure("DISPLACEMENT", PETSC_NULL, it->getMeshsetId(), true);
685 surface_pressure.at(fe_name_str)
686 .methodsOp.push_back(new TimeForceScale());
687 }
688 }
689
690 // edge forces
691 boost::ptr_map<std::string, EdgeForce> edge_forces;
692 {
693 string fe_name_str = "FORCE_FE";
694 edge_forces.insert(fe_name_str, new EdgeForce(m_field));
696 NODESET | FORCESET, it)) {
697 CHKERR edge_forces.at(fe_name_str)
698 .addForce("DISPLACEMENT", PETSC_NULL, it->getMeshsetId(), true);
699 edge_forces.at(fe_name_str).methodsOp.push_back(new TimeForceScale());
700 }
701 }
702
703 // nodal forces
704 boost::ptr_map<std::string, NodalForce> nodal_forces;
705 {
706 string fe_name_str = "FORCE_FE";
707 nodal_forces.insert(fe_name_str, new NodalForce(m_field));
709 NODESET | FORCESET, it)) {
710 CHKERR nodal_forces.at(fe_name_str)
711 .addForce("DISPLACEMENT", F, it->getMeshsetId(), true);
712 nodal_forces.at(fe_name_str).methodsOp.push_back(new TimeForceScale());
713 }
714 }
715
716 MonitorRestart monitor_restart(m_field, ts);
718 m_field, ts, "VELOCITY", "DISPLACEMENT");
719
720 // TS
721 TsCtx ts_ctx(m_field, "DYNAMICS");
722
723 // right hand side
724 // preprocess
725 ts_ctx.getPreProcessIFunction().push_back(&update_and_control);
726 ts_ctx.getPreProcessIFunction().push_back(&my_dirichlet_bc);
727
728 // fe looops
729 auto &loops_to_do_Rhs = ts_ctx.getLoopsIFunction();
730
731 auto add_static_rhs = [&](auto &loops_to_do_Rhs) {
733 loops_to_do_Rhs.push_back(
734 PairNameFEMethodPtr("ELASTIC", &elastic.getLoopFeRhs()));
735 for (auto fit = surface_forces.begin(); fit != surface_forces.end();
736 fit++) {
737 loops_to_do_Rhs.push_back(
738 PairNameFEMethodPtr(fit->first, &fit->second->getLoopFe()));
739 }
740 for (auto fit = surface_pressure.begin(); fit != surface_pressure.end();
741 fit++) {
742 loops_to_do_Rhs.push_back(
743 PairNameFEMethodPtr(fit->first, &fit->second->getLoopFe()));
744 }
745 for (auto fit = edge_forces.begin(); fit != edge_forces.end(); fit++) {
746 loops_to_do_Rhs.push_back(
747 PairNameFEMethodPtr(fit->first, &fit->second->getLoopFe()));
748 }
749 for (auto fit = nodal_forces.begin(); fit != nodal_forces.end(); fit++) {
750 loops_to_do_Rhs.push_back(
751 PairNameFEMethodPtr(fit->first, &fit->second->getLoopFe()));
752 }
753 loops_to_do_Rhs.push_back(PairNameFEMethodPtr(
754 "FLUID_PRESSURE_FE", &fluid_pressure_fe.getLoopFe()));
756 };
757
758 CHKERR add_static_rhs(loops_to_do_Rhs);
759
760 loops_to_do_Rhs.push_back(PairNameFEMethodPtr("DAMPER", &damper.feRhs));
761 loops_to_do_Rhs.push_back(
762 PairNameFEMethodPtr("MASS_ELEMENT", &inertia.getLoopFeMassRhs()));
763
764 // preporcess
765 // calculate residual for velocities
766 ts_ctx.getPreProcessIFunction().push_back(&shell_matrix_residual);
767 // postprocess
768 ts_ctx.getPostProcessIFunction().push_back(&my_dirichlet_bc);
769
770 // left hand side
771 // preprocess
772 ts_ctx.getPreProcessIJacobian().push_back(&update_and_control);
773 ts_ctx.getPreProcessIJacobian().push_back(&shell_matrix_element);
774 ts_ctx.getPostProcessIJacobian().push_back(&update_and_control);
775 // monitor
776 TsCtx::FEMethodsSequence &loopsMonitor =
778 loopsMonitor.push_back(
779 TsCtx::PairNameFEMethodPtr("MASS_ELEMENT", &post_proc));
780 loopsMonitor.push_back(
781 TsCtx::PairNameFEMethodPtr("MASS_ELEMENT", &monitor_restart));
782
783 CHKERR TSSetIFunction(ts, F, TsSetIFunction, &ts_ctx);
784 CHKERR TSSetIJacobian(ts, shell_Aij, shell_Aij, TsSetIJacobian, &ts_ctx);
785
786 CHKERR TSMonitorSet(ts, TsMonitorSet, &ts_ctx, PETSC_NULL);
787
788 double ftime = 1;
789 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
790 CHKERR TSSetSolution(ts, D);
791 CHKERR TSSetFromOptions(ts);
792 // shell matrix pre-conditioner
793 SNES snes;
794 CHKERR TSGetSNES(ts, &snes);
795 // CHKERR SNESSetFromOptions(snes);
796 KSP ksp;
797 CHKERR SNESGetKSP(snes, &ksp);
798 CHKERR KSPSetFromOptions(ksp);
799 PC pc;
800 CHKERR KSPGetPC(ksp, &pc);
801 CHKERR PCSetType(pc, PCSHELL);
802 ConvectiveMassElement::PCShellCtx pc_shell_ctx(shell_Aij);
803 CHKERR PCShellSetContext(pc, (void *)&pc_shell_ctx);
806 CHKERR PCShellSetDestroy(pc, ConvectiveMassElement::PCShellDestroy);
807
808 CHKERR VecZeroEntries(D);
809 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
810 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
811 CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
812 "DYNAMICS", COL, D, INSERT_VALUES, SCATTER_REVERSE);
813
814 // Solve problem at time Zero
815 if (is_solve_at_time_zero) {
816
817 Mat Aij = shellAij_ctx->K;
818 Vec F;
819 CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("Kuu", COL, &F);
820 Vec D;
821 CHKERR VecDuplicate(F, &D);
822
823 // Set vector for Kuu problem from the mesh data
824 CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
825 "Kuu", COL, D, INSERT_VALUES, SCATTER_FORWARD);
826 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
827 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
828
829 SnesCtx snes_ctx(m_field, "Kuu");
830
831 SNES snes;
832 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
833 CHKERR SNESSetApplicationContext(snes, &snes_ctx);
834 CHKERR SNESSetFunction(snes, F, SnesRhs, &snes_ctx);
835 CHKERR SNESSetJacobian(snes, Aij, Aij, SnesMat, &snes_ctx);
836 CHKERR SNESSetFromOptions(snes);
837
838 DirichletDisplacementBc my_dirichlet_bc(m_field, "DISPLACEMENT",
839 PETSC_NULL, D, F);
840
841 SnesCtx::FEMethodsSequence &loops_to_do_Rhs =
842 snes_ctx.getComputeRhs();
843 snes_ctx.getPreProcComputeRhs().push_back(&my_dirichlet_bc);
844 fluid_pressure_fe.getLoopFe().ts_t = 0;
845 CHKERR add_static_rhs(loops_to_do_Rhs);
846 snes_ctx.getPostProcComputeRhs().push_back(&my_dirichlet_bc);
847
848 SnesCtx::FEMethodsSequence &loops_to_do_Mat =
849 snes_ctx.getSetOperators();
850 snes_ctx.getPreProcSetOperators().push_back(&my_dirichlet_bc);
851 loops_to_do_Mat.push_back(
852 SnesCtx::PairNameFEMethodPtr("ELASTIC", &elastic.getLoopFeLhs()));
853 snes_ctx.getPostProcSetOperators().push_back(&my_dirichlet_bc);
854
855 CHKERR m_field.getInterface<FieldBlas>()->fieldScale(0, "VELOCITY");
856 CHKERR m_field.getInterface<FieldBlas>()->fieldScale(0,
857 "DOT_DISPLACEMENT");
858 CHKERR m_field.getInterface<FieldBlas>()->fieldScale(0, "DOT_VELOCITY");
859
860 CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
861 "Kuu", COL, D, INSERT_VALUES, SCATTER_FORWARD);
862
863 CHKERR SNESSolve(snes, PETSC_NULL, D);
864 int its;
865 CHKERR SNESGetIterationNumber(snes, &its);
866 MOFEM_LOG_C("DYNAMIC", Sev::inform, "number of Newton iterations = %d\n",
867 its);
868
869 // Set data on the mesh
870 CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
871 "Kuu", COL, D, INSERT_VALUES, SCATTER_REVERSE);
872
873 CHKERR VecDestroy(&F);
874 CHKERR VecDestroy(&D);
875 CHKERR SNESDestroy(&snes);
876 }
877
878 if (is_solve_at_time_zero) {
879 CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
880 "DYNAMICS", COL, D, INSERT_VALUES, SCATTER_FORWARD);
881 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
882 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
883 CHKERR TSSetSolution(ts, D);
884 }
885
886#if PETSC_VERSION_GE(3, 7, 0)
887 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
888#endif
889 CHKERR TSSolve(ts, D);
890 CHKERR TSGetTime(ts, &ftime);
891
892 PetscInt steps, snesfails, rejects, nonlinits, linits;
893 CHKERR TSGetTimeStepNumber(ts, &steps);
894 CHKERR TSGetSNESFailures(ts, &snesfails);
895 CHKERR TSGetStepRejections(ts, &rejects);
896 CHKERR TSGetSNESIterations(ts, &nonlinits);
897 CHKERR TSGetKSPIterations(ts, &linits);
898 MOFEM_LOG_C("DYNAMIC", Sev::inform,
899 "steps %d (%d rejected, %D SNES fails), ftime %g, nonlinits "
900 "%d, linits %D\n",
901 steps, rejects, snesfails, ftime, nonlinits, linits);
902 CHKERR TSDestroy(&ts);
903
904 CHKERR VecDestroy(&F);
905 CHKERR VecDestroy(&D);
906 CHKERR MatDestroy(&shellAij_ctx->K);
907 CHKERR MatDestroy(&shellAij_ctx->M);
908 CHKERR VecScatterDestroy(&shellAij_ctx->scatterU);
909 CHKERR VecScatterDestroy(&shellAij_ctx->scatterV);
910 CHKERR MatDestroy(&shell_Aij);
911 delete shellAij_ctx;
912 }
914
916
917 return 0;
918}
#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 nme:nme847.
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
char mesh_file_name[255]
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:165
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition TsCtx.cpp:259
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:139
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:27
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
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)
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
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.
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: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.
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:13
Interface for Time Stepping (TS) solver.
Definition TsCtx.hpp:17
BasicMethodsSequence & getPostProcessIJacobian()
Get the postProcess to do IJacobian object.
Definition TsCtx.hpp:128
MoFEM::FEMethodsSequence FEMethodsSequence
Definition TsCtx.hpp:26
FEMethodsSequence & getLoopsMonitor()
Get the loops to do Monitor object.
Definition TsCtx.hpp:98
BasicMethodsSequence & getPostProcessIFunction()
Get the postProcess to do IFunction object.
Definition TsCtx.hpp:112
FEMethodsSequence & getLoopsIFunction()
Get the loops to do IFunction object.
Definition TsCtx.hpp:59
BasicMethodsSequence & getPreProcessIFunction()
Get the preProcess to do IFunction object.
Definition TsCtx.hpp:105
BasicMethodsSequence & getPreProcessIJacobian()
Get the preProcess to do IJacobian object.
Definition TsCtx.hpp:121
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.