rve_fibre_directions_interface.cpp

Go to the documentation of this file.

/* This file is part of MoFEM.
 * MoFEM is free software: you can redistribute it and/or modify it under
 * the terms of the GNU Lesser General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
 * License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
 
//#include "FunctionsForFieldData.hpp"
//#include "cholesky.hpp"
 
extern "C" {
#include <gm_rule.h>
}
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
 
#include <boost/math/special_functions/round.hpp>
#include <MethodForForceScaling.hpp>
#include <DirichletBC.hpp>
 
#include <Projection10NodeCoordsOnField.hpp>
#include <petsctime.h>
 
#include <FEMethod_LowLevelStudent.hpp>
#include <FEMethod_UpLevelStudent.hpp>
 
#include <PostProcVertexMethod.hpp>
#include <PostProcDisplacementAndStrainOnRefindedMesh.hpp>
 
using namespace ObosleteUsersModules;
 
#include <MethodForForceScaling.hpp>
#include "PotentialFlowFEMethod.hpp"
#include "SurfacePressure.hpp"
 
 
static char help[] = "...\n\n";
 
int main(int argc, char *argv[]) {
 
    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
    //Tpressure_sideset_id will comes form salome-meca for fibre_1
    //it is 100101 for neg pressure and 100102 for positive pressure
    int pressure_sideset_id = 10000;
    for (int cc = 0; cc < noOfFibres; cc++) {
        pressure_sideset_id=pressure_sideset_id+100; //increment the fibre
        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() || it->getMeshsetId() == pressure_sideset_id+1){
            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() || it->getMeshsetId() == pressure_sideset_id+2){
            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()
              || it->getMeshsetId()==pressure_sideset_id+1
              || it->getMeshsetId()==pressure_sideset_id+2){
            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();
}