rve_fibre_directions_interface.cpp File Reference

#include <gm_rule.h>
#include <MoFEM.hpp>
#include <boost/math/special_functions/round.hpp>
#include <MethodForForceScaling.hpp>
#include <DirichletBC.hpp>
#include <Projection10NodeCoordsOnField.hpp>
#include <petsctime.h>
#include "SurfacePressure.hpp"

Go to the source code of this file.

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

Variables
static char	help [] = "...\n\n"

Function Documentation

main()

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

Getting No. of Fibres and their index to be used for Potential Flow Problem

Solve problem potential flow problem for fibres one by one

Definition at line 46 of file rve_fibre_directions_interface.cpp.

                                 {
 
    MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
 
    try {
 
    moab::Core mb_instance;
    moab::Interface& moab = mb_instance;
    int rank;
    MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 
    PetscBool flg = PETSC_TRUE;
    char mesh_file_name[255];
    ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
    if(flg != PETSC_TRUE) {
      SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
    }
    PetscInt order;
    ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
    if(flg != PETSC_TRUE) {
      order = 1;
    }
//    cout<<" order "<<order<<endl;
 
 
    PetscInt mesh_refinement_level;
    ierr = PetscOptionsGetInt(PETSC_NULL,"-my_mesh_ref_level",&mesh_refinement_level,&flg); CHKERRQ(ierr);
    if(flg != PETSC_TRUE) {
      mesh_refinement_level = 0;
    }
 
 
    const char *option;
    option = "";//"PARALLEL=BCAST;";//;DEBUG_IO";
    rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
    ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
    if(pcomm == NULL) pcomm =  new ParallelComm(&moab,PETSC_COMM_WORLD);
 
 
    MoFEM::Core core(moab);
    MoFEM::Interface& m_field = core;
    PrismInterface *interface_ptr;
    ierr = m_field.query_interface(interface_ptr); CHKERRQ(ierr);
 
    //=======================================================================================================
    //Seting nodal coordinates on the surface to make sure they are periodic
    //=======================================================================================================
 
    Range SurTrisNeg, SurTrisPos;
    ierr = m_field.get_cubit_msId_entities_by_dimension(101,SIDESET,2,SurTrisNeg,true); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 101 = %d\n",SurTrisNeg.size()); CHKERRQ(ierr);
    ierr = m_field.get_cubit_msId_entities_by_dimension(102,SIDESET,2,SurTrisPos,true); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 102 = %d\n",SurTrisPos.size()); CHKERRQ(ierr);
 
    Range SurNodesNeg,SurNodesPos;
    rval = moab.get_connectivity(SurTrisNeg,SurNodesNeg,true); CHKERRQ_MOAB(rval);
    cout<<" All nodes on negative surfaces " << SurNodesNeg.size()<<endl;
    rval = moab.get_connectivity(SurTrisPos,SurNodesPos,true); CHKERRQ_MOAB(rval);
    cout<<" All nodes on positive surfaces " << SurNodesPos.size()<<endl;
 
 
    double roundfact=1e3;   double coords_nodes[3];
    for(Range::iterator nit = SurNodesNeg.begin(); nit!=SurNodesNeg.end();  nit++) {
      rval = moab.get_coords(&*nit,1,coords_nodes);  CHKERRQ_MOAB(rval);
      //round values to 3 disimal places
      coords_nodes[0]=round(coords_nodes[0]*roundfact)/roundfact;
      coords_nodes[1]=round(coords_nodes[1]*roundfact)/roundfact;
      coords_nodes[2]=round(coords_nodes[2]*roundfact)/roundfact;
      rval = moab.set_coords(&*nit,1,coords_nodes);  CHKERRQ_MOAB(rval);
//      cout<<"   coords_nodes[0]= "<<coords_nodes[0] << "   coords_nodes[1]= "<< coords_nodes[1] << "   coords_nodes[2]= "<< coords_nodes[2] <<endl;
    }
 
    for(Range::iterator nit = SurNodesPos.begin(); nit!=SurNodesPos.end();  nit++) {
      rval = moab.get_coords(&*nit,1,coords_nodes);  CHKERRQ_MOAB(rval);
      //round values to 3 disimal places
 
      coords_nodes[0]=round(coords_nodes[0]*roundfact)/roundfact;
      coords_nodes[1]=round(coords_nodes[1]*roundfact)/roundfact;
      coords_nodes[2]=round(coords_nodes[2]*roundfact)/roundfact;
      rval = moab.set_coords(&*nit,1,coords_nodes);  CHKERRQ_MOAB(rval);
//      cout<<"   coords_nodes[0]= "<<coords_nodes[0] << "   coords_nodes[1]= "<< coords_nodes[1] << "   coords_nodes[2]= "<< coords_nodes[2] <<endl;
    }
    //=======================================================================================================
 
//    string wait;
//    cin>>wait;
 
    //add fields
    ierr = m_field.add_field("POTENTIAL_FIELD",H1,AINSWORTH_LEGENDRE_BASE,1); CHKERRQ(ierr);
    ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
 
    ///Getting No. of Fibres and their index to be used for Potential Flow Problem
    int noOfFibres=0;
    for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,BLOCKSET|UNKNOWNNAME,it)) {
 
      std::size_t found=it->getName().find("PotentialFlow");
      if (found==std::string::npos) continue;
      noOfFibres += 1;
    }
        cout<<"No. of Fibres for Potential Flow : "<<noOfFibres<<endl;
 
    vector<int> fibreList(noOfFibres,0);
    int aaa=0;
 
    for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,BLOCKSET|UNKNOWNNAME,it)) {
 
      std::size_t interfaceFound=it->getName().find("PotentialFlow_");
      if (interfaceFound==std::string::npos) continue;
 
      std::string str2 = it->getName().substr (14,50);
 
      fibreList[aaa] = atoi(str2.c_str()); //atoi converet string to integer
//      cout<<"atoi(str2.c_str()) = " <<atoi(str2.c_str())<<endl;
      aaa += 1;
    }
 
    //add finite elements
    ierr = m_field.add_finite_element( "POTENTIAL_ELEM" ); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_row("POTENTIAL_ELEM" ,"POTENTIAL_FIELD"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_col("POTENTIAL_ELEM" ,"POTENTIAL_FIELD"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_data("POTENTIAL_ELEM" ,"POTENTIAL_FIELD"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_data("POTENTIAL_ELEM" ,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
 
    //add problems
    ierr = m_field.add_problem( "POTENTIAL_PROBLEM" ); CHKERRQ(ierr);
    //define problems and finite elements
    ierr = m_field.modify_problem_add_finite_element( "POTENTIAL_PROBLEM" , "POTENTIAL_ELEM" ); CHKERRQ(ierr);
 
        Tag th_meshset_info;
        int def_meshset_info[2] = {0,0};
        rval = moab.tag_get_handle("MESHSET_INFO",2,MB_TYPE_INTEGER,th_meshset_info,MB_TAG_CREAT|MB_TAG_SPARSE,&def_meshset_info);
 
        int meshset_data[2];
    EntityHandle root = moab.get_root_set();
        rval = moab.tag_get_data(th_meshset_info,&root,1,meshset_data); CHKERRQ_MOAB(rval);
 
    ierr = m_field.seed_ref_level_3D(0,BitRefLevel().set(0)); CHKERRQ(ierr);
    vector<BitRefLevel> bit_levels;
    bit_levels.push_back(BitRefLevel().set(0));
 
        //MESH REFINEMENT
        int ll = 1;
 
 
        //*****INTERFACE INSERTION******
    for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,SIDESET,cit)) {
 
      std::size_t interfaceFound=cit->getName().find("interface");
      if (interfaceFound==std::string::npos) continue;
 
      ierr = PetscPrintf(PETSC_COMM_WORLD,"Insert Interface %d\n",cit->getMeshsetId()); CHKERRQ(ierr);
      EntityHandle cubit_meshset = cit->getMeshset();
      int meshset_data_int[2];
      rval = moab.tag_get_data(th_meshset_info,&cubit_meshset,1,meshset_data_int); CHKERRQ_MOAB(rval);
      if(meshset_data_int[1]==0){
        //get tet enties form back bit_level
        EntityHandle ref_level_meshset = 0;
        rval = moab.create_meshset(MESHSET_SET,ref_level_meshset); CHKERRQ_MOAB(rval);
        ierr = m_field.get_entities_by_type_and_ref_level(bit_levels.back(),BitRefLevel().set(),MBTET,ref_level_meshset); CHKERRQ(ierr);
        ierr = m_field.get_entities_by_type_and_ref_level(bit_levels.back(),BitRefLevel().set(),MBPRISM,ref_level_meshset); CHKERRQ(ierr);
        Range ref_level_tets;
        rval = moab.get_entities_by_handle(ref_level_meshset,ref_level_tets,true); CHKERRQ_MOAB(rval);
        //get faces and test to split
        ierr = interface_ptr->getSides(cubit_meshset,bit_levels.back(),true); CHKERRQ(ierr);
        //set new bit level
        bit_levels.push_back(BitRefLevel().set(ll++));
        //split faces and edges
        ierr = interface_ptr->splitSides(ref_level_meshset,bit_levels.back(),cubit_meshset,true,true,0); CHKERRQ(ierr);
        //clean meshsets
        rval = moab.delete_entities(&ref_level_meshset,1); CHKERRQ_MOAB(rval);
        int meshset_data_ins[2] = {ll,1};
        rval = moab.tag_set_data(th_meshset_info,&cubit_meshset,1,meshset_data_ins); CHKERRQ_MOAB(rval);
      }
      //update cubit meshsets
      for(_IT_CUBITMESHSETS_FOR_LOOP_(m_field,ciit)) {
        EntityHandle cubit_meshset = ciit->meshset;
        ierr = m_field.update_meshset_by_entities_children(cubit_meshset,bit_levels.back(),cubit_meshset,MBVERTEX,true); CHKERRQ(ierr);
        ierr = m_field.update_meshset_by_entities_children(cubit_meshset,bit_levels.back(),cubit_meshset,MBEDGE,true); CHKERRQ(ierr);
        ierr = m_field.update_meshset_by_entities_children(cubit_meshset,bit_levels.back(),cubit_meshset,MBTRI,true); CHKERRQ(ierr);
        ierr = m_field.update_meshset_by_entities_children(cubit_meshset,bit_levels.back(),cubit_meshset,MBTET,true); CHKERRQ(ierr);
      }
    }
 
 
        //End of refinement, save level of refinement
        int meshset_data_root[2]={ll,0};
        rval = moab.tag_set_data(th_meshset_info,&root,1,meshset_data_root); CHKERRQ_MOAB(rval);
 
        /******************TETS TO MESHSET AND SAVING TETS ENTITIES******************/
        EntityHandle out_tet_meshset;
        rval = moab.create_meshset(MESHSET_SET,out_tet_meshset); CHKERRQ_MOAB(rval);
        ierr = m_field.get_entities_by_type_and_ref_level(bit_levels.back(),BitRefLevel().set(),MBTET,out_tet_meshset); CHKERRQ(ierr);
        rval = moab.write_file("out_tets.vtk","VTK","",&out_tet_meshset,1); CHKERRQ_MOAB(rval);
        /*******************************************************/
 
        /******************PRISMS TO MESHSET AND SAVING PRISMS ENTITIES******************/
        EntityHandle out_meshset_prism;
        rval = moab.create_meshset(MESHSET_SET,out_meshset_prism); CHKERRQ_MOAB(rval);
        ierr = m_field.get_entities_by_type_and_ref_level(bit_levels.back(),BitRefLevel().set(),MBPRISM,out_meshset_prism); CHKERRQ(ierr);
        rval = moab.write_file("out_prism.vtk","VTK","",&out_meshset_prism,1); CHKERRQ_MOAB(rval);
        /*******************************************************/
 
        EntityHandle out_meshset;
        rval = moab.create_meshset(MESHSET_SET,out_meshset); CHKERRQ_MOAB(rval);
        ierr = m_field.get_entities_by_ref_level(bit_levels.back(),BitRefLevel().set(),out_meshset); CHKERRQ(ierr);
        rval = moab.write_file("out_all_mesh.vtk","VTK","",&out_meshset,1); CHKERRQ_MOAB(rval);
 
    Range LatestRefinedTets;
        rval = moab.get_entities_by_type(out_meshset, MBTET,LatestRefinedTets,true); CHKERRQ_MOAB(rval);
 
    Range LatestRefinedTris;
        rval = moab.get_entities_by_type(out_meshset, MBTRI,LatestRefinedTris,true); CHKERRQ_MOAB(rval);
 
 
    BitRefLevel problem_bit_level = bit_levels.back();
 
    //add entities to field
    ierr = m_field.add_ents_to_field_by_type(0,MBTET,"POTENTIAL_FIELD"); CHKERRQ(ierr);
    ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
 
 
    ierr = m_field.set_field_order(0,MBVERTEX,"POTENTIAL_FIELD",1); CHKERRQ(ierr);
    ierr = m_field.set_field_order(0,MBEDGE,"POTENTIAL_FIELD",order); CHKERRQ(ierr);
    ierr = m_field.set_field_order(0,MBTRI,"POTENTIAL_FIELD",order); CHKERRQ(ierr);
    ierr = m_field.set_field_order(0,MBTET,"POTENTIAL_FIELD",order); CHKERRQ(ierr);
 
    ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
    ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
    ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",order>1 ? 2 : 1); CHKERRQ(ierr);
    ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
 
    //define flux element
    ierr = m_field.add_finite_element("FLUX_FE" ,MF_ZERO); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_row("FLUX_FE","POTENTIAL_FIELD"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_col("FLUX_FE","POTENTIAL_FIELD"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_data("FLUX_FE","POTENTIAL_FIELD"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_data( "FLUX_FE" ,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
    ierr = m_field.modify_problem_add_finite_element("POTENTIAL_PROBLEM", "FLUX_FE" ); CHKERRQ(ierr);
 
 
    //create matrices and vectors
        Vec F;
        Vec D;
        Mat A;
 
    EntityHandle out_meshset_fibres;
    if(pcomm->rank()==0) {
      rval = moab.create_meshset(MESHSET_SET,out_meshset_fibres); CHKERRQ_MOAB(rval);
    }
 
    ///Solve problem potential flow problem for fibres one by one
    for (int cc = 0; cc < noOfFibres; cc++) {
        ostringstream rrr, ppp1, ppp2;
        rrr << "PotentialFlow_" << fibreList[cc];
        Range new_tets;
        for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,BLOCKSET|UNKNOWNNAME,it)) {
          if(it->getName() ==  rrr.str().c_str() ) {
            //set problem level
            Range TetsInBlock;
            rval = moab.get_entities_by_type(it->meshset, MBTET,TetsInBlock,true); CHKERRQ_MOAB(rval);
            new_tets = intersect(LatestRefinedTets,TetsInBlock);
            //add finite elements entities
            ierr = m_field.add_ents_to_finite_element_by_type(new_tets,MBTET,"POTENTIAL_ELEM"); CHKERRQ(ierr);
          }
        }
 
        Range new_tris_1, new_tris_2;
        //flux boundary conditions
        ppp1 << "PressureIO_" << fibreList[cc] <<"_1";
        ppp2 << "PressureIO_" << fibreList[cc] <<"_2";
        for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,SIDESET|PRESSURESET,it)) {
          if (ppp1.str().c_str()==it->getName()){
            Range tris;
            rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
            new_tris_1 = intersect(LatestRefinedTris,tris);
            ierr = m_field.add_ents_to_finite_element_by_type(new_tris_1,MBTRI, "FLUX_FE" ); CHKERRQ(ierr);
            }
 
 
          if (ppp2.str().c_str()==it->getName()){
            Range tris;
            rval = m_field.get_moab().get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
            new_tris_2 = intersect(LatestRefinedTris,tris);
            ierr = m_field.add_ents_to_finite_element_by_type(new_tris_2,MBTRI, "FLUX_FE" ); CHKERRQ(ierr);
            }
         }
 
        //set problem level
        ierr = m_field.modify_problem_ref_level_add_bit("POTENTIAL_PROBLEM" ,problem_bit_level); CHKERRQ(ierr);
 
        //build fields
        ierr = m_field.build_fields(); CHKERRQ(ierr);
        //build finite elements
        ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
        //build adjacencies
        ierr = m_field.build_adjacencies(problem_bit_level); CHKERRQ(ierr);
        //build problem
        ierr = m_field.build_problems(); CHKERRQ(ierr);
 
        ierr = m_field.partition_problem("POTENTIAL_PROBLEM" ); CHKERRQ(ierr);
        ierr = m_field.partition_finite_elements( "POTENTIAL_PROBLEM" ); CHKERRQ(ierr);
        ierr = m_field.partition_ghost_dofs( "POTENTIAL_PROBLEM" ); CHKERRQ(ierr);
 
        ierr = m_field.VecCreateGhost("POTENTIAL_PROBLEM" ,ROW,&F); CHKERRQ(ierr);
        ierr = m_field.VecCreateGhost("POTENTIAL_PROBLEM" ,COL,&D); CHKERRQ(ierr);
        ierr = m_field.MatCreateMPIAIJWithArrays("POTENTIAL_PROBLEM" ,&A); CHKERRQ(ierr);
 
        Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
        ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
 
        ostringstream zero_pressure;
        zero_pressure << "ZeroPressure_" << fibreList[cc];
       //get nodes and other entities to fix
        Range fix_nodes;
        for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,NODESET|UNKNOWNNAME,it)) {
          if (zero_pressure.str().c_str()==it->getName()){
            rval = moab.get_entities_by_type(it->meshset,MBVERTEX,fix_nodes,true); CHKERRQ_MOAB(rval);
            Range edges;
            rval = moab.get_entities_by_type(it->meshset,MBEDGE,edges,true); CHKERRQ_MOAB(rval);
            Range tris;
            rval = moab.get_entities_by_type(it->meshset,MBTRI,tris,true); CHKERRQ_MOAB(rval);
            Range adj;
            rval = moab.get_connectivity(tris,adj,true); CHKERRQ_MOAB(rval);
            fix_nodes.insert(adj.begin(),adj.end());
            rval = moab.get_connectivity(edges,adj,true); CHKERRQ_MOAB(rval);
            fix_nodes.insert(adj.begin(),adj.end());
            rval = moab.get_adjacencies(tris,1,false,edges,moab::Interface::UNION); CHKERRQ_MOAB(rval);
          }
        }
        FixBcAtEntities fix_dofs(m_field,"POTENTIAL_FIELD",A,D,F,fix_nodes);
        //initialize data structure
        ierr = m_field.problem_basic_method_preProcess("POTENTIAL_PROBLEM",fix_dofs); CHKERRQ(ierr);
 
        ierr = VecZeroEntries(F); CHKERRQ(ierr);
        ierr = VecGhostUpdateBegin(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
        ierr = VecGhostUpdateEnd(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
        ierr = MatZeroEntries(A); CHKERRQ(ierr);
 
        //neuman flux bc elements
        //surface forces
        boost::ptr_map<string,NeummanForcesSurface> neumann_forces;
        string fe_name_str = "FLUX_FE";
        neumann_forces.insert(fe_name_str,new NeummanForcesSurface(m_field));
        for(_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(m_field,SIDESET|PRESSURESET,it)) {
          if (ppp1.str().c_str()==it->getName() || ppp2.str().c_str()==it->getName()){
            bool ho_geometry = m_field.check_field("MESH_NODE_POSITIONS");
            ierr = neumann_forces.at("FLUX_FE").addFlux("POTENTIAL_FIELD",F,it->getMeshsetId(),ho_geometry); CHKERRQ(ierr);
          }
        }
        boost::ptr_map<string,NeummanForcesSurface>::iterator mit = neumann_forces.begin();
        for(;mit!=neumann_forces.end();mit++) {
          ierr = m_field.loop_finite_elements("POTENTIAL_PROBLEM", mit->first, mit->second->getLoopFe()); CHKERRQ(ierr);
        }
 
        LaplacianElem elem(m_field,A,F);
        ierr = m_field.loop_finite_elements("POTENTIAL_PROBLEM", "POTENTIAL_ELEM", elem);  CHKERRQ(ierr);
 
        //post proces fix boundary conditiond
        ierr = m_field.problem_basic_method_postProcess("POTENTIAL_PROBLEM", fix_dofs); CHKERRQ(ierr);
 
        ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
        ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
 
        ierr = VecAssemblyBegin(F); CHKERRQ(ierr);
        ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
 
        //set matrix possitives define and symetric for cholesky and icc preceonditionser
        ierr = MatSetOption(A,MAT_SPD,PETSC_TRUE); CHKERRQ(ierr);
 
        KSP solver;
        ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
        ierr = KSPSetOperators(solver,A,A); CHKERRQ(ierr);
        ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
        ierr = KSPSetUp(solver); CHKERRQ(ierr);
 
        ierr = KSPSolve(solver,F,D); CHKERRQ(ierr);
        ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
        ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
        ierr = m_field.set_global_ghost_vector("POTENTIAL_PROBLEM",ROW,D,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
 
        if(pcomm->rank()==0) {
          ierr = m_field.get_problem_finite_elements_entities("POTENTIAL_PROBLEM","POTENTIAL_ELEM",out_meshset_fibres); CHKERRQ(ierr);
        }
 
        //remove intities from finite elements (to prepare it for next fibres)
        ierr = m_field.remove_ents_from_finite_element("POTENTIAL_ELEM", new_tets);      CHKERRQ(ierr);
        ierr = m_field.remove_ents_from_finite_element("FLUX_FE", new_tris_1);      CHKERRQ(ierr);
        ierr = m_field.remove_ents_from_finite_element("FLUX_FE", new_tris_2);      CHKERRQ(ierr);
        //clear problems (to prepare it for next fibre)
        ierr = m_field.clear_problems(); CHKERRQ(ierr);
 
        ierr = VecDestroy(&F); CHKERRQ(ierr);
        ierr = VecDestroy(&D); CHKERRQ(ierr);
        ierr = MatDestroy(&A); CHKERRQ(ierr);
        ierr = KSPDestroy(&solver); CHKERRQ(ierr);
     }
 
 
    Tag th_phi;
    double def_val = 0;
    rval = moab.tag_get_handle("PHI",1,MB_TYPE_DOUBLE,th_phi,MB_TAG_CREAT|MB_TAG_SPARSE,&def_val); CHKERRQ_MOAB(rval);
    for(_IT_GET_DOFS_FIELD_BY_NAME_FOR_LOOP_(m_field,"POTENTIAL_FIELD",dof)) {
      EntityHandle ent = dof->get()->getEnt();
      double val = dof->get()->getFieldData();
      rval = moab.tag_set_data(th_phi,&ent,1,&val); CHKERRQ_MOAB(rval);
    }
 
    ProjectionFieldOn10NodeTet ent_method_phi_on_10nodeTet(m_field,"POTENTIAL_FIELD",true,false,"PHI");
    ierr = m_field.loop_dofs("POTENTIAL_FIELD",ent_method_phi_on_10nodeTet); CHKERRQ(ierr);
    ent_method_phi_on_10nodeTet.setNodes = false;
    ierr = m_field.loop_dofs("POTENTIAL_FIELD",ent_method_phi_on_10nodeTet); CHKERRQ(ierr);
 
    if(pcomm->rank()==0) {
      rval = moab.write_file("solution_RVE.h5m"); CHKERRQ_MOAB(rval);
    }
 
    if(pcomm->rank()==0) {
      rval = moab.write_file("out_potential_flow.vtk","VTK","",&out_meshset_fibres,1); CHKERRQ_MOAB(rval);
    }
 
    }
    CATCH_ERRORS;
 
    MoFEM::Core::Finalize();
}

Variable Documentation

help

char help[] = "...\n\n"

static

Definition at line 44 of file rve_fibre_directions_interface.cpp.