mofem/tools/mesh_smoothing.cpp

/** \file mesh_smoothing.cpp
 * \brief Improve mesh quality using Volume-length quality measure with barrier
 * \example mofem/tools/mesh_smoothing.cpp
 *
 */
 
 
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <PostProc.hpp>
 
#include <MeshSmoothingOpFactory.hpp>
 
using namespace MoFEM;
 
static char help[] = "mesh smoothing\n\n";
 
int edges_block_set = 1;
int vertex_block_set = 2;
PetscBool output_vtk = PETSC_TRUE;
double tol = 0.1;
double gamma_factor = 0.8;
 
constexpr bool fix_edge_blocks = true;
 
int main(int argc, char *argv[]) {
 
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
 
  try {
 
    PetscOptionsBegin(PETSC_COMM_WORLD, "", "Mesh cut options", "none");
    CHKERR PetscOptionsInt("-edges_block_set", "edges side set", "",
                           edges_block_set, &edges_block_set, PETSC_NULLPTR);
    CHKERR PetscOptionsInt("-vertex_block_set", "vertex side set", "",
                           vertex_block_set, &vertex_block_set, PETSC_NULLPTR);
    CHKERR PetscOptionsBool("-output_vtk", "if true outout vtk file", "",
                            output_vtk, &output_vtk, PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-quality_reduction_tol", "",
                              "Tolerance of quality reduction", tol, &tol,
                              PETSC_NULLPTR);
    CHKERR PetscOptionsScalar("-gamma_factor", "",
                              "Gamma factor", gamma_factor, &gamma_factor,
                              PETSC_NULLPTR);
    PetscOptionsEnd();
 
    // Create MoAB database
    moab::Core moab_core;              // create database
    moab::Interface &moab = moab_core; // create interface to database
    // Create MoFEM database and link it to MoAB
    MoFEM::Core mofem_core(moab);           // create database
    MoFEM::Interface &m_field = mofem_core; // create interface to database
    // Register DM Manager
    CHKERR DMRegister_MoFEM("DMMOFEM"); // register MoFEM DM in PETSc
 
    // Get simple interface is simplified version enabling quick and
    // easy construction of problem.
    Simple *simple_interface = m_field.getInterface<Simple>();
    Range edges, vertices, fixed_vertex;
 
    // Build problem with simple interface
    {
      // Get options for simple interface from command line
      CHKERR simple_interface->getOptions();
      // Load mesh file to database
      CHKERR simple_interface->loadFile();
 
      // Add domain field "U" in space H1 and Legendre base, Ainsworth recipe is
      // used to construct base functions.
      CHKERR simple_interface->addDomainField("MESH_NODE_POSITIONS", H1,
                                              AINSWORTH_LEGENDRE_BASE, 3);
      // Add Lagrange multiplier field on body boundary
      CHKERR simple_interface->addBoundaryField("LAMBDA_SURFACE", H1,
                                                AINSWORTH_LEGENDRE_BASE, 1);
 
      // Set fields order domain and boundary fields.
      CHKERR simple_interface->setFieldOrder("MESH_NODE_POSITIONS",
                                             1); // to approximate function
      CHKERR simple_interface->setFieldOrder("LAMBDA_SURFACE",
                                             1); // to Lagrange multipliers
 
      simple_interface->getDomainFEName() = "SMOOTHING";
      simple_interface->getBoundaryFEName() = "SURFACE_SLIDING";
 
      // Other fields and finite elements added to database directly
        if (m_field.getInterface<MeshsetsManager>()->checkMeshset(
                edges_block_set, BLOCKSET)) {
          CHKERR m_field.getInterface<MeshsetsManager>()
              ->getEntitiesByDimension(edges_block_set, BLOCKSET, 1, edges,
                                       true);
 
        } else {
          CHKERR Smoother::extractMeshEdges(m_field, edges, fixed_vertex);
        }
        if (m_field.getInterface<MeshsetsManager>()->checkMeshset(
                vertex_block_set, BLOCKSET)) {
          CHKERR m_field.getInterface<MeshsetsManager>()
              ->getEntitiesByDimension(vertex_block_set, BLOCKSET, 0,
                                       fixed_vertex, true);
        }
 
        CHKERR m_field.get_moab().get_connectivity(edges, vertices, true);
      
 
      if(!edges.empty() && !fix_edge_blocks) {
          // Declare approximation fields
          CHKERR m_field.add_field("LAMBDA_EDGE", H1, AINSWORTH_LOBATTO_BASE, 2,
                                   MB_TAG_SPARSE, MF_ZERO);
          CHKERR m_field.add_ents_to_field_by_type(edges, MBEDGE,
                                                   "LAMBDA_EDGE");
          CHKERR m_field.getInterface<CommInterface>()
              ->synchroniseFieldEntities("LAMBDA_EDGE");
          CHKERR m_field.set_field_order(0, MBVERTEX, "LAMBDA_EDGE", 1);
 
          CHKERR m_field.add_finite_element("EDGE_SLIDING");
          CHKERR m_field.modify_finite_element_add_field_row(
              "EDGE_SLIDING", "MESH_NODE_POSITIONS");
          CHKERR m_field.modify_finite_element_add_field_col(
              "EDGE_SLIDING", "MESH_NODE_POSITIONS");
          CHKERR m_field.modify_finite_element_add_field_data(
              "EDGE_SLIDING", "MESH_NODE_POSITIONS");
          CHKERR m_field.modify_finite_element_add_field_row("EDGE_SLIDING",
                                                             "LAMBDA_EDGE");
          CHKERR m_field.modify_finite_element_add_field_col("EDGE_SLIDING",
                                                             "LAMBDA_EDGE");
          CHKERR m_field.modify_finite_element_add_field_data("EDGE_SLIDING",
                                                              "LAMBDA_EDGE");
          CHKERR m_field.add_ents_to_finite_element_by_type(edges, MBEDGE,
                                                           "EDGE_SLIDING");
          simple_interface->getOtherFiniteElements().push_back("EDGE_SLIDING");
        }
 
      CHKERR simple_interface->defineFiniteElements();
      CHKERR simple_interface->defineProblem();
      CHKERR simple_interface->buildFields();
      CHKERR simple_interface->buildFiniteElements();
      CHKERR simple_interface->buildProblem();
 
      // Remove vertices form LAMBDA_SURFACE which are on the edges
      if (!edges.empty()) {
        Range edges;
        CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
            edges_block_set, BLOCKSET, 1, edges, true);
        Range verts;
        CHKERR m_field.get_moab().get_connectivity(edges, verts, true);
        Simple *simple_interface = m_field.getInterface<Simple>();
        auto prb_mng = m_field.getInterface<ProblemsManager>();
        CHKERR prb_mng->removeDofsOnEntities(simple_interface->getProblemName(),
                                             "LAMBDA_SURFACE", verts, 0, 1);
        if (fix_edge_blocks) {
          CHKERR prb_mng->removeDofsOnEntities(
              simple_interface->getProblemName(), "MESH_NODE_POSITIONS", verts,
              0, 3);
        }
      }
    }
 
    MeshSmoothingOps::Config smoothing_cfg;
    smoothing_cfg.geometryField = "MESH_NODE_POSITIONS";
    smoothing_cfg.meshReferenceField = "NONE";
    smoothing_cfg.lambdaSurfaceField = "LAMBDA_SURFACE";
    smoothing_cfg.lambdaEdgeField = "LAMBDA_EDGE";
    smoothing_cfg.enableSurfaceSliding = true;
    smoothing_cfg.enableEdgeSliding = !fix_edge_blocks;
    smoothing_cfg.fixEdgeBlocks = fix_edge_blocks;
    smoothing_cfg.qualityType = BARRIER_AND_QUALITY;
    smoothing_cfg.alpha = 1.0;
    smoothing_cfg.gamma = 0.0;
 
    MeshSmoothingOps::Handles smoothing_handles;
    CHKERR MeshSmoothingOps::createOperators(m_field, smoothing_cfg, fixed_vertex,
                                             edges, smoothing_handles);
 
    DM dm;
    CHKERR simple_interface->getDM(&dm);
    CHKERR MeshSmoothingOps::addToSnesDM(dm, smoothing_cfg, smoothing_handles);
 
    struct Solve {
 
      MoFEMErrorCode operator()(DM dm) const {
        MoFEMFunctionBegin;
        auto solver = SmartPetscObj<SNES>();
        CHKERR MeshSmoothingOps::createSnes(dm, solver);
 
        auto hook = [&](SNES snes, Vec x, Vec F,
                                     boost::shared_ptr<CacheTuple> cache_ptr,
                                     void *ctx) -> MoFEMErrorCode {
          MoFEMFunctionBegin;
 
          using Ele = PipelineManager::ElementsAndOpsByDim<3>::BoundaryEle;
          using PostProcEleDomain = PostProcBrokenMeshInMoabBase<Ele>;
 
          auto *m_field_ptr = getInterfacePtr(dm);
          auto post_proc_fe =
              boost::make_shared<PostProcEleDomain>(*m_field_ptr);
          auto &pip = post_proc_fe->getOpPtrVector();
          auto X_ptr = boost::make_shared<MatrixDouble>();
          pip.push_back(new OpCalculateVectorFieldValues<3>(
              "MESH_NODE_POSITIONS", X_ptr));
          auto res_X_ptr = boost::make_shared<MatrixDouble>();
          pip.push_back(new OpCalculateVectorFieldValues<3>(
              "MESH_NODE_POSITIONS", res_X_ptr, SmartPetscObj<Vec>(F, true)));
          auto lambda_ptr = boost::make_shared<VectorDouble>();
          pip.push_back(
              new OpCalculateScalarFieldValues("LAMBDA_SURFACE", lambda_ptr));
          auto res_lambda_ptr = boost::make_shared<VectorDouble>();
          pip.push_back(new OpCalculateScalarFieldValues(
              "LAMBDA_SURFACE", res_lambda_ptr, SmartPetscObj<Vec>(F, true)));
 
          using OpPPMap = OpPostProcMapInMoab<3, 3>;
          pip.push_back(
 
              new OpPPMap(
 
                  post_proc_fe->getPostProcMesh(),
                  post_proc_fe->getMapGaussPts(),
 
                  {{"LAMBDA_SURFACE", lambda_ptr},
                   {"RES_LAMBDA_SURFACE", res_lambda_ptr}},
 
                  {{"MESH_NODE_POSITIONS", X_ptr},
 
                   {"RES_MESH_NODE_POSITIONS", res_X_ptr}},
 
                  {},
 
                  {}
 
                  )
 
          );
 
          CHKERR DMoFEMLoopFiniteElements(dm, "SURFACE_SLIDING", post_proc_fe);
          CHKERR post_proc_fe->writeFile("debug_smoothing.h5m");
          MoFEMFunctionReturn(0);
        };
 
        auto snes_ctx_ptr = getDMSnesCtx(dm);
        snes_ctx_ptr->getRhsHook() = hook;
 
        CHKERR MeshSmoothingOps::solveSnes(dm, solver, {"LAMBDA_SURFACE"});
 
        MoFEMFunctionReturn(0);
      }
 
      MoFEMErrorCode setCoordsFromField(DM dm) const {
        MoFEMFunctionBegin;
        MoFEM::Interface *m_field_ptr;
        CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
        for (_IT_GET_ENT_FIELD_BY_NAME_FOR_LOOP_(*m_field_ptr,
                                                 "MESH_NODE_POSITIONS", it)) {
          if (it->get()->getEntType() != MBVERTEX)
            continue;
          VectorDouble3 coords(3);
          for(int dd = 0;dd!=3;++dd)
            coords[dd] = it->get()->getEntFieldData()[dd];
          EntityHandle ent = it->get()->getEnt();
          CHKERR m_field_ptr->get_moab().set_coords(&ent, 1, &*coords.begin());
        }
        MoFEMFunctionReturn(0);
      }
 
      MoFEMErrorCode setFieldFromCoords(DM dm) const {
        MoFEMFunctionBegin;
        MoFEM::Interface *m_field_ptr;
        CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
        for (_IT_GET_ENT_FIELD_BY_NAME_FOR_LOOP_(*m_field_ptr,
                                                 "MESH_NODE_POSITIONS", it)) {
          if (it->get()->getEntType() != MBVERTEX)
            continue;
          EntityHandle ent = it->get()->getEnt();
          VectorDouble3 coords(3);
          CHKERR m_field_ptr->get_moab().get_coords(&ent, 1, &*coords.begin());
          for(int dd = 0;dd!=3;++dd) 
            it->get()->getEntFieldData()[dd] = coords[dd];
        }
        MoFEMFunctionReturn(0);
      }
 
    };
 
    Solve solve;
    CHKERR solve.setFieldFromCoords(dm);
 
    double min_quality = 0;
    CHKERR MeshSmoothingOps::evaluateMinQuality(
        dm, simple_interface->getDomainFEName(), smoothing_handles, min_quality);
    PetscPrintf(PETSC_COMM_WORLD, "Min quality = %4.3f\n", min_quality);
 
    double gamma = min_quality > 0 ? gamma_factor * min_quality
                                   : min_quality / gamma_factor;
    CHKERR MeshSmoothingOps::updateBarrierGamma(smoothing_handles, gamma);
 
    double min_quality_p, eps;
    do {
 
      min_quality_p = min_quality;
 
      CHKERR solve(dm);
 
      CHKERR solve.setCoordsFromField(dm);
      CHKERR MeshSmoothingOps::evaluateMinQuality(
          dm, simple_interface->getDomainFEName(), smoothing_handles,
          min_quality);
 
      eps = (min_quality - min_quality_p) / min_quality;
      PetscPrintf(PETSC_COMM_WORLD, "Min quality = %4.3f eps = %4.3f\n",
                  min_quality, eps);
 
      double gamma = min_quality > 0 ? gamma_factor * min_quality
                                     : min_quality / gamma_factor;
      CHKERR MeshSmoothingOps::updateBarrierGamma(smoothing_handles, gamma);
 
    } while (eps > tol);
 
    // if (m_field.getInterface<MeshsetsManager>()->checkMeshset(edges_block_set,
    //                                                           BLOCKSET)) {
    //   Range edges;
    //   CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
    //       edges_block_set, BLOCKSET, 1, edges, true);
 
    //   Range tets;
    //   CHKERR moab.get_entities_by_type(0,MBTET,tets);
    //   Skinner skin(&moab);
    //   Range skin_faces; // skin faces from 3d ents
    //   CHKERR skin.find_skin(0, tets, false, skin_faces);
 
    //   CHKERR EdgeSlidingConstrains::CalculateEdgeBase::setTags(moab, edges,
    //                                                          skin_faces);
    //   CHKERR EdgeSlidingConstrains::CalculateEdgeBase::saveEdges(
    //       moab, "out_edges.vtk", edges);
    // }
 
    if (output_vtk)
      CHKERR m_field.getInterface<BitRefManager>()->writeBitLevelByType(
          BitRefLevel().set(0), BitRefLevel().set(), MBTET, "out.vtk", "VTK",
          "");
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
}