v0.9.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  "RHO", vertex_to_fix[0], 0, ROW, dof_ptr);
1135  if (dof_ptr) {
1136  if (dof_ptr->getPart() == m_field.get_comm_rank()) {
1137  nb_dofs_to_fix = 1;
1138  index_to_fix = dof_ptr->getPetscGlobalDofIdx();
1139  cerr << *dof_ptr << endl;
1140  }
1141  }
1142  }
1143  CHKERR MatZeroRowsColumns(D, nb_dofs_to_fix, &index_to_fix,
1144  eps_rho / eps_u, PETSC_NULL, PETSC_NULL);
1145  } else {
1146  std::vector<int> dofs_to_fix;
1147  for (Range::iterator eit = edges_to_fix.begin();
1148  eit != edges_to_fix.end(); eit++) {
1150  problem_ptr, "RHO", *eit, m_field.get_comm_rank(), dof_it)) {
1151  dofs_to_fix.push_back(dof_it->get()->getPetscGlobalDofIdx());
1152  }
1153  }
1154  CHKERR MatZeroRowsColumns(D, dofs_to_fix.size(), &*dofs_to_fix.begin(),
1155  eps_rho / eps_u, PETSC_NULL, PETSC_NULL);
1156  }
1157 
1158  // Matrix View
1159  if (debug) {
1160  cerr << "D" << endl;
1161  MatView(D, PETSC_VIEWER_DRAW_WORLD);
1162  std::string wait;
1163  std::cin >> wait;
1164  }
1165 
1166  boost::shared_ptr<Problem::SubProblemData> sub_data =
1167  problem_ptr->getSubData();
1168  CHKERR sub_data->getRowIs(&nested_is_rows[1]);
1169  CHKERR sub_data->getColIs(&nested_is_cols[1]);
1170  nested_matrices(1, 1) = D;
1171 
1172  CHKERR DMDestroy(&dm_sub_force_penalty);
1173  }
1174 
1175  {
1176  DM dm_sub_force;
1177 
1178  // Craete dm_control instance
1179  CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_force);
1180  CHKERR DMSetType(dm_sub_force, dm_name);
1181  // set dm_sub_force data structure which created mofem data structures
1182  CHKERR DMMoFEMCreateSubDM(dm_sub_force, dm_control, "FORCES_PROB");
1183  CHKERR DMMoFEMSetSquareProblem(dm_sub_force, PETSC_FALSE);
1184  CHKERR DMMoFEMAddSubFieldRow(dm_sub_force, "RHO");
1185  CHKERR DMMoFEMAddSubFieldCol(dm_sub_force, "U");
1186  CHKERR DMMoFEMAddSubFieldCol(dm_sub_force, "UPSILON");
1187  // add elements to dm_sub_force
1188  CHKERR DMMoFEMAddElement(dm_sub_force, "B");
1189  CHKERR DMSetUp(dm_sub_force);
1190 
1191  Mat UB, UPSILONB;
1192  CHKERR DMCreateMatrix(dm_sub_force, &UB);
1193  CHKERR MatZeroEntries(UB);
1194  // note be will be transposed in place latter on
1195  CHKERR DMCreateMatrix(dm_sub_force, &UPSILONB);
1196  CHKERR MatZeroEntries(UPSILONB);
1197  {
1198  CellEngineering::FaceElement face_b_matrices(m_field);
1199  face_b_matrices.getOpPtrVector().push_back(
1200  new OpCalculateInvJacForFace(inv_jac));
1201  if (is_curl) {
1202  face_b_matrices.getOpPtrVector().push_back(
1203  new OpSetInvJacH1ForFace(inv_jac));
1204  face_b_matrices.getOpPtrVector().push_back(
1205  new OpSetInvJacHcurlFace(inv_jac));
1206  face_b_matrices.getOpPtrVector().push_back(new OpCellCurlB(UB, "U"));
1207  face_b_matrices.getOpPtrVector().push_back(
1208  new OpCellCurlB(UPSILONB, "UPSILON"));
1209  } else {
1210  face_b_matrices.getOpPtrVector().push_back(
1211  new OpSetInvJacH1ForFace(inv_jac));
1212  face_b_matrices.getOpPtrVector().push_back(
1213  new OpCellPotentialB(UB, "U"));
1214  face_b_matrices.getOpPtrVector().push_back(
1215  new OpCellPotentialB(UPSILONB, "UPSILON"));
1216  }
1217  CHKERR DMoFEMLoopFiniteElements(dm_sub_force, "B", &face_b_matrices);
1218  }
1219  CHKERR MatAssemblyBegin(UB, MAT_FINAL_ASSEMBLY);
1220  CHKERR MatAssemblyBegin(UPSILONB, MAT_FINAL_ASSEMBLY);
1221  CHKERR MatAssemblyEnd(UB, MAT_FINAL_ASSEMBLY);
1222  CHKERR MatAssemblyEnd(UPSILONB, MAT_FINAL_ASSEMBLY);
1223 
1224  const Problem *problem_ptr;
1225  CHKERR m_field.get_problem("FORCES_PROB", &problem_ptr);
1226 
1227  // Zero rows, force field is given by gradients of potential field, so one
1228  // of the values has to be fixed like for rigid body motion.
1229  if (is_curl == PETSC_FALSE) {
1230  int nb_dofs_to_fix = 0;
1231  int index_to_fix = 0;
1232  if (!vertex_to_fix.empty()) {
1233  boost::shared_ptr<NumeredDofEntity> dof_ptr;
1234  CHKERR problem_ptr->getDofByNameEntAndEntDofIdx(
1235  "RHO", vertex_to_fix[0], 0, ROW, dof_ptr);
1236  if (dof_ptr) {
1237  if (dof_ptr->getPart() == m_field.get_comm_rank()) {
1238  nb_dofs_to_fix = 1;
1239  index_to_fix = dof_ptr->getPetscGlobalDofIdx();
1240  cerr << *dof_ptr << endl;
1241  }
1242  }
1243  }
1244  CHKERR MatZeroRows(UB, nb_dofs_to_fix, &index_to_fix, 0, PETSC_NULL,
1245  PETSC_NULL);
1246  CHKERR MatZeroRows(UPSILONB, nb_dofs_to_fix, &index_to_fix, 0,
1247  PETSC_NULL, PETSC_NULL);
1248  } else {
1249  std::vector<int> dofs_to_fix;
1250  for (Range::iterator eit = edges_to_fix.begin();
1251  eit != edges_to_fix.end(); eit++) {
1253  problem_ptr, "RHO", *eit, m_field.get_comm_rank(), dof_it)) {
1254  dofs_to_fix.push_back(dof_it->get()->getPetscGlobalDofIdx());
1255  }
1256  }
1257  CHKERR MatZeroRows(UB, dofs_to_fix.size(), &*dofs_to_fix.begin(), 0,
1258  PETSC_NULL, PETSC_NULL);
1259  CHKERR MatZeroRows(UPSILONB, dofs_to_fix.size(), &*dofs_to_fix.begin(),
1260  0, PETSC_NULL, PETSC_NULL);
1261  }
1262 
1263  Mat UBT;
1264  CHKERR MatTranspose(UB, MAT_INITIAL_MATRIX, &UBT);
1265  CHKERR MatDestroy(&UB);
1266 
1267  // Matrix View
1268  if (debug) {
1269  cerr << "UBT" << endl;
1270  MatView(UBT, PETSC_VIEWER_DRAW_WORLD);
1271  std::string wait;
1272  std::cin >> wait;
1273  }
1274 
1275  boost::shared_ptr<Problem::SubProblemData> sub_data =
1276  problem_ptr->getSubData();
1277  // CHKERR sub_data->getColIs(&nested_is_rows[0]);
1278  // CHKERR sub_data->getRowIs(&nested_is_cols[1]);
1279  nested_matrices(0, 1) = UBT;
1280 
1281  if (debug) {
1282  cerr << "UPSILONB" << endl;
1283  MatView(UPSILONB, PETSC_VIEWER_DRAW_WORLD);
1284  std::string wait;
1285  std::cin >> wait;
1286  }
1287 
1288  // CHKERR sub_data->getRowIs(&nested_is_rows[1]);
1289  // CHKERR sub_data->getColIs(&nested_is_cols[0]);
1290  nested_matrices(1, 0) = UPSILONB;
1291 
1292  CHKERR DMDestroy(&dm_sub_force);
1293  }
1294 
1295  Mat SubA;
1296  CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &sub_nested_is_rows[0], 2,
1297  &sub_nested_is_cols[0], &sub_nested_matrices(0, 0),
1298  &SubA);
1299  nested_matrices(0, 0) = SubA;
1300 
1301  CHKERR MatAssemblyBegin(SubA, MAT_FINAL_ASSEMBLY);
1302  CHKERR MatAssemblyEnd(SubA, MAT_FINAL_ASSEMBLY);
1303 
1304  if (debug) {
1305  cerr << "Nested SubA" << endl;
1306  MatView(SubA, PETSC_VIEWER_STDOUT_WORLD);
1307  }
1308 
1309  Mat A;
1310  CHKERR MatCreateNest(PETSC_COMM_WORLD, 2, &nested_is_rows[0], 2,
1311  &nested_is_cols[0], &nested_matrices(0, 0), &A);
1312 
1313  CHKERR MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1314  CHKERR MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1315 
1316  if (debug) {
1317  cerr << "Nested A" << endl;
1318  MatView(A, PETSC_VIEWER_STDOUT_WORLD);
1319  }
1320 
1321  Vec D, F;
1322  CHKERR DMCreateGlobalVector(dm_control, &D);
1323  CHKERR DMCreateGlobalVector(dm_control, &F);
1324 
1325  // Asseble the right hand side vector
1326  {
1327  CellEngineering::FaceElement face_element(m_field);
1328  face_element.getOpPtrVector().push_back(new OpGetDispX(common_data));
1329  face_element.getOpPtrVector().push_back(new OpGetDispY(common_data));
1330  face_element.getOpPtrVector().push_back(
1331  new OpCell_g(F, eps_u, common_data));
1332  CHKERR DMoFEMLoopFiniteElements(dm_control, "DISPLACEMENTS_PENALTY",
1333  &face_element);
1334  CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
1335  CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
1336  CHKERR VecAssemblyBegin(F);
1337  CHKERR VecAssemblyEnd(F);
1338  }
1339 
1340  KSP solver;
1341  // KSP solver;
1342  {
1343  CHKERR KSPCreate(PETSC_COMM_WORLD, &solver);
1344  CHKERR KSPSetDM(solver, dm_control);
1345  CHKERR KSPSetFromOptions(solver);
1346  CHKERR KSPSetOperators(solver, A, A);
1347  CHKERR KSPSetDMActive(solver, PETSC_FALSE);
1348  CHKERR KSPSetInitialGuessKnoll(solver, PETSC_FALSE);
1349  CHKERR KSPSetInitialGuessNonzero(solver, PETSC_FALSE);
1350  PC pc;
1351  CHKERR KSPGetPC(solver, &pc);
1352  CHKERR PCSetType(pc, PCFIELDSPLIT);
1353  PetscBool is_pcfs = PETSC_FALSE;
1354  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
1355  if (is_pcfs) {
1356  CHKERR PCSetOperators(pc, A, A);
1357  CHKERR PCFieldSplitSetIS(pc, NULL, nested_is_rows[0]);
1358  CHKERR PCFieldSplitSetIS(pc, NULL, nested_is_rows[1]);
1359  CHKERR PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
1360  CHKERR PCSetUp(pc);
1361  KSP *sub_ksp;
1362  PetscInt n;
1363  CHKERR PCFieldSplitGetSubKSP(pc, &n, &sub_ksp);
1364  {
1365  PC sub_pc_0;
1366  CHKERR KSPGetPC(sub_ksp[0], &sub_pc_0);
1367  CHKERR PCSetOperators(sub_pc_0, SubA, SubA);
1368  CHKERR PCSetType(sub_pc_0, PCFIELDSPLIT);
1369  CHKERR PCFieldSplitSetIS(sub_pc_0, NULL, sub_nested_is_rows[0]);
1370  CHKERR PCFieldSplitSetIS(sub_pc_0, NULL, sub_nested_is_rows[1]);
1371  CHKERR PCFieldSplitSetType(sub_pc_0, PC_COMPOSITE_MULTIPLICATIVE);
1372  // CHKERR
1373  // PCFieldSplitSetSchurFactType(sub_pc_0,PC_FIELDSPLIT_SCHUR_FACT_LOWER);
1374  // CHKERR PCFieldSplitSetType(sub_pc_0,PC_COMPOSITE_SCHUR);
1375  CHKERR PCSetUp(sub_pc_0);
1376  }
1377  } else {
1378  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1379  "Only works with pre-conditioner PCFIELDSPLIT");
1380  }
1381  CHKERR KSPSetUp(solver);
1382  }
1383 
1384  // Solve system of equations
1385  CHKERR KSPSolve(solver, F, D);
1386 
1387  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1388  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1389  CHKERR DMoFEMMeshToLocalVector(dm_control, D, INSERT_VALUES,
1390  SCATTER_REVERSE);
1391 
1392  if (debug) {
1393  CHKERR VecView(D, PETSC_VIEWER_DRAW_WORLD);
1394  std::string wait;
1395  std::cin >> wait;
1396  }
1397 
1398  // Clean sub matrices and sub indices
1399  for (int i = 0; i != 2; i++) {
1400  if (sub_nested_is_rows[i]) {
1401  CHKERR ISDestroy(&sub_nested_is_rows[i]);
1402  }
1403  if (sub_nested_is_cols[i]) {
1404  CHKERR ISDestroy(&sub_nested_is_cols[i]);
1405  }
1406  for (int j = 0; j != 2; j++) {
1407  if (sub_nested_matrices(i, j)) {
1408  CHKERR MatDestroy(&sub_nested_matrices(i, j));
1409  }
1410  }
1411  }
1412  for (int i = 0; i != 2; i++) {
1413  if (nested_is_rows[i]) {
1414  CHKERR ISDestroy(&nested_is_rows[i]);
1415  }
1416  if (nested_is_cols[i]) {
1417  CHKERR ISDestroy(&nested_is_cols[i]);
1418  }
1419  for (int j = 0; j != 2; j++) {
1420  if (nested_matrices(i, j)) {
1421  CHKERR MatDestroy(&nested_matrices(i, j));
1422  }
1423  }
1424  }
1425 
1426  CHKERR MatDestroy(&SubA);
1427  CHKERR MatDestroy(&A);
1428  CHKERR VecDestroy(&D);
1429  CHKERR VecDestroy(&F);
1430 
1431  CHKERR DMDestroy(&dm_sub_volume_control);
1432 
1433  PostProcVolumeOnRefinedMesh post_proc(m_field);
1434  {
1435  CHKERR post_proc.generateReferenceElementMesh();
1436  CHKERR post_proc.addFieldValuesPostProc("U");
1437  CHKERR post_proc.addFieldValuesGradientPostProc("U");
1438  // add postpocessing for sresses
1439  post_proc.getOpPtrVector().push_back(new PostProcHookStress(
1440  m_field, post_proc.postProcMesh, post_proc.mapGaussPts, "U",
1441  post_proc.commonData, &elastic.setOfBlocks));
1442  CHKERR DMoFEMLoopFiniteElements(dm_control, "ELASTIC", &post_proc);
1443  CHKERR post_proc.writeFile("out.h5m");
1444  elastic.getLoopFeEnergy().snes_ctx = SnesMethod::CTX_SNESNONE;
1445  elastic.getLoopFeEnergy().eNergy = 0;
1446  CHKERR DMoFEMLoopFiniteElements(dm_control, "ELASTIC",
1447  &elastic.getLoopFeEnergy());
1448  PetscPrintf(PETSC_COMM_WORLD, "Elastic energy %6.4e\n",
1449  elastic.getLoopFeEnergy().eNergy);
1450  }
1451 
1452  {
1453  PostProcFaceOnRefinedMesh post_proc_face(m_field);
1454  CHKERR post_proc_face.generateReferenceElementMesh();
1455  post_proc_face.getOpPtrVector().push_back(
1456  new OpCalculateInvJacForFace(inv_jac));
1457  if (is_curl) {
1458  post_proc_face.getOpPtrVector().push_back(
1459  new OpSetInvJacHcurlFace(inv_jac));
1460  post_proc_face.getOpPtrVector().push_back(
1461  new OpVirtualCurlRho("RHO", common_data));
1462  } else {
1463  post_proc_face.getOpPtrVector().push_back(
1464  new OpSetInvJacH1ForFace(inv_jac));
1465  post_proc_face.getOpPtrVector().push_back(
1466  new OpVirtualPotentialRho("RHO", common_data));
1467  }
1468  post_proc_face.getOpPtrVector().push_back(
1469  new PostProcTraction(m_field, post_proc_face.postProcMesh,
1470  post_proc_face.mapGaussPts, common_data));
1471  CHKERR DMoFEMLoopFiniteElements(dm_control, "D", &post_proc_face);
1472  CHKERR post_proc_face.writeFile("out_tractions.h5m");
1473  }
1474 
1475  CHKERR DMDestroy(&dm_control);
1476 
1477  }
1478  CATCH_ERRORS;
1479 
1481  return 0;
1482 }
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string &name_row)=0
set field col which finite element use
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode 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
Hook equation.
Definition: Hooke.hpp:31
boost::shared_ptr< SubProblemData > & getSubData() const
Get main problem of sub-problem is.
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:434
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:367
static char help[]
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[], EntityType lo_type=MBVERTEX, EntityType hi_type=MBMAXTYPE)
Definition: DMMMoFEM.cpp:210
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:414
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
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
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.
Core (interface) class.
Definition: Core.hpp:50
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:496
char mesh_file_name[255]
MoFEMErrorCode getDofByNameEntAndEntDofIdx(const string name, 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 Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
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:166
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:187
Operator post-procesing stresses for Hook isotropic material.
DEPRECATED MoFEMErrorCode synchronise_entities(Range &ent, int verb=DEFAULT_VERBOSITY)
virtual int get_comm_rank() const =0
Calculate inverse of jacobian for face element.
#define _IT_NUMEREDDOF_ROW_BY_NAME_ENT_PART_FOR_LOOP_(PROBLEMPTR, NAME, ENT, PART, IT)
use with loops to iterate row DOFs
MoFEMErrorCode printForceSet() const
print meshsets with force boundary conditions data structure
virtual MoFEMErrorCode get_problem(const std::string &problem_name, const Problem **problem_ptr) const =0
Get problem database (data structure)
CHKERRQ(ierr)
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:43
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:53
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:146
MoFEMErrorCode printDisplacementSet() const
print meshsets with displacement boundary conditions data structure
Calculate inverse of jacobian for face element.
Post processing.
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMMoFEM.cpp:457
MoFEMErrorCode printMaterialsSet() const
print meshsets with material data structure set on it
field with continuous tangents
Definition: definitions.h:172
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition: DMMMoFEM.cpp:467
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:109
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:596
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:285
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
Post-process tractions.
Definition: CellForces.hpp:528
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
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:433
DEPRECATED MoFEMErrorCode synchronise_field_entities(const std::string &name, int verb=DEFAULT_VERBOSITY)
continuous field
Definition: definitions.h:171
Calculate and assemble B matrix.
Definition: CellForces.hpp:140
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition: DMMMoFEM.cpp:971
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...
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:61
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.