v0.14.0
Loading...
Searching...
No Matches
Classes | Namespaces | Functions | Variables
cell_forces.cpp File Reference

Calculate cell forces. More...

#include <BasicFiniteElements.hpp>
#include <Hooke.hpp>
#include <CellForces.hpp>
#include <boost/program_options.hpp>

Go to the source code of this file.

Classes

struct  CellEngineering::BlockOptionData
 

Namespaces

namespace  CellEngineering
 

Functions

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

Variables

static char help []
 
static int debug = 0
 

Detailed Description

Calculate cell forces.

Definition in file cell_forces.cpp.

Function Documentation

◆ main()

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

Definition at line 432 of file cell_forces.cpp.

432 {
433
434 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
435
436 try {
437
438 moab::Core mb_instance;
439 moab::Interface &moab = mb_instance;
440
441 PetscBool flg_block_config, flg_file;
442 char mesh_file_name[255];
443 char block_config_file[255];
444 PetscBool flg_order_force;
445 PetscInt order = 2;
446 PetscInt order_force = 2;
447 PetscBool flg_eps_u, flg_eps_rho, flg_eps_l;
448 double eps_u = 1e-6;
449 double eps_rho = 1e-3;
450 double eps_l = 0;
451 PetscBool is_curl = PETSC_TRUE;
452 CHKERR PetscOptionsBegin(PETSC_COMM_WORLD, "", "Elastic Config", "none");
453 CHKERR PetscOptionsString("-my_file", "mesh file name", "", "mesh.h5m",
454 mesh_file_name, 255, &flg_file);
455 CHKERR PetscOptionsInt("-my_order", "default approximation order", "",
456 order, &order, PETSC_NULL);
457 CHKERR PetscOptionsInt(
458 "-my_order_force",
459 "default approximation order for force approximation", "", order_force,
460 &order_force, &flg_order_force);
461 CHKERR PetscOptionsString("-my_block_config", "elastic configure file name",
462 "", "block_conf.in", block_config_file, 255,
463 &flg_block_config);
464 CHKERR PetscOptionsReal("-my_eps_rho", "foce optimality parameter ", "",
465 eps_rho, &eps_rho, &flg_eps_rho);
466 CHKERR PetscOptionsReal("-my_eps_u", "displacement optimality parameter ",
467 "", eps_u, &eps_u, &flg_eps_u);
468 CHKERR PetscOptionsReal("-my_eps_l", "displacement optimality parameter ",
469 "", eps_l, &eps_l, &flg_eps_l);
470 CHKERR PetscOptionsBool("-my_curl", "use Hcurl space to approximate forces",
471 "", is_curl, &is_curl, NULL);
472 ierr = PetscOptionsEnd();
473 CHKERRQ(ierr);
474
475 // Reade parameters from line command
476 if (flg_file != PETSC_TRUE) {
477 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
478 "*** ERROR -my_file (MESH FILE NEEDED)");
479 }
480
481 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
482 if (pcomm == NULL)
483 pcomm = new ParallelComm(&moab, PETSC_COMM_WORLD);
484
485 // Read mesh to MOAB
486 const char *option;
487 option = "PARALLEL=READ_PART;"
488 "PARALLEL_RESOLVE_SHARED_ENTS;"
489 "PARTITION=PARALLEL_PARTITION;";
490 // option = "";
491 CHKERR moab.load_file(mesh_file_name, 0, option);
492
493 // Create MoFEM (Joseph) database
494 MoFEM::Core core(moab);
495 MoFEM::Interface &m_field = core;
496
497 MeshsetsManager *mmanager_ptr;
498 CHKERR m_field.query_interface(mmanager_ptr);
499 // print bcs
500 CHKERR mmanager_ptr->printDisplacementSet();
501 CHKERR mmanager_ptr->printForceSet();
502 // print block sets with materials
503 CHKERR mmanager_ptr->printMaterialsSet();
504
505 // stl::bitset see for more details
506 BitRefLevel bit_level0;
507 bit_level0.set(0);
508 {
509 Range ents3d;
510 CHKERR moab.get_entities_by_dimension(0, 3, ents3d);
511 CHKERR m_field.seed_ref_level(ents3d, bit_level0, false);
512 }
513
514 // set app. order
515
516 std::vector<Range> setOrderToEnts(10);
517
518 // configure blocks by parsing config file
519 // it allow to set approximation order for each block independently
520 Range set_order_ents;
521 std::map<int, BlockOptionData> block_data;
522 if (flg_block_config) {
523 double read_eps_u, read_eps_rho, read_eps_l;
524 try {
525 ifstream ini_file(block_config_file);
526 if (!ini_file.is_open()) {
527 SETERRQ(PETSC_COMM_SELF, 1,
528 "*** -my_block_config does not exist ***");
529 }
530 // std::cerr << block_config_file << std::endl;
531 po::variables_map vm;
532 po::options_description config_file_options;
533 config_file_options.add_options()(
534 "eps_u", po::value<double>(&read_eps_u)->default_value(-1))(
535 "eps_rho", po::value<double>(&read_eps_rho)->default_value(-1))(
536 "eps_l", po::value<double>(&read_eps_l)->default_value(-1));
537 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
538 std::ostringstream str_order;
539 str_order << "block_" << it->getMeshsetId() << ".displacement_order";
540 config_file_options.add_options()(
541 str_order.str().c_str(),
542 po::value<int>(&block_data[it->getMeshsetId()].oRder)
543 ->default_value(order));
544 std::ostringstream str_cond;
545 str_cond << "block_" << it->getMeshsetId() << ".young_modulus";
546 config_file_options.add_options()(
547 str_cond.str().c_str(),
548 po::value<double>(&block_data[it->getMeshsetId()].yOung)
549 ->default_value(-1));
550 std::ostringstream str_capa;
551 str_capa << "block_" << it->getMeshsetId() << ".poisson_ratio";
552 config_file_options.add_options()(
553 str_capa.str().c_str(),
554 po::value<double>(&block_data[it->getMeshsetId()].pOisson)
555 ->default_value(-2));
556 }
557 po::parsed_options parsed =
558 parse_config_file(ini_file, config_file_options, true);
559 store(parsed, vm);
560 po::notify(vm);
561 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
562 if (block_data[it->getMeshsetId()].oRder == -1)
563 continue;
564 if (block_data[it->getMeshsetId()].oRder == order)
565 continue;
566 PetscPrintf(PETSC_COMM_WORLD, "Set block %d order to %d\n",
567 it->getMeshsetId(), block_data[it->getMeshsetId()].oRder);
568 Range block_ents;
569 CHKERR moab.get_entities_by_handle(it->getMeshset(), block_ents,
570 true);
571 // block_ents = block_ents.subset_by_type(MBTET);
572 Range nodes;
573 CHKERR moab.get_connectivity(block_ents, nodes, true);
574 Range ents_to_set_order, ents3d;
575 CHKERR moab.get_adjacencies(nodes, 3, false, ents3d,
576 moab::Interface::UNION);
577 CHKERR moab.get_adjacencies(ents3d, 2, false, ents_to_set_order,
578 moab::Interface::UNION);
579 CHKERR moab.get_adjacencies(ents3d, 1, false, ents_to_set_order,
580 moab::Interface::UNION);
581 ents_to_set_order = subtract(
582 ents_to_set_order, ents_to_set_order.subset_by_type(MBQUAD));
583 ents_to_set_order = subtract(
584 ents_to_set_order, ents_to_set_order.subset_by_type(MBPRISM));
585 set_order_ents.merge(ents3d);
586 set_order_ents.merge(ents_to_set_order);
587 setOrderToEnts[block_data[it->getMeshsetId()].oRder].merge(
588 set_order_ents);
589 }
590 CHKERR m_field.synchronise_entities(set_order_ents, 0);
591 std::vector<std::string> additional_parameters;
592 additional_parameters =
593 collect_unrecognized(parsed.options, po::include_positional);
594 for (std::vector<std::string>::iterator vit =
595 additional_parameters.begin();
596 vit != additional_parameters.end(); vit++) {
597 CHKERR PetscPrintf(PETSC_COMM_WORLD,
598 "** WARNING Unrecognized option %s\n",
599 vit->c_str());
600 }
601 } catch (const std::exception &ex) {
602 std::ostringstream ss;
603 ss << ex.what() << std::endl;
604 SETERRQ(PETSC_COMM_SELF, MOFEM_STD_EXCEPTION_THROW, ss.str().c_str());
605 }
606 if (read_eps_u > 0) {
607 eps_u = read_eps_u;
608 };
609 if (read_eps_rho > 0) {
610 eps_rho = read_eps_rho;
611 }
612 if (read_eps_l > 0) {
613 eps_l = read_eps_l;
614 }
615 }
616
617 PetscPrintf(PETSC_COMM_WORLD, "epsU = %6.4e epsRho = %6.4e\n", eps_u,
618 eps_rho);
619
620 // Fields
621 CHKERR m_field.add_field("U", H1, AINSWORTH_LEGENDRE_BASE, 3, MB_TAG_SPARSE,
622 MF_ZERO);
623 CHKERR m_field.add_field("UPSILON", H1, AINSWORTH_LEGENDRE_BASE, 3,
624 MB_TAG_SPARSE, MF_ZERO);
625 if (is_curl) {
626 CHKERR m_field.add_field("RHO", HCURL, AINSWORTH_LEGENDRE_BASE, 1);
627 } else {
628 CHKERR m_field.add_field("RHO", H1, AINSWORTH_LEGENDRE_BASE, 1);
629 }
630
631 // add entitities (by tets) to the field
632 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "U");
633 CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "UPSILON");
634 CHKERR m_field.add_ents_to_field_by_type(0, MBPRISM, "U");
635 CHKERR m_field.add_ents_to_field_by_type(0, MBPRISM, "UPSILON");
636
637 CHKERR m_field.synchronise_field_entities("U");
638 CHKERR m_field.synchronise_field_entities("UPSILON");
639 CHKERR m_field.synchronise_field_entities("RHO");
640
641 Range vertex_to_fix;
642 Range edges_to_fix;
643 Range ents_1st_layer;
644 // If problem is partitioned meshset culd not exist on this part
645 if (mmanager_ptr->checkMeshset(202, SIDESET)) {
646 CHKERR mmanager_ptr->getEntitiesByDimension(202, SIDESET, 2,
647 ents_1st_layer, true);
648 CHKERR mmanager_ptr->getEntitiesByDimension(202, SIDESET, 0,
649 vertex_to_fix, false);
650 CHKERR mmanager_ptr->getEntitiesByDimension(202, SIDESET, 1, edges_to_fix,
651 false);
652 if (vertex_to_fix.size() != 1 && !vertex_to_fix.empty()) {
653 SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
654 "Should be one vertex only, but is %d", vertex_to_fix.size());
655 }
656 }
657 CHKERR m_field.synchronise_entities(ents_1st_layer, 0);
658 CHKERR m_field.add_ents_to_field_by_type(
659 ents_1st_layer.subset_by_type(MBTRI), MBTRI, "RHO");
660 Range ents_2nd_layer;
661 // If problem is partitioned meshset culd not exist on this part
662 if (mmanager_ptr->checkMeshset(101, SIDESET)) {
663 CHKERR mmanager_ptr->getEntitiesByDimension(101, SIDESET, 2,
664 ents_2nd_layer, true);
665 }
666 CHKERR m_field.synchronise_entities(ents_2nd_layer, 0);
667
668 for (int oo = 2; oo != setOrderToEnts.size(); oo++) {
669 if (setOrderToEnts[oo].size() > 0) {
670 CHKERR m_field.synchronise_entities(setOrderToEnts[oo], 0);
671 CHKERR m_field.set_field_order(setOrderToEnts[oo], "U", oo);
672 CHKERR m_field.set_field_order(setOrderToEnts[oo], "UPSILON", oo);
673 }
674 }
675
676 const int through_thickness_order = 2;
677 {
678 Range ents3d;
679 CHKERR moab.get_entities_by_dimension(0, 3, ents3d);
680 Range ents;
681 CHKERR moab.get_adjacencies(ents3d, 2, false, ents,
682 moab::Interface::UNION);
683 CHKERR moab.get_adjacencies(ents3d, 1, false, ents,
684 moab::Interface::UNION);
685
686 Range prisms;
687 CHKERR moab.get_entities_by_type(0, MBPRISM, prisms);
688 {
689 Range quads;
690 CHKERR moab.get_adjacencies(prisms, 2, false, quads,
691 moab::Interface::UNION);
692 Range prism_tris;
693 prism_tris = quads.subset_by_type(MBTRI);
694 quads = subtract(quads, prism_tris);
695 Range quads_edges;
696 CHKERR moab.get_adjacencies(quads, 1, false, quads_edges,
697 moab::Interface::UNION);
698 Range prism_tris_edges;
699 CHKERR moab.get_adjacencies(prism_tris, 1, false, prism_tris_edges,
700 moab::Interface::UNION);
701 quads_edges = subtract(quads_edges, prism_tris_edges);
702 prisms.merge(quads);
703 prisms.merge(quads_edges);
704 }
705
706 ents.merge(ents3d);
707 ents = subtract(ents, set_order_ents);
708 ents = subtract(ents, prisms);
709
710 CHKERR m_field.synchronise_entities(ents, 0);
711 CHKERR m_field.synchronise_entities(prisms, 0);
712
713 CHKERR m_field.set_field_order(ents, "U", order);
714 CHKERR m_field.set_field_order(ents, "UPSILON", order);
715 // approx. order through thickness to 2
716 CHKERR m_field.set_field_order(prisms, "U", through_thickness_order);
717 CHKERR m_field.set_field_order(prisms, "UPSILON",
718 through_thickness_order);
719 }
720 CHKERR m_field.set_field_order(0, MBVERTEX, "U", 1);
721 CHKERR m_field.set_field_order(0, MBVERTEX, "UPSILON", 1);
722
723 if (is_curl) {
724 CHKERR m_field.set_field_order(0, MBTRI, "RHO", order_force);
725 CHKERR m_field.set_field_order(0, MBEDGE, "RHO", order_force);
726 } else {
727 CHKERR m_field.set_field_order(0, MBTRI, "RHO", order_force);
728 CHKERR m_field.set_field_order(0, MBEDGE, "RHO", order_force);
729 CHKERR m_field.set_field_order(0, MBVERTEX, "RHO", 1);
730 }
731
732 // Add elastic element
733 boost::shared_ptr<Hooke<adouble>> hooke_adouble_ptr(new Hooke<adouble>());
734 boost::shared_ptr<Hooke<double>> hooke_double_ptr(new Hooke<double>());
735 NonlinearElasticElement elastic(m_field, 2);
736 CHKERR elastic.setBlocks(hooke_double_ptr, hooke_adouble_ptr);
737 CHKERR elastic.addElement("ELASTIC", "U");
738 CHKERR elastic.setOperators("U", "MESH_NODE_POSITIONS", false, true);
739 // Add prisms
740 CellEngineering::FatPrism fat_prism_rhs(m_field);
741 CellEngineering::FatPrism fat_prism_lhs(m_field);
742 {
743 CHKERR m_field.add_finite_element("ELASTIC_PRISM", MF_ZERO);
744 CHKERR m_field.modify_finite_element_add_field_row("ELASTIC_PRISM", "U");
745 CHKERR m_field.modify_finite_element_add_field_col("ELASTIC_PRISM", "U");
746 CHKERR m_field.modify_finite_element_add_field_data("ELASTIC_PRISM", "U");
747 Range ents_2nd_layer;
748 CHKERR mmanager_ptr->getEntitiesByDimension(2, BLOCKSET, 3,
749 elastic.setOfBlocks[2].tEts);
750 CHKERR m_field.add_ents_to_finite_element_by_type(
751 elastic.setOfBlocks[2].tEts, MBPRISM, "ELASTIC_PRISM");
752 MatrixDouble inv_jac;
753 // right hand side operators
754 fat_prism_rhs.getOpPtrVector().push_back(
755 new MoFEM::OpCalculateInvJacForFatPrism(inv_jac));
756 fat_prism_rhs.getOpPtrVector().push_back(
757 new MoFEM::OpSetInvJacH1ForFatPrism(inv_jac));
758 fat_prism_rhs.getOpPtrVector().push_back(
759 new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
760 "U", elastic.commonData));
761 fat_prism_rhs.getOpPtrVector().push_back(
762 new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
763 "U", elastic.setOfBlocks[2], elastic.commonData, 2, false, false,
764 true));
765 fat_prism_rhs.getOpPtrVector().push_back(
766 new NonlinearElasticElement::OpRhsPiolaKirchhoff(
767 "U", elastic.setOfBlocks[2], elastic.commonData));
768 // Left hand side operators
769 fat_prism_lhs.getOpPtrVector().push_back(
770 new MoFEM::OpCalculateInvJacForFatPrism(inv_jac));
771 fat_prism_lhs.getOpPtrVector().push_back(
772 new MoFEM::OpSetInvJacH1ForFatPrism(inv_jac));
773 fat_prism_lhs.getOpPtrVector().push_back(
774 new NonlinearElasticElement::OpGetCommonDataAtGaussPts(
775 "U", elastic.commonData));
776 fat_prism_lhs.getOpPtrVector().push_back(
777 new NonlinearElasticElement::OpJacobianPiolaKirchhoffStress(
778 "U", elastic.setOfBlocks[2], elastic.commonData, 2, true, false,
779 true));
780 fat_prism_lhs.getOpPtrVector().push_back(
781 new NonlinearElasticElement::OpLhsPiolaKirchhoff_dx(
782 "U", "U", elastic.setOfBlocks[2], elastic.commonData));
783 }
784
785 // Update material parameters
786 for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
787 int id = it->getMeshsetId();
788 if (block_data[id].yOung > 0) {
789 elastic.setOfBlocks[id].E = block_data[id].yOung;
790 CHKERR PetscPrintf(PETSC_COMM_WORLD,
791 "Block %d set Young modulus %3.4g\n", id,
792 elastic.setOfBlocks[id].E);
793 }
794 if (block_data[id].pOisson >= -1) {
795 elastic.setOfBlocks[id].PoissonRatio = block_data[id].pOisson;
796 CHKERR PetscPrintf(PETSC_COMM_WORLD,
797 "Block %d set Poisson ratio %3.4g\n", id,
798 elastic.setOfBlocks[id].PoissonRatio);
799 }
800 }
801
802 // build field
803 CHKERR m_field.build_fields();
804
805 // Control elements
806 CHKERR m_field.add_finite_element("KUPSUPS");
807 CHKERR m_field.modify_finite_element_add_field_row("KUPSUPS", "UPSILON");
808 CHKERR m_field.modify_finite_element_add_field_col("KUPSUPS", "UPSILON");
809 CHKERR m_field.modify_finite_element_add_field_data("KUPSUPS", "UPSILON");
810 CHKERR m_field.add_ents_to_finite_element_by_type(0, MBTET, "KUPSUPS");
811 CHKERR m_field.add_ents_to_finite_element_by_type(0, MBPRISM, "KUPSUPS");
812
813 CHKERR m_field.add_finite_element("DISPLACEMENTS_PENALTY");
814 CHKERR m_field.modify_finite_element_add_field_row("DISPLACEMENTS_PENALTY",
815 "UPSILON");
816 CHKERR m_field.modify_finite_element_add_field_col("DISPLACEMENTS_PENALTY",
817 "U");
818 CHKERR m_field.modify_finite_element_add_field_data("DISPLACEMENTS_PENALTY",
819 "UPSILON");
820 CHKERR m_field.modify_finite_element_add_field_data("DISPLACEMENTS_PENALTY",
821 "U");
822 CHKERR m_field.modify_finite_element_add_field_data("DISPLACEMENTS_PENALTY",
823 "DISP_X");
824 CHKERR m_field.modify_finite_element_add_field_data("DISPLACEMENTS_PENALTY",
825 "DISP_Y");
826 CHKERR m_field.add_ents_to_finite_element_by_type(ents_2nd_layer, MBTRI,
827 "DISPLACEMENTS_PENALTY");
828
829 // Add element to calculate residual on 1st layer
830 CHKERR m_field.add_finite_element("BT");
831 CHKERR m_field.modify_finite_element_add_field_row("BT", "U");
832 CHKERR m_field.modify_finite_element_add_field_row("BT", "UPSILON");
833 CHKERR m_field.modify_finite_element_add_field_col("BT", "RHO");
834 CHKERR m_field.modify_finite_element_add_field_data("BT", "U");
835 CHKERR m_field.modify_finite_element_add_field_data("BT", "UPSILON");
836 CHKERR m_field.modify_finite_element_add_field_data("BT", "RHO");
837 CHKERR m_field.add_ents_to_finite_element_by_type(ents_1st_layer, MBTRI,
838 "BT");
839
840 CHKERR m_field.add_finite_element("B");
841 CHKERR m_field.modify_finite_element_add_field_row("B", "RHO");
842 CHKERR m_field.modify_finite_element_add_field_col("B", "U");
843 CHKERR m_field.modify_finite_element_add_field_col("B", "UPSILON");
844 CHKERR m_field.modify_finite_element_add_field_data("B", "U");
845 CHKERR m_field.modify_finite_element_add_field_data("B", "UPSILON");
846 CHKERR m_field.modify_finite_element_add_field_data("B", "RHO");
847 CHKERR m_field.add_ents_to_finite_element_by_type(ents_1st_layer, MBTRI,
848 "B");
849
850 // Add element to calculate residual on 1st layer
851 CHKERR m_field.add_finite_element("D");
852 CHKERR m_field.modify_finite_element_add_field_row("D", "RHO");
853 CHKERR m_field.modify_finite_element_add_field_col("D", "RHO");
854 CHKERR m_field.modify_finite_element_add_field_data("D", "RHO");
855 CHKERR m_field.add_ents_to_finite_element_by_type(ents_1st_layer, MBTRI,
856 "D");
857
858 // build finite elemnts
859 CHKERR m_field.build_finite_elements();
860 // build adjacencies
861 CHKERR m_field.build_adjacencies(bit_level0);
862
863 // Register MOFEM DM
864 DMType dm_name = "MOFEM";
865 CHKERR DMRegister_MoFEM(dm_name);
866
867 DM dm_control;
868 {
869 // Craete dm_control instance
870 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_control);
871 CHKERR DMSetType(dm_control, dm_name);
872 // set dm_control data structure which created mofem datastructures
873 CHKERR DMMoFEMCreateMoFEM(dm_control, &m_field, "CONTROL_PROB",
874 bit_level0);
875 CHKERR DMSetFromOptions(dm_control);
876 CHKERR DMMoFEMSetSquareProblem(dm_control, PETSC_TRUE);
877 CHKERR DMMoFEMSetIsPartitioned(dm_control, PETSC_TRUE);
878 // add elements to dm_control
879 CHKERR DMMoFEMAddElement(dm_control, "ELASTIC");
880 CHKERR DMMoFEMAddElement(dm_control, "ELASTIC_PRISM");
881 CHKERR DMMoFEMAddElement(dm_control, "KUPSUPS");
882 CHKERR DMMoFEMAddElement(dm_control, "DISPLACEMENTS_PENALTY");
883 CHKERR DMMoFEMAddElement(dm_control, "B");
884 CHKERR DMMoFEMAddElement(dm_control, "BT");
885 CHKERR DMMoFEMAddElement(dm_control, "D");
886 CHKERR DMSetUp(dm_control);
887 }
888
889 MatrixDouble inv_jac;
890 CellEngineering::CommonData common_data;
891
892 ublas::matrix<Mat> nested_matrices(2, 2);
893 ublas::vector<IS> nested_is_rows(2);
894 ublas::vector<IS> nested_is_cols(2);
895 for (int i = 0; i != 2; i++) {
896 nested_is_rows[i] = PETSC_NULL;
897 nested_is_cols[i] = PETSC_NULL;
898 for (int j = 0; j != 2; j++) {
899 nested_matrices(i, j) = PETSC_NULL;
900 }
901 }
902
903 ublas::matrix<Mat> sub_nested_matrices(2, 2);
904 ublas::vector<IS> sub_nested_is_rows(2);
905 ublas::vector<IS> sub_nested_is_cols(2);
906 for (int i = 0; i != 2; i++) {
907 sub_nested_is_rows[i] = PETSC_NULL;
908 sub_nested_is_cols[i] = PETSC_NULL;
909 for (int j = 0; j != 2; j++) {
910 sub_nested_matrices(i, j) = PETSC_NULL;
911 }
912 }
913
914 DM dm_sub_volume_control;
915 {
916 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_volume_control);
917 CHKERR DMSetType(dm_sub_volume_control, dm_name);
918 // set dm_sub_volume_control data structure which created mofem data
919 // structures
920 CHKERR DMMoFEMCreateSubDM(dm_sub_volume_control, dm_control,
921 "SUB_CONTROL_PROB");
922 CHKERR DMMoFEMSetSquareProblem(dm_sub_volume_control, PETSC_TRUE);
923 CHKERR DMMoFEMAddSubFieldRow(dm_sub_volume_control, "U");
924 CHKERR DMMoFEMAddSubFieldRow(dm_sub_volume_control, "UPSILON");
925 CHKERR DMMoFEMAddSubFieldCol(dm_sub_volume_control, "U");
926 CHKERR DMMoFEMAddSubFieldCol(dm_sub_volume_control, "UPSILON");
927 // add elements to dm_sub_volume_control
928 CHKERR DMSetUp(dm_sub_volume_control);
929
930 const Problem *prb_ptr;
931 CHKERR m_field.get_problem("SUB_CONTROL_PROB", &prb_ptr);
932 boost::shared_ptr<Problem::SubProblemData> sub_data =
933 prb_ptr->getSubData();
934
935 CHKERR sub_data->getRowIs(&nested_is_rows[0]);
936 CHKERR sub_data->getColIs(&nested_is_cols[0]);
937 // That will be filled at the end
938 nested_matrices(0, 0) = PETSC_NULL;
939 }
940
941 {
942 DM dm_sub_sub_elastic;
943
944 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_sub_elastic);
945 CHKERR DMSetType(dm_sub_sub_elastic, dm_name);
946 // set dm_sub_sub_elastic data structure which created mofem data
947 // structures
948 CHKERR DMMoFEMCreateSubDM(dm_sub_sub_elastic, dm_sub_volume_control,
949 "ELASTIC_PROB");
950 CHKERR DMMoFEMSetSquareProblem(dm_sub_sub_elastic, PETSC_TRUE);
951 CHKERR DMMoFEMAddElement(dm_sub_sub_elastic, "ELASTIC");
952 CHKERR DMMoFEMAddElement(dm_sub_sub_elastic, "ELASTIC_PRISM");
953 CHKERR DMMoFEMAddSubFieldRow(dm_sub_sub_elastic, "U");
954 CHKERR DMMoFEMAddSubFieldCol(dm_sub_sub_elastic, "U");
955 // add elements to dm_sub_sub_elastic
956 CHKERR DMSetUp(dm_sub_sub_elastic);
957
958 Mat Kuu;
959 Vec Du, Fu;
960 CHKERR DMCreateMatrix(dm_sub_sub_elastic, &Kuu);
961 CHKERR DMCreateGlobalVector(dm_sub_sub_elastic, &Du);
962 CHKERR DMCreateGlobalVector(dm_sub_sub_elastic, &Fu);
963 CHKERR MatZeroEntries(Kuu);
964 CHKERR VecZeroEntries(Du);
965 CHKERR VecZeroEntries(Fu);
966 CHKERR DMoFEMMeshToLocalVector(dm_sub_sub_elastic, Du, INSERT_VALUES,
967 SCATTER_REVERSE);
968
969 // Manage Dirichlet bC
970 boost::shared_ptr<FEMethod> dirihlet_bc_ptr;
971 dirihlet_bc_ptr = boost::shared_ptr<FEMethod>(
972 new DirichletDisplacementBc(m_field, "U", Kuu, Du, Fu));
973 dirihlet_bc_ptr->snes_ctx = FEMethod::CTX_SNESNONE;
974 dirihlet_bc_ptr->ts_ctx = FEMethod::CTX_TSNONE;
975 // preproc
976 CHKERR DMoFEMPreProcessFiniteElements(dm_sub_sub_elastic,
977 dirihlet_bc_ptr.get());
978 // internal force vector (to take into account Dirichlet boundary
979 // conditions)
980 elastic.getLoopFeRhs().snes_f = Fu;
981 fat_prism_rhs.snes_f = Fu;
982 CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC",
983 &elastic.getLoopFeRhs());
984 CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC_PRISM",
985 &fat_prism_rhs);
986 // elastic element matrix
987 elastic.getLoopFeLhs().snes_B = Kuu;
988 fat_prism_lhs.snes_B = Kuu;
989 CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC",
990 &elastic.getLoopFeLhs());
991 CHKERR DMoFEMLoopFiniteElements(dm_sub_sub_elastic, "ELASTIC_PRISM",
992 &fat_prism_lhs);
993 // postproc
994 CHKERR DMoFEMPostProcessFiniteElements(dm_sub_sub_elastic,
995 dirihlet_bc_ptr.get());
996 CHKERR VecGhostUpdateBegin(Fu, ADD_VALUES, SCATTER_REVERSE);
997 CHKERR VecGhostUpdateEnd(Fu, ADD_VALUES, SCATTER_REVERSE);
998 CHKERR VecAssemblyBegin(Fu);
999 CHKERR VecAssemblyEnd(Fu);
1000 CHKERR VecScale(Fu, -1);
1001 CHKERR VecDestroy(&Du);
1002 CHKERR VecDestroy(&Fu);
1003
1004 const Problem *prb_ptr;
1005 CHKERR m_field.get_problem("ELASTIC_PROB", &prb_ptr);
1006 boost::shared_ptr<Problem::SubProblemData> sub_data =
1007 prb_ptr->getSubData();
1008
1009 CHKERR sub_data->getRowIs(&sub_nested_is_rows[0]);
1010 CHKERR sub_data->getColIs(&sub_nested_is_cols[0]);
1011 sub_nested_matrices(0, 0) = Kuu;
1012 IS isUpsilon;
1013 CHKERR m_field.query_interface<ISManager>()
1014 ->isCreateFromProblemFieldToOtherProblemField(
1015 "ELASTIC_PROB", "U", ROW, "SUB_CONTROL_PROB", "UPSILON", ROW,
1016 PETSC_NULL, &isUpsilon);
1017 sub_nested_is_rows[1] = isUpsilon;
1018 sub_nested_is_cols[1] = isUpsilon;
1019 sub_nested_matrices(1, 1) = Kuu;
1020 PetscObjectReference((PetscObject)Kuu);
1021 PetscObjectReference((PetscObject)isUpsilon);
1022
1023 // Matrix View
1024 if (debug) {
1025 cerr << "Kuu" << endl;
1026 MatView(Kuu, PETSC_VIEWER_DRAW_WORLD);
1027 std::string wait;
1028 std::cin >> wait;
1029 }
1030
1031 CHKERR DMDestroy(&dm_sub_sub_elastic);
1032 }
1033
1034 {
1035 DM dm_sub_disp_penalty;
1036
1037 // Craete dm_control instance
1038 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_disp_penalty);
1039 CHKERR DMSetType(dm_sub_disp_penalty, dm_name);
1040 // set dm_sub_disp_penalty data structure which created mofem
1041 // datastructures
1042 CHKERR DMMoFEMCreateSubDM(dm_sub_disp_penalty, dm_sub_volume_control,
1043 "S_PROB");
1044 CHKERR DMMoFEMSetSquareProblem(dm_sub_disp_penalty, PETSC_FALSE);
1045 CHKERR DMMoFEMAddSubFieldRow(dm_sub_disp_penalty, "UPSILON");
1046 CHKERR DMMoFEMAddSubFieldCol(dm_sub_disp_penalty, "U");
1047 // add elements to dm_sub_disp_penalty
1048 CHKERR DMMoFEMAddElement(dm_sub_disp_penalty, "DISPLACEMENTS_PENALTY");
1049 CHKERR DMSetUp(dm_sub_disp_penalty);
1050
1051 Mat S;
1052 CHKERR DMCreateMatrix(dm_sub_disp_penalty, &S);
1053 CHKERR MatZeroEntries(S);
1054
1055 CellEngineering::FaceElement face_element(m_field);
1056 face_element.getOpPtrVector().push_back(new OpCellS(S, eps_u));
1057 CHKERR DMoFEMLoopFiniteElements(dm_sub_disp_penalty,
1058 "DISPLACEMENTS_PENALTY", &face_element);
1059 CHKERR MatAssemblyBegin(S, MAT_FLUSH_ASSEMBLY);
1060 CHKERR MatAssemblyEnd(S, MAT_FLUSH_ASSEMBLY);
1061
1062 // // Matrix View
1063 if (debug) {
1064 cerr << "S" << endl;
1065 MatView(S, PETSC_VIEWER_DRAW_WORLD);
1066 std::string wait;
1067 std::cin >> wait;
1068 }
1069
1070 const Problem *problem_ptr;
1071 CHKERR m_field.get_problem("S_PROB", &problem_ptr);
1072 boost::shared_ptr<Problem::SubProblemData> sub_data =
1073 problem_ptr->getSubData();
1074 // CHKERR sub_data->getRowIs(&sub_nested_is_rows[1]);
1075 // CHKERR sub_data->getColIs(&sub_nested_is_cols[0]);
1076 sub_nested_matrices(1, 0) = S;
1077
1078 CHKERR DMDestroy(&dm_sub_disp_penalty);
1079 }
1080
1081 // Calculate penalty matrix
1082 {
1083 DM dm_sub_force_penalty;
1084
1085 // Craete dm_control instance
1086 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_force_penalty);
1087 CHKERR DMSetType(dm_sub_force_penalty, dm_name);
1088 // set dm_sub_force_penalty data structure which created mofem
1089 // datastructures
1090 CHKERR DMMoFEMCreateSubDM(dm_sub_force_penalty, dm_control, "D_PROB");
1091 CHKERR DMMoFEMSetSquareProblem(dm_sub_force_penalty, PETSC_TRUE);
1092 CHKERR DMMoFEMAddSubFieldRow(dm_sub_force_penalty, "RHO");
1093 CHKERR DMMoFEMAddSubFieldCol(dm_sub_force_penalty, "RHO");
1094 // add elements to dm_sub_force_penalty
1095 CHKERR DMMoFEMAddElement(dm_sub_force_penalty, "D");
1096 CHKERR DMSetUp(dm_sub_force_penalty);
1097
1098 Mat D;
1099 CHKERR DMCreateMatrix(dm_sub_force_penalty, &D);
1100 CHKERR MatZeroEntries(D);
1101
1102 {
1103 CellEngineering::FaceElement face_d_matrix(m_field);
1104 face_d_matrix.getOpPtrVector().push_back(
1105 new OpCalculateInvJacForFace(inv_jac));
1106 if (is_curl) {
1107 face_d_matrix.getOpPtrVector().push_back(
1108 new OpSetInvJacHcurlFace(inv_jac));
1109 face_d_matrix.getOpPtrVector().push_back(
1110 new OpCellCurlD(D, eps_rho / eps_u, eps_l * eps_u));
1111 } else {
1112 face_d_matrix.getOpPtrVector().push_back(
1113 new OpSetInvJacH1ForFace(inv_jac));
1114 face_d_matrix.getOpPtrVector().push_back(
1115 new OpCellPotentialD(D, eps_rho / eps_u));
1116 }
1117 CHKERR DMoFEMLoopFiniteElements(dm_sub_force_penalty, "D",
1118 &face_d_matrix);
1119 }
1120 CHKERR MatAssemblyBegin(D, MAT_FINAL_ASSEMBLY);
1121 CHKERR MatAssemblyEnd(D, MAT_FINAL_ASSEMBLY);
1122
1123 const Problem *problem_ptr;
1124 CHKERR m_field.get_problem("D_PROB", &problem_ptr);
1125
1126 // Zero rows, force field is given by gradients of potential field, so one
1127 // of the values has to be fixed like for rigid body motion.
1128 if (is_curl == PETSC_FALSE) {
1129 int nb_dofs_to_fix = 0;
1130 int index_to_fix = 0;
1131 if (!vertex_to_fix.empty()) {
1132 boost::shared_ptr<NumeredDofEntity> dof_ptr;
1133 CHKERR problem_ptr->getDofByNameEntAndEntDofIdx(
1134 m_field.get_field_bit_number("RHO"), vertex_to_fix[0], 0, ROW,
1135 dof_ptr);
1136 if (dof_ptr) {
1137 if (dof_ptr->getPart() == m_field.get_comm_rank()) {
1138 nb_dofs_to_fix = 1;
1139 index_to_fix = dof_ptr->getPetscGlobalDofIdx();
1140 cerr << *dof_ptr << endl;
1141 }
1142 }
1143 }
1144 CHKERR MatZeroRowsColumns(D, nb_dofs_to_fix, &index_to_fix,
1145 eps_rho / eps_u, PETSC_NULL, PETSC_NULL);
1146 } else {
1147 std::vector<int> dofs_to_fix;
1148 for (auto p_eit = edges_to_fix.pair_begin();
1149 p_eit != edges_to_fix.pair_end(); ++p_eit) {
1150 auto bit_number = m_field.get_field_bit_number("RHO");
1151 auto row_dofs = problem_ptr->numeredRowDofsPtr;
1152 auto lo = row_dofs->lower_bound(
1153 FieldEntity::getLoLocalEntityBitNumber(bit_number, p_eit->first));
1154 auto hi =
1155 row_dofs->lower_bound(FieldEntity::getHiLocalEntityBitNumber(
1156 bit_number, p_eit->second));
1157 for (; lo != hi; ++lo)
1158 dofs_to_fix.push_back((*lo)->getPetscGlobalDofIdx());
1159 }
1160 CHKERR MatZeroRowsColumns(D, dofs_to_fix.size(), &*dofs_to_fix.begin(),
1161 eps_rho / eps_u, PETSC_NULL, PETSC_NULL);
1162 }
1163
1164 // Matrix View
1165 if (debug) {
1166 cerr << "D" << endl;
1167 MatView(D, PETSC_VIEWER_DRAW_WORLD);
1168 std::string wait;
1169 std::cin >> wait;
1170 }
1171
1172 boost::shared_ptr<Problem::SubProblemData> sub_data =
1173 problem_ptr->getSubData();
1174 CHKERR sub_data->getRowIs(&nested_is_rows[1]);
1175 CHKERR sub_data->getColIs(&nested_is_cols[1]);
1176 nested_matrices(1, 1) = D;
1177
1178 CHKERR DMDestroy(&dm_sub_force_penalty);
1179 }
1180
1181 {
1182 DM dm_sub_force;
1183
1184 // Craete dm_control instance
1185 CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_force);
1186 CHKERR DMSetType(dm_sub_force, dm_name);
1187 // set dm_sub_force data structure which created mofem data structures
1188 CHKERR DMMoFEMCreateSubDM(dm_sub_force, dm_control, "FORCES_PROB");
1189 CHKERR DMMoFEMSetSquareProblem(dm_sub_force, PETSC_FALSE);
1190 CHKERR DMMoFEMAddSubFieldRow(dm_sub_force, "RHO");
1191 CHKERR DMMoFEMAddSubFieldCol(dm_sub_force, "U");
1192 CHKERR DMMoFEMAddSubFieldCol(dm_sub_force, "UPSILON");
1193 // add elements to dm_sub_force
1194 CHKERR DMMoFEMAddElement(dm_sub_force, "B");
1195 CHKERR DMSetUp(dm_sub_force);
1196
1197 Mat UB, UPSILONB;
1198 CHKERR DMCreateMatrix(dm_sub_force, &UB);
1199 CHKERR MatZeroEntries(UB);
1200 // note be will be transposed in place latter on
1201 CHKERR DMCreateMatrix(dm_sub_force, &UPSILONB);
1202 CHKERR MatZeroEntries(UPSILONB);
1203 {
1204 CellEngineering::FaceElement face_b_matrices(m_field);
1205 face_b_matrices.getOpPtrVector().push_back(
1206 new OpCalculateInvJacForFace(inv_jac));
1207 if (is_curl) {
1208 face_b_matrices.getOpPtrVector().push_back(
1209 new OpSetInvJacH1ForFace(inv_jac));
1210 face_b_matrices.getOpPtrVector().push_back(
1211 new OpSetInvJacHcurlFace(inv_jac));
1212 face_b_matrices.getOpPtrVector().push_back(new OpCellCurlB(UB, "U"));
1213 face_b_matrices.getOpPtrVector().push_back(
1214 new OpCellCurlB(UPSILONB, "UPSILON"));
1215 } else {
1216 face_b_matrices.getOpPtrVector().push_back(
1217 new OpSetInvJacH1ForFace(inv_jac));
1218 face_b_matrices.getOpPtrVector().push_back(
1219 new OpCellPotentialB(UB, "U"));
1220 face_b_matrices.getOpPtrVector().push_back(
1221 new OpCellPotentialB(UPSILONB, "UPSILON"));
1222 }
1223 CHKERR DMoFEMLoopFiniteElements(dm_sub_force, "B", &face_b_matrices);
1224 }
1225 CHKERR MatAssemblyBegin(UB, MAT_FINAL_ASSEMBLY);
1226 CHKERR MatAssemblyBegin(UPSILONB, MAT_FINAL_ASSEMBLY);
1227 CHKERR MatAssemblyEnd(UB, MAT_FINAL_ASSEMBLY);
1228 CHKERR MatAssemblyEnd(UPSILONB, MAT_FINAL_ASSEMBLY);
1229
1230 const Problem *problem_ptr;
1231 CHKERR m_field.get_problem("FORCES_PROB", &problem_ptr);
1232
1233 // Zero rows, force field is given by gradients of potential field, so one
1234 // of the values has to be fixed like for rigid body motion.
1235 if (is_curl == PETSC_FALSE) {
1236 int nb_dofs_to_fix = 0;
1237 int index_to_fix = 0;
1238 if (!vertex_to_fix.empty()) {
1239 boost::shared_ptr<NumeredDofEntity> dof_ptr;
1240 CHKERR problem_ptr->getDofByNameEntAndEntDofIdx(
1241 m_field.get_field_bit_number("RHO"), vertex_to_fix[0], 0, ROW,
1242 dof_ptr);
1243 if (dof_ptr) {
1244 if (dof_ptr->getPart() == m_field.get_comm_rank()) {
1245 nb_dofs_to_fix = 1;
1246 index_to_fix = dof_ptr->getPetscGlobalDofIdx();
1247 cerr << *dof_ptr << endl;
1248 }
1249 }
1250 }
1251 CHKERR MatZeroRows(UB, nb_dofs_to_fix, &index_to_fix, 0, PETSC_NULL,
1252 PETSC_NULL);
1253 CHKERR MatZeroRows(UPSILONB, nb_dofs_to_fix, &index_to_fix, 0,
1254 PETSC_NULL, PETSC_NULL);
1255 } else {
1256 std::vector<int> dofs_to_fix;
1257 for (auto p_eit = edges_to_fix.pair_begin();
1258 p_eit != edges_to_fix.pair_end(); ++p_eit) {
1259 auto bit_number = m_field.get_field_bit_number("RHO");
1260 auto row_dofs = problem_ptr->numeredRowDofsPtr;
1261 auto lo = row_dofs->lower_bound(
1262 FieldEntity::getLoLocalEntityBitNumber(bit_number, p_eit->first));
1263 auto hi =
1264 row_dofs->lower_bound(FieldEntity::getHiLocalEntityBitNumber(
1265 bit_number, p_eit->second));
1266 for (; lo != hi; ++lo)
1267 dofs_to_fix.push_back((*lo)->getPetscGlobalDofIdx());
1268 }
1269 CHKERR MatZeroRows(UB, dofs_to_fix.size(), &*dofs_to_fix.begin(), 0,
1270 PETSC_NULL, PETSC_NULL);
1271 CHKERR MatZeroRows(UPSILONB, dofs_to_fix.size(), &*dofs_to_fix.begin(),
1272 0, PETSC_NULL, PETSC_NULL);
1273 }
1274
1275 Mat UBT;
1276 CHKERR MatTranspose(UB, MAT_INITIAL_MATRIX, &UBT);
1277 CHKERR MatDestroy(&UB);
1278
1279 // Matrix View
1280 if (debug) {
1281 cerr << "UBT" << endl;
1282 MatView(UBT, PETSC_VIEWER_DRAW_WORLD);
1283 std::string wait;
1284 std::cin >> wait;
1285 }
1286
1287 boost::shared_ptr<Problem::SubProblemData> sub_data =
1288 problem_ptr->getSubData();
1289 // CHKERR sub_data->getColIs(&nested_is_rows[0]);
1290 // CHKERR sub_data->getRowIs(&nested_is_cols[1]);
1291 nested_matrices(0, 1) = UBT;
1292
1293 if (debug) {
1294 cerr << "UPSILONB" << endl;
1295 MatView(UPSILONB, PETSC_VIEWER_DRAW_WORLD);
1296 std::string wait;
1297 std::cin >> wait;
1298 }
1299
1300 // CHKERR sub_data->getRowIs(&nested_is_rows[1]);
1301 // CHKERR sub_data->getColIs(&nested_is_cols[0]);
1302 nested_matrices(1, 0) = UPSILONB;
1303
1304 CHKERR DMDestroy(&dm_sub_force);
1305 }
1306
1307 Mat SubA;
1308 CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &sub_nested_is_rows[0], 2,
1309 &sub_nested_is_cols[0], &sub_nested_matrices(0, 0),
1310 &SubA);
1311 nested_matrices(0, 0) = SubA;
1312
1313 CHKERR MatAssemblyBegin(SubA, MAT_FINAL_ASSEMBLY);
1314 CHKERR MatAssemblyEnd(SubA, MAT_FINAL_ASSEMBLY);
1315
1316 if (debug) {
1317 cerr << "Nested SubA" << endl;
1318 MatView(SubA, PETSC_VIEWER_STDOUT_WORLD);
1319 }
1320
1321 Mat A;
1322 CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &nested_is_rows[0], 2,
1323 &nested_is_cols[0], &nested_matrices(0, 0), &A);
1324
1325 CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1326 CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1327
1328 if (debug) {
1329 cerr << "Nested A" << endl;
1330 MatView(A, PETSC_VIEWER_STDOUT_WORLD);
1331 }
1332
1333 Vec D, F;
1334 CHKERR DMCreateGlobalVector(dm_control, &D);
1335 CHKERR DMCreateGlobalVector(dm_control, &F);
1336
1337 // Asseble the right hand side vector
1338 {
1339 CellEngineering::FaceElement face_element(m_field);
1340 face_element.getOpPtrVector().push_back(new OpGetDispX(common_data));
1341 face_element.getOpPtrVector().push_back(new OpGetDispY(common_data));
1342 face_element.getOpPtrVector().push_back(
1343 new OpCell_g(F, eps_u, common_data));
1344 CHKERR DMoFEMLoopFiniteElements(dm_control, "DISPLACEMENTS_PENALTY",
1345 &face_element);
1346 CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
1347 CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
1348 CHKERR VecAssemblyBegin(F);
1349 CHKERR VecAssemblyEnd(F);
1350 }
1351
1352 KSP solver;
1353 // KSP solver;
1354 {
1355 CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
1356 CHKERR KSPSetDM(solver, dm_control);
1357 CHKERR KSPSetFromOptions(solver);
1358 CHKERR KSPSetOperators(solver, A, A);
1359 CHKERR KSPSetDMActive(solver, PETSC_FALSE);
1360 CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
1361 CHKERR KSPSetInitialGuessNonzero(solver, PETSC_FALSE);
1362 PC pc;
1363 CHKERR KSPGetPC(solver, &pc);
1364 CHKERR PCSetType(pc, PCFIELDSPLIT);
1365 PetscBool is_pcfs = PETSC_FALSE;
1366 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1367 if (is_pcfs) {
1368 CHKERR PCSetOperators(pc, A, A);
1369 CHKERR PCFieldSplitSetIS(pc, NULL, nested_is_rows[0]);
1370 CHKERR PCFieldSplitSetIS(pc, NULL, nested_is_rows[1]);
1371 CHKERR PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
1372 CHKERR PCSetUp(pc);
1373 KSP *sub_ksp;
1374 PetscInt n;
1375 CHKERR PCFieldSplitGetSubKSP(pc, &n, &sub_ksp);
1376 {
1377 PC sub_pc_0;
1378 CHKERR KSPGetPC(sub_ksp[0], &sub_pc_0);
1379 CHKERR PCSetOperators(sub_pc_0, SubA, SubA);
1380 CHKERR PCSetType(sub_pc_0, PCFIELDSPLIT);
1381 CHKERR PCFieldSplitSetIS(sub_pc_0, NULL, sub_nested_is_rows[0]);
1382 CHKERR PCFieldSplitSetIS(sub_pc_0, NULL, sub_nested_is_rows[1]);
1383 CHKERR PCFieldSplitSetType(sub_pc_0, PC_COMPOSITE_MULTIPLICATIVE);
1384 // CHKERR
1385 // PCFieldSplitSetSchurFactType(sub_pc_0,PC_FIELDSPLIT_SCHUR_FACT_LOWER);
1386 // CHKERR PCFieldSplitSetType(sub_pc_0,PC_COMPOSITE_SCHUR);
1387 CHKERR PCSetUp(sub_pc_0);
1388 }
1389 } else {
1390 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1391 "Only works with pre-conditioner PCFIELDSPLIT");
1392 }
1393 CHKERR KSPSetUp(solver);
1394 }
1395
1396 // Solve system of equations
1397 CHKERR KSPSolve(solver, F, D);
1398
1399 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1400 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1401 CHKERR DMoFEMMeshToLocalVector(dm_control, D, INSERT_VALUES,
1402 SCATTER_REVERSE);
1403
1404 if (debug) {
1405 CHKERR VecView(D, PETSC_VIEWER_DRAW_WORLD);
1406 std::string wait;
1407 std::cin >> wait;
1408 }
1409
1410 // Clean sub matrices and sub indices
1411 for (int i = 0; i != 2; i++) {
1412 if (sub_nested_is_rows[i]) {
1413 CHKERR ISDestroy(&sub_nested_is_rows[i]);
1414 }
1415 if (sub_nested_is_cols[i]) {
1416 CHKERR ISDestroy(&sub_nested_is_cols[i]);
1417 }
1418 for (int j = 0; j != 2; j++) {
1419 if (sub_nested_matrices(i, j)) {
1420 CHKERR MatDestroy(&sub_nested_matrices(i, j));
1421 }
1422 }
1423 }
1424 for (int i = 0; i != 2; i++) {
1425 if (nested_is_rows[i]) {
1426 CHKERR ISDestroy(&nested_is_rows[i]);
1427 }
1428 if (nested_is_cols[i]) {
1429 CHKERR ISDestroy(&nested_is_cols[i]);
1430 }
1431 for (int j = 0; j != 2; j++) {
1432 if (nested_matrices(i, j)) {
1433 CHKERR MatDestroy(&nested_matrices(i, j));
1434 }
1435 }
1436 }
1437
1438 CHKERR MatDestroy(&SubA);
1439 CHKERR MatDestroy(&A);
1440 CHKERR VecDestroy(&D);
1441 CHKERR VecDestroy(&F);
1442
1443 CHKERR DMDestroy(&dm_sub_volume_control);
1444
1445 PostProcVolumeOnRefinedMesh post_proc(m_field);
1446 {
1447 CHKERR post_proc.generateReferenceElementMesh();
1448 CHKERR post_proc.addFieldValuesPostProc("U");
1449 CHKERR post_proc.addFieldValuesGradientPostProc("U");
1450 // add postpocessing for sresses
1451 post_proc.getOpPtrVector().push_back(new PostProcHookStress(
1452 m_field, post_proc.postProcMesh, post_proc.mapGaussPts, "U",
1453 post_proc.commonData, &elastic.setOfBlocks));
1454 CHKERR DMoFEMLoopFiniteElements(dm_control, "ELASTIC", &post_proc);
1455 CHKERR post_proc.writeFile("out.h5m");
1456 elastic.getLoopFeEnergy().snes_ctx = SnesMethod::CTX_SNESNONE;
1457 elastic.getLoopFeEnergy().eNergy = 0;
1458 CHKERR DMoFEMLoopFiniteElements(dm_control, "ELASTIC",
1459 &elastic.getLoopFeEnergy());
1460 PetscPrintf(PETSC_COMM_WORLD, "Elastic energy %6.4e\n",
1461 elastic.getLoopFeEnergy().eNergy);
1462 }
1463
1464 {
1465 PostProcFaceOnRefinedMesh post_proc_face(m_field);
1466 CHKERR post_proc_face.generateReferenceElementMesh();
1467 post_proc_face.getOpPtrVector().push_back(
1468 new OpCalculateInvJacForFace(inv_jac));
1469 if (is_curl) {
1470 post_proc_face.getOpPtrVector().push_back(
1471 new OpSetInvJacHcurlFace(inv_jac));
1472 post_proc_face.getOpPtrVector().push_back(
1473 new OpVirtualCurlRho("RHO", common_data));
1474 } else {
1475 post_proc_face.getOpPtrVector().push_back(
1476 new OpSetInvJacH1ForFace(inv_jac));
1477 post_proc_face.getOpPtrVector().push_back(
1478 new OpVirtualPotentialRho("RHO", common_data));
1479 }
1480 post_proc_face.getOpPtrVector().push_back(
1481 new PostProcTraction(m_field, post_proc_face.postProcMesh,
1482 post_proc_face.mapGaussPts, common_data));
1483 CHKERR DMoFEMLoopFiniteElements(dm_control, "D", &post_proc_face);
1484 CHKERR post_proc_face.writeFile("out_tractions.h5m");
1485 }
1486
1487 CHKERR DMDestroy(&dm_control);
1488
1489 }
1490 CATCH_ERRORS;
1491
1492 MoFEM::Core::Finalize();
1493 return 0;
1494}
help
static char help[]
Definition cell_forces.cpp:409
ROW
@ ROW
Definition definitions.h:136
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition definitions.h:385
MF_ZERO
@ MF_ZERO
Definition definitions.h:111
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60
H1
@ H1
continuous field
Definition definitions.h:85
HCURL
@ HCURL
field with continuous tangents
Definition definitions.h:86
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition definitions.h:228
SIDESET
@ SIDESET
Definition definitions.h:160
BLOCKSET
@ BLOCKSET
Definition definitions.h:161
MOFEM_STD_EXCEPTION_THROW
@ MOFEM_STD_EXCEPTION_THROW
Definition definitions.h:39
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition definitions.h:36
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:548
order
constexpr int order
Definition dg_projection.cpp:18
F
@ F
Definition free_surface.cpp:394
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition DMMoFEM.cpp:1123
MoFEM::DMMoFEMCreateSubDM
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition DMMoFEM.cpp:215
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:497
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition DMMoFEM.cpp:456
MoFEM::DMMoFEMCreateMoFEM
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition DMMoFEM.cpp:114
MoFEM::DMoFEMPostProcessFiniteElements
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:556
MoFEM::DMoFEMMeshToLocalVector
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:523
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition DMMoFEM.cpp:238
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:586
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition DMMoFEM.cpp:280
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:546
MoFEM::CoreInterface::get_problem
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MoFEM::CoreInterface::modify_finite_element_add_field_col
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
MoFEM::CoreInterface::add_ents_to_finite_element_by_type
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
MoFEM::CoreInterface::modify_finite_element_add_field_row
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
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEM::DeprecatedCoreInterface::synchronise_field_entities
DEPRECATED MoFEMErrorCode synchronise_field_entities(const std::string &name, int verb=DEFAULT_VERBOSITY)
Definition DeprecatedCoreInterface.cpp:47
MoFEM::CoreInterface::set_field_order
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.
MoFEM::CoreInterface::add_ents_to_field_by_type
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.
MoFEM::MeshsetsManager::checkMeshset
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
Definition MeshsetsManager.cpp:360
MoFEM::MeshsetsManager::getEntitiesByDimension
MoFEMErrorCode getEntitiesByDimension(const int ms_id, const unsigned int cubit_bc_type, const int dimension, Range &entities, const bool recursive=true) const
get entities from CUBIT/meshset of a particular entity dimension
Definition MeshsetsManager.cpp:669
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Definition MeshsetsManager.hpp:71
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
D
double D
Definition initial_diffusion.cpp:44
n
const double n
refractive index of diffusive medium
Definition initial_diffusion.cpp:38
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:19
mesh_file_name
char mesh_file_name[255]
Definition mesh_smoothing.cpp:22
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition Exceptions.hpp:76
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
MoFEM::debug
static const bool debug
Definition ConstrainMatrixCtx.cpp:9
MoFEM::OpSetInvJacHcurlFace
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
Definition UserDataOperators.hpp:2978
A
constexpr AssemblyType A
Definition operators_tests.cpp:30
CellEngineering::OpCell_g
Calculate and assemble g vector.
Definition CellForces.hpp:268
CellEngineering::OpCellCurlB
Calculate and assemble Z matrix.
Definition CellForces.hpp:393
CellEngineering::OpCellCurlD
Calculate and assemble Z matrix.
Definition CellForces.hpp:316
CellEngineering::OpCellPotentialB
Calculate and assemble B matrix.
Definition CellForces.hpp:141
CellEngineering::OpCellPotentialD
Calculate and assemble D matrix.
Definition CellForces.hpp:74
CellEngineering::OpCellS
Calculate and assemble S matrix.
Definition CellForces.hpp:202
CellEngineering::OpVirtualCurlRho
Post-process tractions.
Definition CellForces.hpp:529
CellEngineering::PostProcTraction
Shave results on mesh tags for post-processing.
Definition CellForces.hpp:601
DirichletDisplacementBc
Set Dirichlet boundary conditions on displacements.
Definition DirichletBC.hpp:57
Hooke
Hook equation.
Definition Hooke.hpp:21
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
MoFEM::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEM::CoreInterface::add_field
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.
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition Core.hpp:82
MoFEM::CoreTmp< 0 >::Initialize
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
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::DeprecatedCoreInterface::synchronise_entities
DEPRECATED MoFEMErrorCode synchronise_entities(Range &ent, int verb=DEFAULT_VERBOSITY)
Definition DeprecatedCoreInterface.cpp:41
MoFEM::DeprecatedCoreInterface::seed_ref_level
DEPRECATED MoFEMErrorCode seed_ref_level(const Range &ents, const BitRefLevel &bit, const bool only_tets=true, int verb=-1)
seed entities in the range and their adjacencies in a particular BitRefLevel
Definition DeprecatedCoreInterface.cpp:24
MoFEM::FieldEntity::getLoLocalEntityBitNumber
static UId getLoLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
Definition FieldEntsMultiIndices.hpp:247
MoFEM::FieldEntity::getHiLocalEntityBitNumber
static UId getHiLocalEntityBitNumber(const char bit_number, const EntityHandle ent)
Definition FieldEntsMultiIndices.hpp:258
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition ISManager.hpp:23
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition MeshsetsManager.hpp:104
MoFEM::MeshsetsManager::printForceSet
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
Definition MeshsetsManager.cpp:304
MoFEM::MeshsetsManager::printMaterialsSet
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
Definition MeshsetsManager.cpp:325
MoFEM::MeshsetsManager::printDisplacementSet
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Definition MeshsetsManager.cpp:290
MoFEM::OpCalculateInvJacForFatPrism
Calculate inverse of jacobian for face element.
Definition UserDataOperators.hpp:3146
MoFEM::OpSetInvJacH1ForFatPrism
Transform local reference derivatives of shape functions to global derivatives.
Definition UserDataOperators.hpp:3173
MoFEM::Problem
keeps basic data about problem
Definition ProblemsMultiIndices.hpp:54
MoFEM::Problem::getDofByNameEntAndEntDofIdx
MoFEMErrorCode getDofByNameEntAndEntDofIdx(const int field_bit_number, const EntityHandle ent, const int ent_dof_idx, const RowColData row_or_col, boost::shared_ptr< NumeredDofEntity > &dof_ptr) const
get DOFs from problem
Definition ProblemsMultiIndices.cpp:132
MoFEM::Problem::getSubData
boost::shared_ptr< SubProblemData > & getSubData() const
Get main problem of sub-problem is.
Definition ProblemsMultiIndices.hpp:119
MoFEM::Problem::numeredRowDofsPtr
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem
Definition ProblemsMultiIndices.hpp:73
MoFEM::UnknownInterface::query_interface
virtual MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const =0
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx
Definition NonLinearElasticElement.hpp:557
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition HookeElement.hpp:27
PostProcFaceOnRefinedMesh
Postprocess on face.
Definition PostProcOnRefMesh.hpp:1033
PostProcHookStress
Operator post-procesing stresses for Hook isotropic material.
Definition PostProcHookStresses.hpp:41
PostProcVolumeOnRefinedMesh
Post processing.
Definition PostProcOnRefMesh.hpp:957

Variable Documentation

◆ debug

int debug = 0
static

Definition at line 430 of file cell_forces.cpp.

◆ help

char help[]
static
Initial value:
= "-my_block_config set block data\n"
"\n"

Definition at line 409 of file cell_forces.cpp.

Generated by Doxygen 1.12.0 and hosted at University of Glasgow