geometry_surface_fit.cpp

Go to the documentation of this file.

/** \file geometry_surface_fit.cpp
 *
 * \brief Fit GEOMETRY mid-edge nodes to a fine surface using closest points.
 *
 * Algorithm (summary):
 * - Build coarse skin edges and a KD-tree over fine triangles.
 * - Pre-filter fine triangles to keep feature regions (curved edges) and
 *   optionally drop outer slivers; build a feature-only KD-tree too.
 * - For each coarse skin edge, skip planar edges (adjacent normals nearly
 *   aligned). For curved edges, project the mid-edge point to the closest
 *   feature triangle (fallback to full fine mesh if needed) and clamp the move
 *   by a max displacement factor.
 */
 
#include <MoFEM.hpp>
#include <PostProc.hpp>
 
#include <cmath>
#include <cstdio>
#include <vector>
 
using namespace MoFEM;
 
static char help[] = "Fit GEOMETRY mid-edge nodes to a fine surface.\n\n";
 
int main(int argc, char *argv[]) {
  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
 
  try {
    DMType dm_name = "DMMOFEM";
    CHKERR DMRegister_MoFEM(dm_name);
    
    moab::Core coarse_core;
    moab::Interface &coarse_moab = coarse_core;
    MoFEM::Core core(coarse_moab);
    MoFEM::Interface &m_field = core;
 
    moab::Core fine_core;
    moab::Interface &fine_moab = fine_core;
 
    // Options (inputs/outputs and optimization parameters).
    char coarse_file[255] = "coarse.h5m";
    char surface_file[255] = "surface.h5m";
    char output_file[255] = "out_geometry_fit.h5m";
    char output_postproc[255] = "geometry_fit_postproc.h5m";
    double planar_angle_deg = 2.0;
    double max_disp_factor = 0.3;
    PetscBool write_postproc = PETSC_TRUE;
    PetscBool debug_surfaces = PETSC_FALSE;
 
    // Parse CLI options.
    PetscOptionsBegin(m_field.get_comm(), "", "geometry fit tool", "none");
    CHKERR PetscOptionsString("-coarse_file", "coarse 3d mesh file name", "",
                              coarse_file, coarse_file, 255, PETSC_NULLPTR);
    CHKERR PetscOptionsString("-surface_file", "fine surface mesh file name",
                              "", surface_file, surface_file, 255,
                              PETSC_NULLPTR);
    CHKERR PetscOptionsString("-output_file", "output mesh file name", "",
                              output_file, output_file, 255, PETSC_NULLPTR);
    CHKERR PetscOptionsString("-output_postproc",
                              "output postprocess mesh file name", "",
                              output_postproc, output_postproc, 255,
                              PETSC_NULLPTR);
    CHKERR PetscOptionsReal(
        "-planar_angle_deg",
        "angle threshold in degrees to skip planar edges", "",
        planar_angle_deg, &planar_angle_deg, PETSC_NULLPTR);
    CHKERR PetscOptionsReal(
        "-max_midnode_disp_factor",
        "skip adjustment if |delta| > factor * edge length", "",
        max_disp_factor, &max_disp_factor, PETSC_NULLPTR);
    CHKERR PetscOptionsBool("-write_postproc",
                            "write postprocess mesh output", "",
                            write_postproc, &write_postproc, PETSC_NULLPTR);
    CHKERR PetscOptionsBool("-debug",
                            "dump intermediate fine-surface feature VTKs", "",
                            debug_surfaces, &debug_surfaces, PETSC_NULLPTR);
    PetscOptionsEnd();
 
    MOFEM_LOG_CHANNEL("WORLD");
    MOFEM_LOG_TAG("WORLD", "geometry_surface_fit");
 
    // Load coarse (volume) and fine (surface) meshes.
 
    auto simple = m_field.getInterface<Simple>();
    CHKERR simple->loadFile("", std::string(coarse_file));
    CHKERR fine_moab.load_file(surface_file);
 
    // Coarse volume mesh must exist to define skin edges.
    Range coarse_vols;
    CHKERR coarse_moab.get_entities_by_dimension(0, 3, coarse_vols);
    if (coarse_vols.empty()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "no 3d entities found in coarse mesh");
    }
 
 
    // Build fine surface triangles (use skin if volume exists).
    Range fine_vols;
    Range fine_tris;
    CHKERR fine_moab.get_entities_by_dimension(0, 3, fine_vols);
    if (!fine_vols.empty()) {
      Skinner fine_skinner(&fine_moab);
      CHKERR fine_skinner.find_skin(0, fine_vols, false, fine_tris);
    } else {
      CHKERR fine_moab.get_entities_by_dimension(0, 2, fine_tris);
    }
    if (fine_tris.empty()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "no surface triangles found in fine mesh");
    }
 
    bool has_geometry = false;
    auto base = AINSWORTH_LEGENDRE_BASE;
    if (m_field.check_field("GEOMETRY")) {
      has_geometry = true;
      base = m_field.get_field_structure("GEOMETRY")->getApproxBase();
    }
    simple->addDomainField("GEOMETRY", H1, base, 3);
    CHKERR simple->setFieldOrder("GEOMETRY", 2);
    CHKERR simple->setUp(PETSC_FALSE);
 
    if (!has_geometry) {
      Projection10NodeCoordsOnField ent_method(m_field, "GEOMETRY");
      CHKERR m_field.loop_dofs("GEOMETRY", ent_method);
    }
 
    // Project mid-edge DOFs to fine surface using closest-point tagging.
    Projection10NodeCoordsFromFineSurfaceOnField fit_method(
        m_field, fine_moab, fine_tris, "GEOMETRY", planar_angle_deg,
        max_disp_factor, debug_surfaces);
    CHKERR m_field.loop_entities("GEOMETRY", fit_method);
 
    // Write updated mesh + optional outputs.
    CHKERR coarse_moab.write_file(output_file);
    MOFEM_LOG("WORLD", Sev::inform) << "Wrote file " << output_file;
    if (write_postproc) {
      auto pip = m_field.getInterface<PipelineManager>();
 
      auto post_proc_mesh = boost::make_shared<moab::Core>();
      auto post_proc_begin =
          boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(m_field,
                                                                post_proc_mesh);
      auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
          m_field, post_proc_mesh);
 
      using DomainEle = PipelineManager::VolEle;
      using PostProcEleDomain = PostProcBrokenMeshInMoabBaseCont<DomainEle>;
      auto post_proc_fe =
          boost::make_shared<PostProcEleDomain>(m_field, post_proc_mesh);
      CHKERR AddHOOps<3, 3, 3>::add(post_proc_fe->getOpPtrVector(), {H1},
                                    "GEOMETRY");
      auto x_ptr = boost::make_shared<MatrixDouble>();
      post_proc_fe->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<3>("GEOMETRY", x_ptr));
      using OpPPMap = OpPostProcMapInMoab<3, 3>;
      post_proc_fe->getOpPtrVector().push_back(new OpPPMap(
          post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(), {},
          {{"GEOMETRY", x_ptr}}, {}, {}));
 
      pip->getDomainRhsFE() = post_proc_fe;
 
      CHKERR DMoFEMPreProcessFiniteElements(simple->getDM(),
                                            post_proc_begin->getFEMethod());
      CHKERR pip->loopFiniteElements();
      CHKERR DMoFEMPostProcessFiniteElements(simple->getDM(),
                                             post_proc_end->getFEMethod());
      CHKERR post_proc_end->writeFile(output_postproc);
      MOFEM_LOG("WORLD", Sev::inform)
          << "Wrote postprocess mesh " << output_postproc;
    }
  }
  CATCH_ERRORS;
 
  MoFEM::Core::Finalize();
  return 0;
}