Blob.h

Go to the documentation of this file.

//
//  This file is part of the FFEA simulation package
//
//  Copyright (c) by the Theory and Development FFEA teams,
//  as they appear in the README.md file.
//
//  FFEA is free software: you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation, either version 3 of the License, or
//  (at your option) any later version.
//
//  FFEA is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License
//  along with FFEA.  If not, see <http://www.gnu.org/licenses/>.
//
//  To help us fund FFEA development, we humbly ask that you cite
//  the research papers on the package.
//

/*
 *      Blob.h
 */

#ifndef BLOB_H_INCLUDED
#define BLOB_H_INCLUDED

#include "SparsityPattern.h"
#include "ConnectivityTypes.h"

#include <stdio.h>
#include <math.h>
#include <set>
#include <omp.h>
#include <algorithm>  // std::find
#include <Eigen/Sparse>
#include "FFEA_return_codes.h"
#include "mat_vec_types.h"
#include "mat_vec_fns.h"
#include "mesh_node.h"
#include "tetra_element_linear.h"
#include "SimulationParams.h"
#include "Solver.h"
#include "SparseSubstitutionSolver.h"
#include "ConjugateGradientSolver.h"
#include "MassLumpedSolver.h"
#include "NoMassCGSolver.h"
#include "Face.h"
#include "CG_solver.h"
#include "BEM_Poisson_Boltzmann.h"
#include "LJ_matrix.h"
#include "BindingSite.h"
#include "PreComp_solver.h"
#include "dimensions.h"

#ifdef USE_DOUBLE_LESS
typedef Eigen::MatrixXf Eigen_MatrixX;
typedef Eigen::VectorXf Eigen_VectorX;
typedef Eigen::Matrix3f Eigen_Matrix3;
typedef Eigen::Vector3f Eigen_Vector3;
#else
typedef Eigen::MatrixXd Eigen_MatrixX;
typedef Eigen::VectorXd Eigen_VectorX;
typedef Eigen::Matrix3d Eigen_Matrix3;
typedef Eigen::Vector3d Eigen_Vector3;
#endif


struct Blob_conf {

   // initial placement, orientation and velocity:
   int set_centroid;
   scalar centroid[3];
   int set_velocity; 
   scalar velocity[3];
   int set_rotation; 
   int rotation_type;
   scalar rotation[9];

   // kinetics:
   string states, rates;
   vector<string> maps; 
   vector<int> maps_conf_index_to, maps_conf_index_from; 
   
};


/*
 * The "Blob" class
 */
class Blob {
public:

    Blob();

    ~Blob();

    int config(const int blob_index, const int conformation_index, string node_filename, 
             string topology_filename, string surface_filename, string material_params_filename,
             string stokes_filename, string ssint_filename, string pin_filename, 
             string binding_filename, string beads_filename, scalar scale, scalar calc_compress,
             scalar compress, int linear_solver, int blob_state, SimulationParams *params,
             PreComp_params *pc_params, SSINT_matrix *ssint_matrix,
             BindingSite_matrix *binding_matrix, RngStream rng[]);
    int init();

    int update_internal_forces();

    int check_inversion();

    int update_positions();

    int reset_solver();

    void translate_linear(vector3 *vec);

    void rotate(float r11, float r12, float r13, float r21, float r22, float r23, float r31, float r32, float r33, bool beads=false);

    void rotate(float xang, float yang, float zang, bool beads=false);

    vector3 position(scalar x, scalar y, scalar z);

    void position_beads(scalar x, scalar y, scalar z);

    int forget_beads();

     void add_steric_nodes();

    void move(scalar dx, scalar dy, scalar dz);

    void get_CoM(vector3 *com);

    void get_centroid(vector3 *com);
    void calc_and_store_centroid(vector3 &com);
    vector3 calc_centroid();

    void set_pos_0();
    void kinetically_set_faces(bool state);
    int create_viewer_node_file(const char *node_filename, scalar scale);

    void write_nodes_to_file(FILE *trajectory_out);

    void pre_print(); 
    void write_pre_print_to_file(FILE *trajectory_out); 
    int toBePrinted_conf[2]; 
    int toBePrinted_state[2]; 

    int read_nodes_from_file(FILE *trajectory_out);

    void make_measurements();

    void write_measurements_to_file(FILE *fout);

    int calculate_deformation();

    scalar calc_volume();

    void make_stress_measurements(FILE *stress_out, int blob_number);

    /* DEPRECATED
     *   Will be removed.
    void calculate_vdw_bb_interaction_with_another_blob(FILE *vdw_measurement_out, int other_blob_index);
     */

    void calc_centroids_and_normals_of_all_faces();

    void calc_all_centroids();

    /*
     *
     */
    int get_num_faces();

    Face *get_face(int i);

    Face *absolutely_get_face(int i);

    tetra_element_linear *get_element(int i);

    void get_bead_position(int i, arr3 &v);

    int *get_bead_type_ptr();

    int get_bead_type(int i);

    vector<int> get_bead_assignment(int i);


    scalar get_ssint_area();

    //      /*
    //       *
    //       */
    //      int get_num_surface_elements();


    int solve_poisson(scalar *phi_gamma_IN, scalar *J_Gamma_OUT);

    int apply_ctforces();


    void zero_force();

    void set_forces_to_zero();

    // vector3 get_node(int index);
    // std::array<scalar,3> get_node(int index);
    void get_node(int index, arr3 &v);
    
    void get_node_0(int index, arr3 &v);

    void copy_node_positions(vector3 *nodes);

    vector3 ** get_actual_node_positions();

    void set_node_positions(vector3 *node_pos);

    void add_force_to_node(vector3 f, int index);

    // void zero_vdw_bb_measurement_data(); // DEPRECATED

    // void zero_vdw_xz_measurement_data(); // DEPRECATED

    void velocity_all(scalar vel_x, scalar vel_y, scalar vel_z);

    void build_poisson_matrices();

    //      int elements_are_connected(int e1, int e2);

    int build_linear_node_viscosity_matrix(Eigen::SparseMatrix<scalar> *K);

    int build_linear_node_rp_diffusion_matrix(Eigen_MatrixX *D);

    int build_linear_node_elasticity_matrix(Eigen::SparseMatrix<scalar> *A);

    int build_linear_node_mass_matrix(Eigen::SparseMatrix<scalar> *M);

    scalar get_mass();

    void enforce_box_boundaries(vector3 *box_dim);

    /* DEPRECTATED
     *   Will be removed.
     * Set the interaction flag for all faces of this Blob back to false (i.e. not interacting)
    void reset_all_faces();
     */

    void linearise_elements();

    void linearise_force();



    void compress_blob(scalar compress);

    int get_num_nodes();

    int get_num_elements();

    int get_motion_state();

    scalar get_scale();
    
    scalar get_RandU01(); 

    int get_num_linear_nodes();

    int get_num_beads();
    bool is_using_beads();

    scalar get_rmsd();

    int get_linear_solver();

    // std::array<scalar,3> get_CoG();
    // arr3 get_CoG();
    void get_stored_centroid(arr3 &cog); 

    int get_conformation_index();
    int get_previous_conformation_index();
    void set_previous_conformation_index(int index);
    int get_state_index();
    void set_state_index(int index);
    int get_previous_state_index();
    void set_previous_state_index(int index);
    BindingSite* get_binding_site(int index);

    scalar calculate_strain_energy();

    void get_min_max(vector3 *blob_min, vector3 *blob_max);

    /* Blob, conformation and state indices */
    int blob_index;
    int conformation_index, previous_conformation_index;
    int state_index, previous_state_index;

    int num_binding_sites;
    BindingSite *binding_site;

    /*
     *
     */
    //void kinetic_bind(int site_index);
    //void kinetic_unbind(int site_index);

    void pin_binding_site(set<int> node_indices);

    void unpin_binding_site(set<int> node_indices);

    void print_node_positions();
    void print_bead_positions();
    bool there_is_mass();
    void set_springs_on_blob(bool state);
    bool there_are_springs();
    bool there_are_beads();
    bool there_is_ssint();

    scalar get_kinetic_energy();
    scalar get_strain_energy();

    int pbc_count[3];

private:

    int num_nodes;

    int num_elements;

    int num_surface_elements;

    int num_interior_elements;

    int num_surface_faces;

    int num_surface_nodes;

    int num_interior_nodes;

    int num_pinned_nodes;

    int num_beads;

    int num_l_ctf;
    int num_r_ctf;
    int num_sltotal_ctf;
    int num_slsets_ctf; 

    int *ctf_slsurf_ndx;

    int blob_state;

    scalar mass;

    scalar ssint_bb_energy;

    mesh_node *node;

    vector3 **node_position;

    tetra_element_linear *elem;

    Face *surface;

    int *pinned_nodes_list;

    set<int> bsite_pinned_nodes_list;

    scalar *bead_position;

    vector <vector<int>> bead_assignment;

    int *bead_type;

    int *ctf_l_nodes;
    int *ctf_r_nodes;
    int *ctf_sl_faces;
    int *ctf_sl_surfsize;

    scalar *ctf_l_forces;
    scalar *ctf_r_forces;
    scalar *ctf_r_axis;
    char *ctf_r_type;
    scalar *ctf_sl_forces;

    string s_node_filename, s_topology_filename, s_surface_filename, 
          s_material_params_filename, s_stokes_filename, s_ssint_filename, 
          s_pin_filename, s_binding_filename, s_beads_filename;


    scalar scale;

    scalar calc_compress, compress; 

    SimulationParams *params;
    PreComp_params *pc_params;

    BindingSite_matrix *binding_matrix;

    SSINT_matrix *ssint_matrix;

    Solver *solver;

    int linear_solver;

    bool mass_in_blob;

    bool springs_on_blob;

    bool ssint_on_blob;

    bool beads_on_blob;

    vector3 *force;

    RngStream *rng;


    scalar kineticenergy, strainenergy;

    vector3 L;


    vector3 CoM, CoG, CoM_0, CoG_0;
    scalar rmsd;

    CG_solver *poisson_solver;
    SparseMatrixFixedPattern *poisson_surface_matrix;
    SparseMatrixFixedPattern *poisson_interior_matrix;
    scalar *phi_Omega;
    scalar *phi_Gamma;
    scalar *q;
    scalar *nodal_q;
    scalar *poisson_rhs;

    int *num_contributing_faces;

    connectivity_entry *element_connectivity_table;

    /*
     *
     */
    SparseMatrixFixedPattern *M;


    /*
     */
    scalar *toBePrinted_nodes; 

    int load_nodes(const char *node_filename, scalar scale);

    int load_topology(const char *topology_filename);

    int load_surface(const char *surface_filename, SimulationParams* params);


    int load_surface_no_topology(const char *surface_filename, SimulationParams *params);

    int load_material_params(const char *material_params_filename);

    int load_stokes_params(const char *stokes_filename, scalar scale);


    int load_ssint(const char *ssint_filename, int num_ssint_face_types, string ssint_method);

    int load_beads(const char *beads_filename, PreComp_params *pc_params, scalar scale);


    int load_ctforces(string ctforces_fname);


    int load_binding_sites(); // const char *binding_filename, int num_binding_site_types);

    int load_pinned_nodes(const char *pin_filename);

    int create_pinned_nodes(set<int> list);

    void calc_rest_state_info();

    /*
     *
     */
    int aggregate_forces_and_solve();

    /*
     *
     */
    void euler_integrate();

    /*
     *
     */
    int calculate_node_element_connectivity();

    int build_mass_matrix();
};

#endif