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;
 
static MoFEMErrorCode extract_mesh_edges(MoFEM::Interface &m_field,
                                       Range &output_edges,
                                       Range &output_vertices);
 
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 extract_mesh_edges(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 = AdolCOps::BARRIER_AND_QUALITY;
    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();
}
 
static MoFEMErrorCode extract_mesh_edges(MoFEM::Interface &m_field,
                                       Range &output_edges,
                                       Range &output_vertices) {
  MoFEMFunctionBegin;
 
  Range fineFeatureTris;
  FTENSOR_INDEX(3, i);
  Range vols;
  CHKERR m_field.get_moab().get_entities_by_dimension(0, 3, vols);
  if (vols.empty()) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "no 3d entities found in the mesh");
  }
 
  Skinner skinner(&m_field.get_moab());
  Range skin_tris;
  CHKERR skinner.find_skin(0, vols, false, skin_tris);
  if (skin_tris.empty()) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "no skin triangles found in the mesh");
  }
 
  auto tri_normal = [&](const EntityHandle tri,
                        FTensor::Tensor1<double, 3> &normal) {
    MoFEMFunctionBeginHot;
    CHKERR m_field.getInterface<Tools>()->getTriNormal(tri, &normal(0));
    normal.normalize();
    MoFEMFunctionReturnHot(0);
  };
 
  const double fine_planar_cos = std::cos(1e-6); //FIXME: parameter?
 
  Range skin_edges;
  CHKERR m_field.get_moab().get_adjacencies(skin_tris, 1, true, skin_edges,
                                            moab::Interface::UNION);
 
  for (auto edge : skin_edges) {
    Range adj_tris;
    CHKERR m_field.get_moab().get_adjacencies(&edge, 1, 2, false, adj_tris,
                                              moab::Interface::UNION);
    adj_tris = intersect(adj_tris, skin_tris);
    if (adj_tris.size() != 2) {
      continue;
    }
 
    FTensor::Tensor1<double, 3> n0;
    FTensor::Tensor1<double, 3> n1;
    CHKERR tri_normal(adj_tris[0], n0);
    CHKERR tri_normal(adj_tris[1], n1);
    const double dot = n0(i) * n1(i);
    if (std::abs(dot) <= fine_planar_cos) {
      fineFeatureTris.insert(adj_tris[0]);
      fineFeatureTris.insert(adj_tris[1]);
      output_edges.insert(edge);
    }
  }
 
  std::map<EntityHandle, std::vector<EntityHandle>> vert_edges;
  for (auto edge : output_edges) {
    const EntityHandle *conn = nullptr;
    int nconn = 0;
    CHKERR m_field.get_moab().get_connectivity(edge, conn, nconn, false);
    if (nconn != 2)
      continue;
    vert_edges[conn[0]].push_back(edge);
    vert_edges[conn[1]].push_back(edge);
  }
 
  Range straight_edges;
  for (auto edge : output_edges) {
    const EntityHandle *conn = nullptr;
    int nconn = 0;
    CHKERR m_field.get_moab().get_connectivity(edge, conn, nconn, false);
    if (nconn != 2)
      continue;
 
    bool keep_edge = true;
    for (int vv = 0; vv < 2; ++vv) {
      const auto v = conn[vv];
      const auto it = vert_edges.find(v);
      if (it == vert_edges.end()) {
        keep_edge = false;
        break;
      }
      const auto &inc_edges = it->second;
      if (inc_edges.size() > 2) {
        keep_edge = false;
        break;
      }
      if (inc_edges.size() == 2) {
        const auto other_edge =
            (inc_edges[0] == edge) ? inc_edges[1] : inc_edges[0];
        const EntityHandle *conn_other = nullptr;
        int nconn_other = 0;
        CHKERR m_field.get_moab().get_connectivity(other_edge, conn_other,
                                                   nconn_other, false);
        if (nconn_other != 2) {
          keep_edge = false;
          break;
        }
        const EntityHandle v_other =
            (conn_other[0] == v) ? conn_other[1] : conn_other[0];
 
        double pv[3], pa[3], pb[3];
        CHKERR m_field.get_moab().get_coords(&v, 1, pv);
        EntityHandle va = conn[1 - vv];
        CHKERR m_field.get_moab().get_coords(&va, 1, pa);
        CHKERR m_field.get_moab().get_coords(&v_other, 1, pb);
 
        FTensor::Tensor1<double, 3> a;
        FTensor::Tensor1<double, 3> b;
        FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_pv(&pv[0], &pv[1],
                                                                &pv[2]);
        FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_pa(&pa[0], &pa[1],
                                                                &pa[2]);
        FTensor::Tensor1<FTensor::PackPtr<double *, 3>, 3> t_pb(&pb[0], &pb[1],
                                                                &pb[2]);
        a(i) = t_pa(i) - t_pv(i);
        b(i) = t_pb(i) - t_pv(i);
        const double na = a.l2();
        const double nb = b.l2();
        if (na == 0 || nb == 0) {
          keep_edge = false;
          break;
        }
        const double dot = (a(i) * b(i)) / (na * nb);
        if (std::abs(dot) < fine_planar_cos) {
          keep_edge = false;
          break;
        }
      }
    }
 
    if (keep_edge) {
      straight_edges.insert(edge);
    }
  }
 
  output_edges = straight_edges;
  output_vertices.clear();
  if (!output_edges.empty()) {
    std::map<EntityHandle, int> vert_count;
    for (auto edge : output_edges) {
      const EntityHandle *conn = nullptr;
      int nconn = 0;
      CHKERR m_field.get_moab().get_connectivity(edge, conn, nconn, false);
      if (nconn != 2)
        continue;
      vert_count[conn[0]]++;
      vert_count[conn[1]]++;
    }
    Range filtered_edges;
    for (auto edge : output_edges) {
      const EntityHandle *conn = nullptr;
      int nconn = 0;
      CHKERR m_field.get_moab().get_connectivity(edge, conn, nconn, false);
      if (nconn != 2)
        continue;
      const int d0 = vert_count[conn[0]];
      const int d1 = vert_count[conn[1]];
      if (!(d0 == 1 && d1 == 1)) {
        filtered_edges.insert(edge);
        for (int vv = 0; vv < 2; ++vv) {
          const auto v = conn[vv];
          auto it = vert_count.find(v);
          if (it != vert_count.end() && it->second == 1) {
            output_vertices.insert(v);
          }
        }
      }
    }
    output_edges = filtered_edges;
  }
 
  if (!fineFeatureTris.empty())
    CHKERR CutMeshInterface::SaveData(m_field.get_moab())(
        "tri_features_test.vtk", fineFeatureTris);
  if (!output_edges.empty())
    CHKERR CutMeshInterface::SaveData(m_field.get_moab())("edges_test.vtk",
                                                          output_edges);
  if (!output_vertices.empty())
    CHKERR CutMeshInterface::SaveData(m_field.get_moab())("vertices_test.vtk",
                                                          output_vertices);
 
  MoFEMFunctionReturn(0);
}