mass_calculation.cpp

Go to the documentation of this file.

/* This file is part of MoFEM.
* \ingroup bone_remodelling
 
* 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 <MoFEM.hpp>
using namespace MoFEM;
#include <Projection10NodeCoordsOnField.hpp>
 
#ifdef WITH_METAIO
#include <metaImage.h>
#endif
 
#include <moab/Skinner.hpp>
#include <moab/AdaptiveKDTree.hpp>
 
#include <cholesky.hpp>
#include <DensityMaps.hpp>
using namespace BoneRemodeling;
 
#include <PostProcOnRefMesh.hpp>
 
static char help[] = "\n";
 
 
 
#ifndef WITH_METAIO
  #error "You need compile users modules with MetaIO -DWITH_METAIO=1"
#endif
 
int main(int argc, char *argv[]) {
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
 
  try {
 
    moab::Core mb_instance;
    moab::Interface& moab = mb_instance;
 
    char mesh_file_name[255];
    PetscBool flg_file = PETSC_FALSE;
    char meta_file_name[255];
    PetscBool meta_flg_file = PETSC_FALSE;
    PetscInt order = 1;
    PetscBool is_atom_test = PETSC_FALSE;
    // char approximation_base[255];
 
    ierr = PetscOptionsBegin(PETSC_COMM_WORLD,"","Bone remodeling","none"); CHKERRQ(ierr);
    ierr = PetscOptionsString(
      "-my_file",
      "mesh file name","",
      "mesh.h5m",mesh_file_name,255,&flg_file
    ); CHKERRQ(ierr);
 
    ierr = PetscOptionsString(
      "-meta_image",
      "meta image name","",
      "file.mhd",meta_file_name,255,&meta_flg_file
    ); CHKERRQ(ierr);
 
    // ierr = PetscOptionsString(
    //   "-base",
    //   "approximation base name","",
    //   "",approximation_base,255,PETSC_NULL
    // ); CHKERRQ(ierr);
 
    ierr = PetscOptionsBool(
      "-my_is_atom_test",
      "is used with testing, exit with error when diverged","",
      PETSC_FALSE,&is_atom_test,PETSC_NULL
    ); CHKERRQ(ierr);
 
    ierr = PetscOptionsInt(
      "-my_order",
      "default approximation order","",
      1,&order,PETSC_NULL
    ); CHKERRQ(ierr);
 
    ierr = PetscOptionsEnd(); CHKERRQ(ierr);
 
    // #ifdef WITH_METAIO
 
    MetaImage meta_file(meta_file_name);
    cout << "Quantity " << meta_file.Quantity() << endl;
    const int *dim_size = meta_file.DimSize();
    cout << "Image dimension: " << dim_size[0] << " " << dim_size[1] << " " << dim_size[2] << endl;
    const double *elemet_size = meta_file.ElementSize();
    cout << "Image element size: " << elemet_size[0] << " " << elemet_size[1] << " " << elemet_size[2] << endl;
    const double *offset = meta_file.Offset();
    cout << "Image offset: " << offset[0] << " " << offset[1] << " " << offset[2] << endl;
    const double * transform_matrix = meta_file.TransformMatrix();
    cout << "Image Transform Matrix:\n";
    for(int rr = 0;rr<3;rr++) {
      for(int cc = 0;cc<3;cc++) {
        cout << transform_matrix[rr*3+cc] << " ";
      }
      cout << endl;
    }
 
    //Read parameters from line command
    if(flg_file != PETSC_TRUE) {
      SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
    }
    const char *option;
    option = "";
    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;
    //meshset consisting all entities in mesh
    EntityHandle root_set = moab.get_root_set();
    // Seed all mesh entities to MoFEM database, those entities can be potentially used as finite elements
    // or as entities which carry some approximation field.
    BitRefLevel bit_level0;
    bit_level0.set(0);
    ierr = m_field.seed_ref_level_3D(root_set,bit_level0); CHKERRQ(ierr);
    // Add field
    ierr = m_field.add_field("RHO",H1,AINSWORTH_LEGENDRE_BASE,1); CHKERRQ(ierr);
    ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3); CHKERRQ(ierr);
 
    // Add meshsets with material and boundary conditions
    auto *meshsets_manager_ptr = m_field.getInterface<MeshsetsManager>();
    ierr = meshsets_manager_ptr->setMeshsetFromFile(); CHKERRQ(ierr);
 
    PetscPrintf(PETSC_COMM_WORLD,"Read meshset. Added meshsets for bc.cfg\n");
    for(_IT_CUBITMESHSETS_FOR_LOOP_(m_field,it)) {
      PetscPrintf(
        PETSC_COMM_WORLD,
        "%s",static_cast<std::ostringstream&>(std::ostringstream().seekp(0) << *it << endl).str().c_str()
      );
      cerr << *it << endl;
    }
 
    Range tets_all;
    rval = m_field.get_moab().get_entities_by_type(0,MBTET,tets_all); CHKERRQ_MOAB(rval);
 
    for(_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field,BLOCKSET,it)) {
      string name = it->getName();
      if (name.compare(0,14,"NO_REMODELLING") == 0) {
        //nOremodellingBlock = true;
        EntityHandle meshset = it->getMeshset();
        Range tets_block;
        rval = m_field.get_moab().get_entities_by_type(
          meshset,MBTET,tets_block,true
        ); CHKERRQ_MOAB(rval);
        tets_all = subtract(tets_all,tets_block);
      }
    }
 
 
    // if ( strcmp(approximation_base, "lobatto") == 0){
    //   ierr = m_field.add_field("RHO",L2,AINSWORTH_LOBATTO_BASE,1); CHKERRQ(ierr);
    // } else if ( strcmp(approximation_base, "berstein") == 0){                 //not implemented yet
    //   ierr = m_field.add_field("RHO",L2,AINSOWRTH_BERNSTEIN_BEZIER_BASE,1); CHKERRQ(ierr);
    // }
    // else {
    //   ierr = m_field.add_field("RHO",L2,AINSWORTH_LEGENDRE_BASE,1); CHKERRQ(ierr);
    // }
    // Add entities to field (root_mesh, i.e. on all mesh entities fields are approx.)
    // On those entities and lower dimension entities approximation is spaned,
    ierr = m_field.add_ents_to_field_by_type(tets_all,MBTET,"RHO"); CHKERRQ(ierr);
    ierr = m_field.add_ents_to_field_by_type(root_set,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
 
 
    {
      ierr = m_field.add_field("DISTANCE_FROM_SURFACE",H1,AINSWORTH_LEGENDRE_BASE,1); CHKERRQ(ierr);
      Range vertices;
      rval = moab.get_entities_by_type(root_set,MBVERTEX,vertices); CHKERRQ_MOAB(rval);
      ierr = m_field.add_ents_to_field_by_type(vertices,MBTET,"DISTANCE_FROM_SURFACE"); CHKERRQ(ierr);
      Range tets;
      rval = moab.get_entities_by_type(root_set,MBTET,tets); CHKERRQ_MOAB(rval);
      Range edges;
      rval = moab.get_adjacencies(tets,1,true,edges,moab::Interface::UNION); CHKERRQ_MOAB(rval);
      ierr = m_field.add_ents_to_field_by_type(edges,MBEDGE,"DISTANCE_FROM_SURFACE"); CHKERRQ(ierr);
    }
 
    // Set approximation oder
    ierr = m_field.set_field_order(root_set,MBVERTEX,"RHO",1); CHKERRQ(ierr);
    ierr = m_field.set_field_order(root_set,MBEDGE,"RHO",order); CHKERRQ(ierr);
    ierr = m_field.set_field_order(root_set,MBTRI,"RHO",order); CHKERRQ(ierr);
    ierr = m_field.set_field_order(root_set,MBTET,"RHO",order); CHKERRQ(ierr);
    ierr = m_field.set_field_order(root_set,MBEDGE,"DISTANCE_FROM_SURFACE",2); CHKERRQ(ierr);
    ierr = m_field.set_field_order(root_set,MBVERTEX,"DISTANCE_FROM_SURFACE",1); CHKERRQ(ierr);
 
    // Apply 2nd order geometry approximation only on skin
    {
 
      Skinner skin(&moab);
      Range faces,tets;
      rval = moab.get_entities_by_type(0,MBTET,tets); CHKERRQ_MOAB(rval);
      rval = skin.find_skin(0,tets,false,faces); CHKERRQ_MOAB(rval);
      Range edges;
      rval = moab.get_adjacencies(
        faces,1,false,edges,moab::Interface::UNION
      ); CHKERRQ_MOAB(rval);
      ierr = m_field.set_field_order(edges,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
    }
    ierr = m_field.set_field_order(root_set,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
 
    ierr = m_field.build_fields(); CHKERRQ(ierr);
 
    {
      Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
      ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
    }
 
    SurfaceKDTree surface_distance(m_field);
    {
      Range tets;
      rval = moab.get_entities_by_type(root_set,MBTET,tets); CHKERRQ_MOAB(rval);
      ierr = surface_distance.takeASkin(tets); CHKERRQ(ierr);
      ierr = surface_distance.buildTree(); CHKERRQ(ierr);
      ierr = surface_distance.setDistanceFromSurface("DISTANCE_FROM_SURFACE"); CHKERRQ(ierr);
    }
 
    // Add finite element (this defines element declaration comes later)
    ierr = m_field.add_finite_element("DENSITY_MAP_ELEMENT"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_row("DENSITY_MAP_ELEMENT","RHO"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_col("DENSITY_MAP_ELEMENT","RHO"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_data("DENSITY_MAP_ELEMENT","RHO"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_data("DENSITY_MAP_ELEMENT","DISTANCE_FROM_SURFACE"); CHKERRQ(ierr);
    ierr = m_field.modify_finite_element_add_field_data("DENSITY_MAP_ELEMENT","MESH_NODE_POSITIONS"); CHKERRQ(ierr);
 
    // Add entities to that element
    ierr = m_field.add_ents_to_finite_element_by_type(tets_all,MBTET,"DENSITY_MAP_ELEMENT"); CHKERRQ(ierr);
 
    //build finite elements
    ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
    //build adjacencies between elements and degrees of freedom
    ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
 
    DM dm_rho;
    DMType dm_name = "DM_RHO";
    // Register DM problem
    ierr = DMRegister_MoFEM(dm_name); CHKERRQ(ierr);
    ierr = DMCreate(PETSC_COMM_WORLD,&dm_rho);CHKERRQ(ierr);
    ierr = DMSetType(dm_rho,dm_name);CHKERRQ(ierr);
    // Create DM instance
    ierr = DMMoFEMCreateMoFEM(dm_rho,&m_field,dm_name,bit_level0); CHKERRQ(ierr);
    // Configure DM form line command options (DM itself, solvers, pre-conditioners, ... )
    ierr = DMSetFromOptions(dm_rho); CHKERRQ(ierr);
    // Add elements to dm (only one here)
    ierr = DMMoFEMAddElement(dm_rho,"DENSITY_MAP_ELEMENT"); CHKERRQ(ierr);
    // Set up problem (number dofs, partition mesh, etc.)
    //ierr = DMSetUp(dm_rho); CHKERRQ(ierr);
 
    {
      ierr = DMSetUp(dm_rho); CHKERRQ(ierr);
      // Do serial matrix
      // const MoFEM::Problem *problem_ptr;
      // ierr = DMMoFEMGetProblemPtr(dm_rho,&problem_ptr); CHKERRQ(ierr);
      // ierr = m_field.build_problem(problem_ptr->getName()); CHKERRQ(ierr);
      // // FIXME: That is not optimal partitioning
      // ierr = m_field.partition_simple_problem(problem_ptr->getName()); CHKERRQ(ierr);
      // ierr = m_field.partition_finite_elements(problem_ptr->getName());
      // ierr = m_field.partition_ghost_dofs(problem_ptr->getName()); CHKERRQ(ierr);
    }
 
    Vec volume_vec;
    int volume_vec_ghost[] = { 0 };
    ierr = VecCreateGhost(
      PETSC_COMM_WORLD,(!m_field.get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
    );  CHKERRQ(ierr);
    ierr = VecZeroEntries(volume_vec); CHKERRQ(ierr);
    Vec mass_vec;
    ierr = VecDuplicate(volume_vec,&mass_vec); CHKERRQ(ierr);
    Vec mass_vec_approx;
    ierr = VecDuplicate(volume_vec,&mass_vec_approx); CHKERRQ(ierr);
 
    DensityMapFe fe_volume(m_field);
    auto &list_of_operators = fe_volume.getOpPtrVector();
    list_of_operators.push_back(new OpVolumeCalculation("RHO",volume_vec));
    DataFromMetaIO data_from_meta_io(meta_file);
    ierr = data_from_meta_io.fromOptions(); CHKERRQ(ierr);
    list_of_operators.push_back(new OpMassCalculation("RHO",mass_vec,data_from_meta_io));
    ierr = DMoFEMLoopFiniteElements(dm_rho,"DENSITY_MAP_ELEMENT",&fe_volume); CHKERRQ(ierr);
 
    ierr = VecAssemblyBegin(volume_vec); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(volume_vec); CHKERRQ(ierr);
    double volume;
    ierr = VecSum(volume_vec,&volume);  CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Volume %0.32g\n",volume); CHKERRQ(ierr);
 
    ierr = VecAssemblyBegin(mass_vec); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(mass_vec); CHKERRQ(ierr);
    double mass;
    ierr = VecSum(mass_vec,&mass);  CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Mass %0.32g\n",mass); CHKERRQ(ierr);
 
    ierr = VecDestroy(&volume_vec); CHKERRQ(ierr);
    ierr = VecDestroy(&mass_vec); CHKERRQ(ierr);
 
 
    Mat A;
    Vec F,D;
    ierr = DMCreateMatrix(dm_rho,&A); CHKERRQ(ierr);
    ierr = m_field.getInterface<VecManager>()->vecCreateGhost("DM_RHO",ROW,&F); CHKERRQ(ierr);
    ierr = m_field.getInterface<VecManager>()->vecCreateGhost("DM_RHO",COL,&D); CHKERRQ(ierr);
 
    ierr = MatZeroEntries(A); CHKERRQ(ierr);
    ierr = VecZeroEntries(F); CHKERRQ(ierr);
 
    CommonData common_data;
    common_data.globA = A;
    common_data.globF = F;
 
    common_data.distanceValue = boost::shared_ptr<VectorDouble>(new VectorDouble());
    common_data.distanceGradient = boost::shared_ptr<MatrixDouble>(new MatrixDouble());
    boost::shared_ptr<VectorDouble> dist_value = common_data.distanceValue;
    boost::shared_ptr<MatrixDouble> grad_value = common_data.distanceGradient;
 
    DensityMapFe fe_density_approximation(m_field);
 
    fe_density_approximation.getOpPtrVector().push_back
    (new OpCalculateScalarFieldValues("DISTANCE_FROM_SURFACE",dist_value));
    fe_density_approximation.getOpPtrVector().push_back
    (new OpCalculateScalarFieldGradient<3>("DISTANCE_FROM_SURFACE",grad_value));
    fe_density_approximation.getOpPtrVector().push_back
    (new OpCalculateLhs("RHO","RHO",common_data,data_from_meta_io));
    fe_density_approximation.getOpPtrVector().push_back
    (new OpCalulatefRhoAtGaussPts("RHO",common_data,data_from_meta_io,surface_distance));
    fe_density_approximation.getOpPtrVector().push_back
    (new OpAssmbleRhs("RHO",common_data));
    ierr = DMoFEMLoopFiniteElements(dm_rho,"DENSITY_MAP_ELEMENT",&fe_density_approximation); CHKERRQ(ierr);
 
    ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
    ierr = VecGhostUpdateBegin(F,ADD_VALUES,SCATTER_REVERSE);
    ierr = VecGhostUpdateEnd(F,ADD_VALUES,SCATTER_REVERSE);
    ierr = VecAssemblyBegin(F); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
 
    // MatView(A,PETSC_VIEWER_DRAW_WORLD);
    // std::string wait;
    // std::cin >> wait;
 
    {
      KSP ksp;
      ierr = KSPCreate(PETSC_COMM_WORLD,&ksp); CHKERRQ(ierr);
      ierr = KSPSetOperators(ksp,A,A); CHKERRQ(ierr);
      ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);
      ierr = KSPSolve(ksp,F,D); CHKERRQ(ierr); CHKERRQ(ierr);
    }
 
    ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = m_field.getInterface<VecManager>()->setGlobalGhostVector("DM_RHO",COL,D,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
 
    DensityMapFe fe_post_process_calulations(m_field);
 
    fe_post_process_calulations.getOpPtrVector().push_back(
      new OpMassCalculationFromApprox("RHO",mass_vec_approx)
    );
    ierr = DMoFEMLoopFiniteElements(dm_rho,"DENSITY_MAP_ELEMENT",&fe_post_process_calulations); CHKERRQ(ierr);
 
    ierr = VecAssemblyBegin(mass_vec_approx); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(mass_vec_approx); CHKERRQ(ierr);
    double mass_approx;
    ierr = VecSum(mass_vec_approx,&mass_approx);  CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Mass after approximation %0.32g\n",mass_approx); CHKERRQ(ierr);
    if(is_atom_test) {
      if( abs(mass_approx - volume) > 1.0e-2) {
        SETERRQ(PETSC_COMM_SELF,MOFEM_ATOM_TEST_INVALID,"atom test diverged");
      }
    }
    ierr = VecDestroy(&mass_vec_approx); CHKERRQ(ierr);
 
    PostProcVolumeOnRefinedMesh post_proc(m_field);
    ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
    ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
    ierr = post_proc.addFieldValuesPostProc("RHO"); CHKERRQ(ierr);
    ierr = post_proc.addFieldValuesPostProc("DISTANCE_FROM_SURFACE"); CHKERRQ(ierr);
    ierr = post_proc.addFieldValuesGradientPostProc("RHO"); CHKERRQ(ierr);
    ierr = post_proc.addFieldValuesGradientPostProc("DISTANCE_FROM_SURFACE"); CHKERRQ(ierr);
    ierr = DMoFEMLoopFiniteElements(dm_rho,"DENSITY_MAP_ELEMENT",&post_proc); CHKERRQ(ierr);
    ierr = post_proc.writeFile("out.h5m"); CHKERRQ(ierr);
 
    ierr = VecDestroy(&F); CHKERRQ(ierr);
    ierr = VecDestroy(&D); CHKERRQ(ierr);
    ierr = MatDestroy(&A); CHKERRQ(ierr);
    ierr = DMDestroy(&dm_rho);
 
    // Delete element which will be no longer used
    ierr = m_field.delete_finite_element("DENSITY_MAP_ELEMENT",2); CHKERRQ(ierr);
    if (m_field.get_comm_rank() == 0)
      CHKERR moab.write_file("analysis_mesh.h5m");
 
 
  } CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
 
  return (0);
}