Mesh.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <math.h>
#include <cmath>
#include <vector>
#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <stdexcept>
#include <array>
#include <unordered_map>
#include <omp.h>
#include <algorithm>
// model setup
//#include "Edge.hpp"
#include "mesh/Face.hpp"
#include "mesh/Vertex.hpp"
#include "energy_force/Energy.hpp"
#include "energy_force/Force.hpp"
#include "mesh/Gauss_quadrature.hpp"
// matrix math
#include "linalg/Linear_algebra.hpp"
// parameters
#include "Parameters.hpp"
 
 
using namespace std;
 
class Mesh
{
public:
    struct GagReactionTarget
    {
        std::string interface1;
        std::string interface2;
        double sigma = 0.0;
        std::array<double, 5> assocAngles = {{0.0, 0.0, 0.0, 0.0, 0.0}};
    };
 
    struct GagBond
    {
        int point1 = -1;
        int point2 = -1;
        std::string interface1;
        std::string interface2;
        double expectedLength = 0.0;
        double thetaAtPoint1 = 0.0;
        double thetaAtPoint2 = 0.0;
        double expectedPhi = 0.0;
        double expectedOmega = 0.0;
    };
 
    struct GagSubunit
    {
        int pointIndex = -1;
        int moleculeId = -1;
        Matrix rotation = Matrix(3, 3, true);
        std::unordered_map<std::string, Matrix> localInterfaces;
        Matrix localReferenceVector = Matrix(3, 1);
        bool hasReferenceVector = false;
    };
 
    struct GagInteraction
    {
        int subunit1 = -1;
        int subunit2 = -1;
        std::string interface1;
        std::string interface2;
        double expectedSigma = 0.0;
        std::array<double, 5> assocAngles = {{0.0, 0.0, 0.0, 0.0, 0.0}};
    };
 
    struct GagAngle
    {
        int point1 = -1;
        int point2 = -1;
        int point3 = -1;
        double expectedAngle = 0.0;
    };
 
    struct GagTorsion
    {
        int point1 = -1;
        int point2 = -1;
        int point3 = -1;
        int point4 = -1;
        double expectedDihedral = 0.0;
        double expectedPhi = 0.0;
        double expectedOmega = 0.0;
    };
 
    std::vector<Vertex> vertices; 
    std::vector<Face> faces;      
    Param& param;                  
    // Scaffolding points, see Scaffolding_points.cpp
    Matrix centerScaffoldingSphere; 
    Matrix forceTotalOnScaffolding; 
    Matrix scaffoldingMovementVector; 
    std::vector<Matrix> forceOnScaffoldingPoints; 
    bool gagScaffoldingTopologyInitialized = false; 
    std::vector<GagBond> gagBonds; 
    std::vector<GagAngle> gagAngles; 
    std::vector<GagTorsion> gagTorsions; 
    std::vector<GagSubunit> gagSubunits; 
    std::vector<GagInteraction> gagInteractions; 
    Matrix gagInitialAlignmentRotation = Matrix(3, 3, true); 
 
    // New members... for halfedge mesh
    //std::vector<Edge> edges; ///< Vector to store all edges in the mesh
    //std::vector<Halfedge> halfedges; ///< Vector to store all halfedges in the mesh
 
 
    Mesh(Param &srcParam);
 
    Mesh(const std::vector<Vertex> &srcVertices, const std::vector<Face> &srcFaces, Param &srcParam)
    __attribute__((deprecated("Initialize with Mesh(Param &srcParam) and setup with setup_from_vertices_faces instead.")));
 
 
    //void initializeHalfedges();
 
    //Halfedge* createHalfedge(int vertexIndex, int faceIndex);
 
    // setup
 
    void setup_flat();
 
    void setup_from_vertices_faces(const std::vector<std::vector<double>>& verticesData, 
                                   const std::vector<std::vector<int>>& facesData);
 
    void set_axes_division_flat();
 
    void set_vertices_faces_flat();
 
    void set_adjacent_faces_of_vertices_sorted();
 
    bool faces_share_edge(const Face& face1, const Face& face2);
 
    bool faces_share_edge(const Face& face1, const Face& face2, std::vector<int>& commonElements);
 
    void set_adjacent_faces_of_faces();
 
    void set_adjacent_vertices_of_vertices_sorted();
 
    int find_opposite_node_index(const int &node1, const int &node2, const int &node3);
 
    void set_one_ring_vertices_sorted();
 
    void determine_boundary_vertices_faces();
 
    void determine_ghost_vertices_faces();
 
    void sort_vertices_on_faces();
 
    void set_insertion_patch(const vector<vector<int>> &insertionPatch);
 
    void set_spontaneous_curvature_for_face(const double &insertCurv, const double &spontCurv);
 
    void calculate_element_area_volume();
 
    void sum_membrane_area_and_volume(double &area, double &volume);
 
    void Compute_Energy_And_Force();
 
    void element_energy_force_regular(const std::vector<Matrix> &coordOneRingVertices,
                                      Face& face,
                                      const double spontCurv,
                                      double &meanCurv,
                                      Matrix &normVector,
                                      double &eBend,
                                      Matrix &fBend,
                                      Matrix &fArea,
                                      Matrix &fVolume);
 
    void energy_force_regularization();
 
    void manage_force_for_boundary_ghost_vertex();
 
    double get_max_force_magnitude();
 
    double calculate_mean_force();
 
    // defines in update.cpp
 
    void update_previous_coord_for_vertex();
 
    void update_reference_coord_from_previous_coord();
 
    void update_previous_force_for_vertex();
 
    void update_previous_energy_for_face();
 
    void clear_force_on_vertices_and_energy_on_faces();
 
    // implemented in scaffolding_points.cpp
 
    int findClosestRcap(double membraneCapRadius);
 
    double interpolateHeight(double r, const std::vector<double>& r_vals, const std::vector<double>& h_vals);
 
    // --- 3. Load and use in case3 ---
    void loadMembraneShapeProfile(double membraneCapRadius,
                                std::vector<double>& r_vals,
                                std::vector<double>& h_vals);
 
    bool move_vertices_based_on_scaffolding(bool fixDir = true);
 
    Matrix& find_center_of_scaffolding_sphere(bool use_default);
 
 
    double approximate_scaffolding_cap_radius(bool use_default);
 
    void set_scaffolding_vertices_correspondence();
 
    void set_scaffolding_insertion_curvature();
 
    double calculate_scaffolding_energy_force(bool doLocalSearch);
 
    bool initialize_gag_scaffolding_topology();
 
    bool orient_scaffolding_plane_to_membrane();
 
    double calculate_gag_scaffolding_internal_energy() const;
 
    double calculate_gag_scaffolding_internal_energy(const std::vector<Matrix> &scaffoldingPoints) const;
 
    double calculate_gag_scaffolding_internal_energy(const std::vector<Matrix> &scaffoldingPoints,
                                                     const std::vector<GagSubunit> &subunits) const;
 
    std::vector<Matrix> calculate_gag_scaffolding_forces_fd() const;
 
    std::vector<Matrix> calculate_gag_subunit_translation_forces_fd();
 
    std::vector<Matrix> calculate_gag_subunit_rotation_torques_fd();
 
    bool propagate_scaffolding();
 
    bool pre_relax_gag_scaffolding();
 
    // io
 
    void write_faces_csv(const std::string &outfile_name);
 
    void write_vertices_csv(const std::string &outfile_name);
 
    void write_vertices_csv_with_type(const std::string &outfile_name);
 
    void write_gag_scaffolding_state_dat(const std::string &outfile_name) const;
 
    friend std::ostream &operator<<(std::ostream &stream, const Mesh &mesh);
 
protected:
    void enumerate_gauss_quadrature_point_area_volume(const Matrix &dots,
                                                      double &area,
                                                      double &volume);
 
    Matrix get_one_ring_vertex_matrix(const Face &face);
 
    double get_squared_distance_sp_and_v(const Matrix &scaffoldingPoint, const Vertex &vertex);
};