SimWildMeshTri.hpp

#pragma once
 
#include <wmtk/TriOptimizerMesh.h>
 
#include <wmtk/SurfaceTagAttributes.h>
#include <limits>
#include <unordered_set>
 
#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>
 
// clang-format off
#include <wmtk/utils/DisableWarnings.hpp>
#include <igl/write_triangle_mesh.h>
#include <wmtk/threading/enumerable_thread_specific.hpp>
#include <fastenvelope/FastEnvelope.h>
#include <wmtk/utils/EnableWarnings.hpp>
// clang-format on
 
#include "ConnectedComponent.hpp"
#include "Parameters.h"
#include "expression_parser/Expression.hpp"
 
 
namespace wmtk::components::simwild::tri {
 
using VertexAttributes = wmtk::TriOptimizerMesh::VertexAttributes;
using EdgeAttributes = wmtk::TriOptimizerMesh::EdgeAttributes;
using FaceAttributes = wmtk::TriOptimizerMesh::FaceAttributes;
 
class SimWildMeshTri : public wmtk::TriOptimizerMesh
{
public:
    using ExprPtr = expression_parser::ExpressionPtr;
 
    Parameters& m_sim_params;
 
    std::vector<Vector2d> m_V_envelope;
    std::vector<Vector2i> m_E_envelope;
 
    std::vector<std::tuple<ExprPtr, double>> m_sizing_field;
    std::vector<std::tuple<ExprPtr, double>> m_quality_field;
 
    bool m_collapse_check_link_condition = false; // classical link condition
    bool m_collapse_check_topology = false; // sanity check
    bool m_collapse_check_manifold = false; // manifoldness check after collapse
 
    bool smoothing_position_is_allowed(const size_t, const Vector2d&) const override
    {
        return true;
    }
 
    SimWildMeshTri(Parameters& _m_params, double envelope_eps, int _num_threads = 0)
        : wmtk::TriOptimizerMesh(_m_params)
        , m_sim_params(_m_params)
    {
        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);
    }
 
    ~SimWildMeshTri() {}
 
public:
    void init_from_image(
        const MatrixXd& V,
        const MatrixXi& T,
        const MatrixSi& T_tags,
        const std::vector<std::string>& tag_names);
 
    void init_from_image(
        const MatrixXr& V,
        const MatrixXi& T,
        const MatrixSi& T_tags,
        const std::vector<std::string>& tag_names);
 
    void init_surfaces_and_boundaries();
 
    void init_envelope(const MatrixXd& V, const MatrixXi& F);
 
    CellTag string_set_to_cell_tag(const std::set<std::string>& str_set);
 
    void set_sizing_field(const nlohmann::json& sizing_field_json);
 
    void set_quality_field(const nlohmann::json& quality_field_json);
 
    double target_quality(const size_t tid) const;
    double target_quality(const Tuple& t) const;
    double quality_rel(const size_t tid) const override;
    double quality_rel(const Tuple& t) const;
    bool check_mesh_quality(double& max_rel_quality, const bool verbose = false) const;
    std::vector<size_t> active_vertices() const override;
 
    size_t refine_sizing_around_worst(double max_metric = 0.) override;
 
    void write_msh(std::string file, const bool write_envelope = true);
 
    void write_vtu(const std::string& path) const;
    void write_vtu_with_energies(const std::string& path) const;
 
public:
    std::vector<Vector2d> get_surface_assembles(const Tuple& t) const;
    std::shared_ptr<polysolve::nonlinear::Problem> get_envelope_energy(const Tuple& t) const;
 
    std::vector<std::array<double, 6>> get_amips_assembles(const Tuple& t) const;
    std::shared_ptr<polysolve::nonlinear::Problem> get_amips_energy(const Tuple& t) const;
 
    void log_total_surface_energy();
 
 
    double triangle_area(const size_t fid) const;
 
    std::vector<ConnectedComponent> compute_connected_components(const CellTag& tag_in) const;
    std::vector<ConnectedComponent> compute_connected_components(const ExprPtr& expr) const;
 
    std::vector<ConnectedComponent> find_holes(const std::vector<CellTag>& tag_in) const;
 
    void compute_tag_boundary(const CellTag& tag, MatrixXd& V, MatrixXi& E) const;
 
    void keep_largest_connected_component(
        const std::vector<CellTag>& lcc_tags,
        const size_t n_lcc = 1);
 
    void fill_holes_topo(
        const std::vector<CellTag>& fill_holes_tags,
        double threshold = std::numeric_limits<double>::infinity());
 
    void seal_connected_components(
        const std::vector<CellTag>& tag_sets,
        const std::vector<ConnectedComponent>& components);
 
    void tight_seal_topo(
        const std::vector<std::vector<CellTag>>& tight_seal_tag_sets,
        double threshold = std::numeric_limits<double>::infinity());
 
    void resolve_overlaps(const std::vector<std::array<ExprPtr, 2>>& intersecting_tags);
 
    void replace_tags(const std::vector<CellTag>& tags_in, const std::vector<CellTag>& tags_out);
 
    void tag_priority(const std::vector<int64_t>& tags_order);
 
private:
    // When set, split_edge_after binary-searches vmid onto the zero-crossing of this function.
    // Negative = stays on v1 side, positive = stays on v2 side.
    // Set before split_edge(), cleared immediately after.
    std::function<double(const Vector2d&)> m_voronoi_split_fn = nullptr;
    size_t m_last_split_vertex = 0;
 
protected:
    std::tuple<double, double> optimization_quality_stats() override;
    double optimization_stop_metric() const override { return 1.; }
    bool optimization_stop_at_float() const override { return m_sim_params.stop_at_float; }
 
    void write_smoothing_debug_output(const std::string& path) const override { write_vtu(path); }
 
    bool collapse_quality_allowed(size_t v1, size_t fid, double q, double ring_max) const override;
    void collapse_after_vertex(size_t v1, size_t v2) override;
    bool split_adjust_position(size_t v_new, const std::vector<Tuple>& children) override;
    void split_after_vertex(size_t v_new) override { m_last_split_vertex = v_new; }
};
 
} // namespace wmtk::components::simwild::tri