v0.10.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));
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);
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);
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);
751  elastic.setOfBlocks[2].tEts, MBPRISM, "ELASTIC_PRISM");
752  MatrixDouble inv_jac;
753  // right hand side operators
754  fat_prism_rhs.getOpPtrVector().push_back(
756  fat_prism_rhs.getOpPtrVector().push_back(
757  new MoFEM::OpSetInvJacH1ForFatPrism(inv_jac));
758  fat_prism_rhs.getOpPtrVector().push_back(
760  "U", elastic.commonData));
761  fat_prism_rhs.getOpPtrVector().push_back(
763  "U", elastic.setOfBlocks[2], elastic.commonData, 2, false, false,
764  true));
765  fat_prism_rhs.getOpPtrVector().push_back(
767  "U", elastic.setOfBlocks[2], elastic.commonData));
768  // Left hand side operators
769  fat_prism_lhs.getOpPtrVector().push_back(
771  fat_prism_lhs.getOpPtrVector().push_back(
772  new MoFEM::OpSetInvJacH1ForFatPrism(inv_jac));
773  fat_prism_lhs.getOpPtrVector().push_back(
775  "U", elastic.commonData));
776  fat_prism_lhs.getOpPtrVector().push_back(
778  "U", elastic.setOfBlocks[2], elastic.commonData, 2, true, false,
779  true));
780  fat_prism_lhs.getOpPtrVector().push_back(
782  "U", "U", elastic.setOfBlocks[2], elastic.commonData));
783  }
784 
785  // Update material parameters
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");
843  CHKERR m_field.modify_finite_element_add_field_col("B", "UPSILON");
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 
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 
1493  return 0;
1494 }
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
Deprecated interface functions.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
Calculate and assemble Z matrix.
Definition: CellForces.hpp:315
Calculate and assemble Z matrix.
Definition: CellForces.hpp:392
CoreTmp< 0 > Core
Definition: Core.hpp:1129
Hook equation.
Definition: Hooke.hpp:31
boost::shared_ptr< SubProblemData > & getSubData() const
Get main problem of sub-problem is.
const double D
diffusivity
Section manager is used to create indexes and sectionsFIXME: ISManager is not properly testsed by ato...
Definition: ISManager.hpp:36
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:445
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problemIt if true is assumed that matrix has the same indexing on rows and columns....
Definition: DMMMoFEM.cpp:378
auto post_proc
static char help[]
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
Definition: DMMMoFEM.cpp:221
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 dimensionNodeset can contain nodes,...
PetscErrorCode DMMoFEMAddElement(DM dm, const char fe_name[])
add element to dm
Definition: DMMMoFEM.cpp:425
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
Calculate and assemble S matrix.
Definition: CellForces.hpp:201
DEPRECATED MoFEMErrorCode query_interface(IFace *&ptr) const
Calculate and assemble g vector.
Definition: CellForces.hpp:267
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:60
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.
Interface for managing meshsets containing materials and boundary conditions.
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.
boost::shared_ptr< NumeredDofEntity_multiIndex > numeredRowDofsPtr
store DOFs on rows for this problem
keeps basic data about problemThis is low level structure with information about problem,...
static int debug
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
char mesh_file_name[255]
static Index< 'n', 3 > n
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition: DMMMoFEM.cpp:177
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
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
Definition: DMMMoFEM.cpp:198
Operator post-procesing stresses for Hook isotropic material.
DEPRECATED MoFEMErrorCode synchronise_entities(Range &ent, int verb=DEFAULT_VERBOSITY)
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:87
virtual int get_comm_rank() const =0
constexpr int order
Calculate inverse of jacobian for face element.
static Index< 'i', 3 > i
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
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.
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:69
static Index< 'j', 3 > j
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Calculate inverse of jacobian for face element.
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
Post processing.
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMMoFEM.cpp:468
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
field with continuous tangents
Definition: definitions.h:178
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMMoFEM.cpp:478
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: DMMMoFEM.cpp:105
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
#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.
#define CHKERR
Inline error check.
Definition: definitions.h:604
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
Transform local reference derivatives of shape functions to global derivatives.
Transform local reference derivatives of shape functions to global derivatives.
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:292
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
Postprocess on face.
brief Transform local reference derivatives of shape function to global derivatives for face
virtual const Problem * get_problem(const std::string &problem_name) const =0
Get the problem object.
Post-process tractions.
Definition: CellForces.hpp:528
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
Shave results on mesh tags for post-processing.
Definition: CellForces.hpp:600
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
Calculate and assemble D matrix.
Definition: CellForces.hpp:73
Core (interface) class.
Definition: Core.hpp:77
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:441
DEPRECATED MoFEMErrorCode synchronise_field_entities(const std::string &name, int verb=DEFAULT_VERBOSITY)
continuous field
Definition: definitions.h:177
Calculate and assemble B matrix.
Definition: CellForces.hpp:140
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:982
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
structure grouping operators and data used for calculation of nonlinear elastic elementIn order to as...

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.