v0.14.0
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

 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[] 
)
Examples
cell_forces.cpp.

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 }

Variable Documentation

◆ debug

int debug = 0
static
Examples
cell_forces.cpp.

Definition at line 430 of file cell_forces.cpp.

◆ help

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

Definition at line 409 of file cell_forces.cpp.

PostProcHookStress
Operator post-procesing stresses for Hook isotropic material.
Definition: PostProcHookStresses.hpp:40
SIDESET
@ SIDESET
Definition: definitions.h:147
debug
static int debug
Definition: cell_forces.cpp:430
CellEngineering::PostProcTraction
Shave results on mesh tags for post-processing.
Definition: CellForces.hpp:600
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
MoFEM::OpSetInvJacHcurlFace
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
Definition: UserDataOperators.hpp:2982
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::DMMoFEMAddSubFieldCol
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:284
MOFEM_STD_EXCEPTION_THROW
@ MOFEM_STD_EXCEPTION_THROW
Definition: definitions.h:39
NonlinearElasticElement::OpJacobianPiolaKirchhoffStress
Operator performs automatic differentiation.
Definition: NonLinearElasticElement.hpp:370
MoFEM::UnknownInterface::query_interface
virtual MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const =0
CellEngineering::FaceElement
Definition: CellForces.hpp:48
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
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
MoFEM::DMMoFEMSetSquareProblem
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:460
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
CellEngineering::OpCellCurlB
Calculate and assemble Z matrix.
Definition: CellForces.hpp:392
MoFEM::CoreInterface::get_field_bit_number
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::CoreInterface::get_problem
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
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:527
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
DirichletDisplacementBc
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:57
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::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:501
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::OpCalculateInvJacForFatPrism
Calculate inverse of jacobian for face element.
Definition: UserDataOperators.hpp:3150
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
ROW
@ ROW
Definition: definitions.h:123
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
help
static char help[]
Definition: cell_forces.cpp:409
CellEngineering::OpVirtualCurlRho
Post-process tractions.
Definition: CellForces.hpp:528
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
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::MeshsetsManager::printForceSet
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
Definition: MeshsetsManager.cpp:301
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::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::OpSetInvJacH1ForFatPrism
Transform local reference derivatives of shape functions to global derivatives.
Definition: UserDataOperators.hpp:3177
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
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:666
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::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:219
CellEngineering::OpCellCurlD
Calculate and assemble Z matrix.
Definition: CellForces.hpp:315
CellEngineering::OpCellPotentialD
Calculate and assemble D matrix.
Definition: CellForces.hpp:73
MoFEM::Problem::getSubData
boost::shared_ptr< SubProblemData > & getSubData() const
Get main problem of sub-problem is.
Definition: ProblemsMultiIndices.hpp:119
NonlinearElasticElement::OpLhsPiolaKirchhoff_dx
Definition: NonLinearElasticElement.hpp:556
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
MoFEM::DeprecatedCoreInterface::synchronise_entities
DEPRECATED MoFEMErrorCode synchronise_entities(Range &ent, int verb=DEFAULT_VERBOSITY)
Definition: DeprecatedCoreInterface.cpp:41
CellEngineering::FatPrism
Definition: CellForces.hpp:21
MoFEM::DMoFEMPreProcessFiniteElements
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:550
CellEngineering::OpCellS
Calculate and assemble S matrix.
Definition: CellForces.hpp:201
MoFEM::DMoFEMPostProcessFiniteElements
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMoFEM.cpp:560
CellEngineering::OpGetDispX
Definition: CellForces.hpp:57
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:118
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::MeshsetsManager::printMaterialsSet
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
Definition: MeshsetsManager.cpp:322
PostProcFaceOnRefinedMesh
Postprocess on face.
Definition: PostProcOnRefMesh.hpp:1032
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
convert.n
n
Definition: convert.py:82
MoFEM::Problem::numeredRowDofsPtr
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem
Definition: ProblemsMultiIndices.hpp:73
MoFEM::OpSetInvJacH1ForFace
Definition: UserDataOperators.hpp:2922
Range
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
MF_ZERO
@ MF_ZERO
Definition: definitions.h:98
NonlinearElasticElement::OpGetCommonDataAtGaussPts
Definition: NonLinearElasticElement.hpp:362
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
NonlinearElasticElement::OpRhsPiolaKirchhoff
Definition: NonLinearElasticElement.hpp:520
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
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:955
MoFEM::MeshsetsManager::printDisplacementSet
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Definition: MeshsetsManager.cpp:287
_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
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1094
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
j
FTensor::Index< 'j', 3 > j
Definition: matrix_function.cpp:19
CellEngineering::OpCellPotentialB
Calculate and assemble B matrix.
Definition: CellForces.hpp:140
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
CellEngineering::CommonData
Definition: CellForces.hpp:31
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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_filed)=0
set finite element field data
CellEngineering::OpGetDispY
Definition: CellForces.hpp:65
MoFEM::DeprecatedCoreInterface::synchronise_field_entities
DEPRECATED MoFEMErrorCode synchronise_field_entities(const std::string &name, int verb=DEFAULT_VERBOSITY)
Definition: DeprecatedCoreInterface.cpp:47
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
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::CoreInterface::build_adjacencies
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
MoFEM::Problem
keeps basic data about problem
Definition: ProblemsMultiIndices.hpp:54
MoFEM::DMMoFEMAddSubFieldRow
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:242
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::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::MeshsetsManager::checkMeshset
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
Definition: MeshsetsManager.cpp:357
CellEngineering::OpCell_g
Calculate and assemble g vector.
Definition: CellForces.hpp:267
CellEngineering::OpVirtualPotentialRho
Post-process tractions.
Definition: CellForces.hpp:460
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:590
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
F
@ F
Definition: free_surface.cpp:394
Hooke
Hook equation.
Definition: Hooke.hpp:19
MoFEM::DMMoFEMSetIsPartitioned
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMoFEM.cpp:1127
Generated by Doxygen 1.8.17 and hosted at University of Glasgow