TriWildMesh.h

#pragma once
 
#include <wmtk/TriOptimizerMesh.h>
 
#include <wmtk/SurfaceTagAttributes.h>
#include <cstdlib>
 
#include <wmtk/utils/PartitionMesh.h>
#include <wmtk/utils/VectorUtils.h>
#include <polysolve/nonlinear/Problem.hpp>
#include <wmtk/AttributeCollection.hpp>
#include <wmtk/Types.hpp>
#include <wmtk/envelope/Envelope.hpp>
#include <wmtk/optimization/SmoothVertex.hpp>
#include <wmtk/optimization/solver.hpp>
#include <wmtk/threading/enumerable_thread_specific.hpp>
 
// clang-format off
#include <wmtk/utils/DisableWarnings.hpp>
#include <fastenvelope/FastEnvelope.h>
#include <wmtk/utils/EnableWarnings.hpp>
// clang-format on
 
#include "Parameters.h"
 
#include <atomic>
#include <memory>
#include <set>
 
namespace wmtk::components::triwild {
 
using VertexAttributes = wmtk::TriOptimizerMesh::VertexAttributes;
using EdgeAttributes = wmtk::TriOptimizerMesh::EdgeAttributes;
using FaceAttributes = wmtk::TriOptimizerMesh::FaceAttributes;
 
inline constexpr size_t NO_FEATURE = std::numeric_limits<size_t>::max();
 
class TriWildMesh : public wmtk::TriOptimizerMesh
{
public:
    struct VertexExtras
    {
        size_t m_feature_id = std::numeric_limits<size_t>::max();
    };
    wmtk::AttributeCollection<VertexExtras> m_vertex_extra;
 
    Parameters& m_tri_params;
 
    std::vector<Vector2d> m_V_envelope;
    std::vector<Vector2i> m_E_envelope;
    double m_feature_eps = -1;
 
    std::vector<Vector2d> m_feature_points;
    std::atomic<size_t> m_feature_rejects = 0;
 
    bool smoothing_position_is_allowed(const size_t vid, const Vector2d& p) const override;
    void split_after_vertex(const size_t vid) override
    {
        m_vertex_extra[vid].m_feature_id = NO_FEATURE;
    }
 
    std::pair<size_t, size_t> feature_retention(double* worst_ratio = nullptr) const;
 
    bool collapse_breaks_feature(const size_t v1_id, const size_t v2_id) const;
 
 
    TriWildMesh(Parameters& _m_params, double envelope_eps, int _num_threads = 0)
        : wmtk::TriOptimizerMesh(_m_params)
        , m_tri_params(_m_params)
        , m_feature_eps(envelope_eps)
    {
        m_vertex_attr_group.add(&m_vertex_extra);
        m_envelope_eps = envelope_eps;
        NUM_THREADS = _num_threads;
 
        optimization::deactivate_opt_logger();
 
        double& wa = m_params.w_amips;
        double& we = m_params.w_envelope;
        we = 1 - wa;
        logger().info("w_envelope = {}", we);
    }
 
    ~TriWildMesh() {}
 
public:
    void init_mesh(
        const MatrixXd& V,
        const std::vector<Vector2r>& V_rational,
        const MatrixXi& F,
        const MatrixXi& E,
        const std::vector<std::string>& tag_names,
        const MatrixXd& V_env,
        const MatrixXi& E_env);
 
    void init_surfaces_and_boundaries();
 
    void init_envelope(const MatrixXd& V, const MatrixXi& F);
 
    size_t refine_sizing_around_worst(double max_energy) override;
 
    void write_msh_groups(std::string file, const bool write_envelope = true);
 
    void write_vtu(const std::string& path) const;
 
public:
    void compute_winding_numbers(const std::vector<MatrixXd>& Vs, const std::vector<MatrixXi>& Es);
 
    void filter_with_input_winding_number();
    void filter_with_flood_fill();
 
    int flood_fill();
 
protected:
    void write_smoothing_debug_output(const std::string& path) const override { write_vtu(path); }
 
    void collapse_pass_begin() override { m_feature_rejects = 0; }
    void collapse_pass_end(size_t) override
    {
        if (const size_t n = m_feature_rejects.load(); n > 0) {
            logger().info(
                "[feature] {} collapses refused to keep a polyline endpoint or junction "
                "within {:.6} of its input position",
                n,
                m_envelope_eps);
        }
    }
    bool collapse_before_vertex(size_t v1, size_t v2) override
    {
        if (!collapse_breaks_feature(v1, v2)) return true;
        m_feature_rejects.fetch_add(1, std::memory_order_relaxed);
        return false;
    }
    void collapse_after_vertex(size_t v1, size_t v2) override
    {
        if (m_vertex_extra[v2].m_feature_id == NO_FEATURE) {
            m_vertex_extra[v2].m_feature_id = m_vertex_extra[v1].m_feature_id;
        }
    }
};
 
 
} // namespace wmtk::components::triwild