|
@ -8,6 +8,8 @@ |
|
|
#include <vector> |
|
|
#include <vector> |
|
|
#include <set> |
|
|
#include <set> |
|
|
#include <memory> |
|
|
#include <memory> |
|
|
|
|
|
#include <fstream> |
|
|
|
|
|
#include <iomanip> |
|
|
#include <Eigen/Eigen> |
|
|
#include <Eigen/Eigen> |
|
|
#include <igl/readOBJ.h> |
|
|
#include <igl/readOBJ.h> |
|
|
#include <spdlog/spdlog.h> |
|
|
#include <spdlog/spdlog.h> |
|
@ -25,7 +27,7 @@ |
|
|
#include "Config.hpp" |
|
|
#include "Config.hpp" |
|
|
#include <igl/read_triangle_mesh.h> |
|
|
#include <igl/read_triangle_mesh.h> |
|
|
|
|
|
|
|
|
namespace ipc::rigid{ |
|
|
namespace ipc::rigid { |
|
|
class UIStaticSimState; |
|
|
class UIStaticSimState; |
|
|
} |
|
|
} |
|
|
namespace ssim { |
|
|
namespace ssim { |
|
@ -42,149 +44,178 @@ namespace ssim { |
|
|
}; |
|
|
}; |
|
|
}; |
|
|
}; |
|
|
|
|
|
|
|
|
struct MaterialProperty{ |
|
|
struct MaterialProperty { |
|
|
float Youngs_Modulus=1.0f; |
|
|
float Youngs_Modulus = 1.0f; |
|
|
float Poisson_ratio=0.3f; |
|
|
float Poisson_ratio = 0.3f; |
|
|
float density=1.0f; |
|
|
float density = 1.0f; |
|
|
}; |
|
|
}; |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
class StaticSim { |
|
|
class StaticSim { |
|
|
using MeshModel = Model; |
|
|
using MeshModel = Model; |
|
|
using Mesh = Model; |
|
|
using Mesh = Model; |
|
|
|
|
|
|
|
|
friend class ipc::rigid::UIStaticSimState; |
|
|
friend class ipc::rigid::UIStaticSimState; |
|
|
|
|
|
|
|
|
public: |
|
|
public: |
|
|
|
|
|
|
|
|
StaticSim(const SimTargetOption &option, const std::string &jsonPath); |
|
|
StaticSim(const SimTargetOption &option, const std::string &jsonPath); |
|
|
|
|
|
|
|
|
~StaticSim() {} |
|
|
~StaticSim() {} |
|
|
|
|
|
|
|
|
void simulation(); |
|
|
void simulation(); |
|
|
|
|
|
|
|
|
/**
|
|
|
/**
|
|
|
* update BC's absBBox |
|
|
* update BC's absBBox |
|
|
*/ |
|
|
*/ |
|
|
void updateBC(); |
|
|
void updateBC(); |
|
|
|
|
|
|
|
|
|
|
|
/**
|
|
|
/**
|
|
|
* Given query points Q, compute their displacements and stress |
|
|
* Given query points Q, compute their displacements and stress |
|
|
* @param Q (nQ, 3), query points |
|
|
* @param Q (nQ, 3), query points |
|
|
* @param QU return value, (nQ), norm of displacements |
|
|
* @param QU return value, (nQ), norm of displacements |
|
|
* @param Qstress return value, (nQ, 6), stress size is (6) on each query point |
|
|
* @param Qstress return value, (nQ, 6), stress size is (6) on each query point |
|
|
*/ |
|
|
*/ |
|
|
void postprocess(Eigen::MatrixXd &Q, |
|
|
void postprocess(Eigen::MatrixXd &Q, |
|
|
Eigen::VectorXd &QU, |
|
|
Eigen::VectorXd &QU, |
|
|
Eigen::MatrixXd &Qstress); |
|
|
Eigen::MatrixXd &Qstress); |
|
|
|
|
|
|
|
|
Eigen::MatrixXd EvaluateTarget(SimTargetOption::Target target) { |
|
|
Eigen::MatrixXd EvaluateTarget(SimTargetOption::Target target) { |
|
|
if (!option_.is_option_set(target) || map_target_to_evaluated_[target].size()==0) { |
|
|
if (!option_.is_option_set(target)) { |
|
|
// If no cache, update option_ and map_
|
|
|
// If no cache, update option_ and map_
|
|
|
MapAppendTarget(target); |
|
|
MapAppendTarget(target); |
|
|
option_.set_option(target); |
|
|
option_.set_option(target); |
|
|
} |
|
|
|
|
|
// Return cache
|
|
|
|
|
|
return map_target_to_evaluated_[target]; |
|
|
} |
|
|
} |
|
|
// Return cache
|
|
|
|
|
|
return map_target_to_evaluated_[target]; |
|
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
// return surf tri mesh of tet mesh
|
|
|
// return surf tri mesh of tet mesh
|
|
|
Model get_mesh(); |
|
|
Model get_mesh(); |
|
|
|
|
|
|
|
|
private: |
|
|
private: |
|
|
|
|
|
|
|
|
void computeFeatures(); |
|
|
void computeFeatures(); |
|
|
|
|
|
|
|
|
void computeK(); |
|
|
void computeK(); |
|
|
|
|
|
|
|
|
void solve(); |
|
|
void solve(); |
|
|
|
|
|
|
|
|
void setBC(); |
|
|
void setBC(); |
|
|
|
|
|
|
|
|
void prepare_surf_result(); |
|
|
void prepare_surf_result(); |
|
|
|
|
|
|
|
|
void MapAppendTarget(int target); |
|
|
void MapAppendTarget(int target); |
|
|
|
|
|
|
|
|
MaterialProperty& get_material_property(){ |
|
|
Eigen::VectorXi getSurfTriForBox(const Eigen::Vector3d &minBox, |
|
|
return material_property_; |
|
|
const Eigen::Vector3d &maxBox); |
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
MaterialProperty &get_material_property() { |
|
|
|
|
|
return material_property_; |
|
|
|
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
// Return: #mesh.V.rows() x 1
|
|
|
|
|
|
Eigen::MatrixXd EvaluateUNorm() const; |
|
|
|
|
|
|
|
|
// Return: #mesh.V.rows() x 1
|
|
|
// Return: #mesh.V.rows() x 1
|
|
|
Eigen::MatrixXd EvaluateUNorm() const; |
|
|
Eigen::MatrixXd EvaluateUX() const; |
|
|
|
|
|
|
|
|
// Return: #mesh.V.rows() x 1
|
|
|
// Return: #mesh.V.rows() x 1
|
|
|
Eigen::MatrixXd EvaluateUX() const; |
|
|
Eigen::MatrixXd EvaluateUY() const; |
|
|
|
|
|
|
|
|
// Return: #mesh.V.rows() x 1
|
|
|
// Return: #mesh.V.rows() x 1
|
|
|
Eigen::MatrixXd EvaluateUY() const; |
|
|
Eigen::MatrixXd EvaluateUZ() const; |
|
|
|
|
|
|
|
|
// Return: #mesh.V.rows() x 1
|
|
|
// Return: #mesh.V.rows() x 1
|
|
|
Eigen::MatrixXd EvaluateUZ() const; |
|
|
Eigen::MatrixXd EvaluateSNorm() const; |
|
|
|
|
|
|
|
|
// Return: #mesh.V.rows() x 1
|
|
|
// Return: #mesh.V.rows() x 1
|
|
|
Eigen::MatrixXd EvaluateSNorm() const; |
|
|
Eigen::MatrixXd EvaluateSVonMises() const; |
|
|
|
|
|
|
|
|
// Return: #mesh.V.rows() x 1
|
|
|
// Return: #mesh.V.rows() x 1
|
|
|
Eigen::MatrixXd EvaluateSVonMises() const; |
|
|
Eigen::MatrixXd EvaluateSX() const; |
|
|
|
|
|
|
|
|
// Return: #mesh.V.rows() x 1
|
|
|
// Return: #mesh.V.rows() x 1
|
|
|
Eigen::MatrixXd EvaluateSX() const; |
|
|
Eigen::MatrixXd EvaluateSY() const; |
|
|
|
|
|
|
|
|
// Return: #mesh.V.rows() x 1
|
|
|
// Return: #mesh.V.rows() x 1
|
|
|
Eigen::MatrixXd EvaluateSY() const; |
|
|
Eigen::MatrixXd EvaluateSZ() const; |
|
|
|
|
|
|
|
|
// Return: #mesh.V.rows() x 1
|
|
|
// Return: 1 x 1
|
|
|
Eigen::MatrixXd EvaluateSZ() const; |
|
|
Eigen::MatrixXd EvaluateCompliance() const; |
|
|
|
|
|
|
|
|
// Return: 1 x 1
|
|
|
void writePntVTK(const std::string &path, const Eigen::MatrixXd &V) { |
|
|
Eigen::MatrixXd EvaluateCompliance() const; |
|
|
std::ofstream out(path); |
|
|
|
|
|
out << "# vtk DataFile Version 3.0\n" |
|
|
|
|
|
"Volume Mesh\n" |
|
|
|
|
|
"ASCII\n" |
|
|
|
|
|
"DATASET UNSTRUCTURED_GRID" << std::endl; |
|
|
|
|
|
out << "POINTS " << V.rows() << " float" << std::endl; |
|
|
|
|
|
for (int i = 0; i < V.rows(); ++i) { |
|
|
|
|
|
out << std::setprecision(4) << V.row(i).x() << " " << V.row(i).y() << " " << V.row(i).z() << std::endl; |
|
|
|
|
|
} |
|
|
|
|
|
out << "CELLS " << V.rows() << " " << V.rows() * (1 + 1) << std::endl; |
|
|
|
|
|
for (int i = 0; i < V.rows(); ++i) { |
|
|
|
|
|
out << "1 " << i << std::endl; |
|
|
|
|
|
} |
|
|
|
|
|
out << "CELL_TYPES " << V.rows() << std::endl; |
|
|
|
|
|
for (int i = 0; i < V.rows(); ++i) { |
|
|
|
|
|
out << 1 << std::endl; |
|
|
|
|
|
} |
|
|
|
|
|
} |
|
|
|
|
|
|
|
|
private: |
|
|
private: |
|
|
MeshModel mesh_; |
|
|
MeshModel mesh_; |
|
|
MaterialProperty material_property_; |
|
|
MaterialProperty material_property_; |
|
|
SimTargetOption option_; |
|
|
SimTargetOption option_; |
|
|
std::vector<Eigen::MatrixXd> map_target_to_evaluated_; |
|
|
std::vector<Eigen::MatrixXd> map_target_to_evaluated_; |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
private: |
|
|
private: |
|
|
std::vector<DirichletBC> DirichletBCs; |
|
|
std::vector<DirichletBC> DirichletBCs; |
|
|
std::vector<NeumannBC> NeumannBCs; |
|
|
std::vector<NeumannBC> NeumannBCs; |
|
|
Eigen::Matrix<double, 6, 6> D; // constitutive matrix
|
|
|
Eigen::Matrix<double, 6, 6> D; // constitutive matrix
|
|
|
|
|
|
|
|
|
// owned data
|
|
|
// owned data
|
|
|
int nN, nEle, nDof; |
|
|
int nN, nEle, nDof; |
|
|
Eigen::MatrixXd TV; // vertices coordinates
|
|
|
Eigen::MatrixXd TV; // vertices coordinates
|
|
|
Eigen::MatrixXd TV1; // deformed vertices coordinates
|
|
|
Eigen::MatrixXd TV1; // deformed vertices coordinates
|
|
|
Eigen::MatrixXi TT; // vertice index of each tetrahedron
|
|
|
Eigen::MatrixXi TT; // vertice index of each tetrahedron
|
|
|
Eigen::MatrixXi SF; |
|
|
Eigen::MatrixXi SF; |
|
|
|
|
|
|
|
|
int eleNodeNum; |
|
|
int eleNodeNum; |
|
|
int eleDofNum; |
|
|
int eleDofNum; |
|
|
std::vector<Eigen::VectorXi> eDof; |
|
|
std::vector<Eigen::VectorXi> eDof; |
|
|
|
|
|
|
|
|
// owned features
|
|
|
// owned features
|
|
|
Eigen::VectorXd load; // load of each dof
|
|
|
Eigen::VectorXd load; // load of each dof
|
|
|
Eigen::VectorXd U; // dofs' displacement to be computed
|
|
|
Eigen::VectorXd U; // dofs' displacement to be computed
|
|
|
|
|
|
|
|
|
Eigen::VectorXi DBC_nI; // vertex in DBC
|
|
|
Eigen::VectorXi DBC_nI; // vertex in DBC
|
|
|
Eigen::VectorXi isDBC; // 0: not in DBC, 1: in DBC
|
|
|
Eigen::VectorXi isDBC; // 0: not in DBC, 1: in DBC
|
|
|
|
|
|
|
|
|
Eigen::VectorXi SVI; // vertice indices of surface nodes
|
|
|
Eigen::VectorXi SVI; // vertice indices of surface nodes
|
|
|
Eigen::MatrixXi F_surf; // boundary vertice indices in surface triangles mesh
|
|
|
Eigen::MatrixXi F_surf; // boundary vertice indices in surface triangles mesh
|
|
|
|
|
|
std::unordered_map<int, int> vI2SVI; |
|
|
// indices for fast access
|
|
|
|
|
|
std::vector<std::set<int>> vNeighbor; // records all vertices' indices adjacent to each vertice
|
|
|
// indices for fast access
|
|
|
std::vector<std::set<std::pair<int, int>>> vFLoc; |
|
|
std::vector<std::set<int>> vNeighbor; // records all vertices' indices adjacent to each vertice
|
|
|
std::shared_ptr<LinSysSolver<Eigen::VectorXi, Eigen::VectorXd>> linSysSolver; |
|
|
std::vector<std::set<std::pair<int, int>>> vFLoc; |
|
|
|
|
|
std::shared_ptr<LinSysSolver<Eigen::VectorXi, Eigen::VectorXd>> linSysSolver; |
|
|
// resulted data to evaluate
|
|
|
|
|
|
double compliance_; |
|
|
// resulted data to evaluate
|
|
|
Eigen::MatrixXd surf_U_; |
|
|
double compliance_; |
|
|
Eigen::MatrixXd surf_stress_; |
|
|
Eigen::MatrixXd surf_U_; |
|
|
Eigen::VectorXd surf_vonstress_; |
|
|
Eigen::MatrixXd surf_stress_; |
|
|
|
|
|
Eigen::VectorXd surf_vonstress_; |
|
|
}; |
|
|
|
|
|
|
|
|
Eigen::VectorXi DBC_faceIdx_; |
|
|
|
|
|
Eigen::VectorXi DBC_vertexIdx_; |
|
|
|
|
|
Eigen::VectorXi NBC_faceIdx_; |
|
|
|
|
|
Eigen::VectorXi NBC_vertexIdx_; |
|
|
|
|
|
|
|
|
|
|
|
}; |
|
|
|
|
|
|
|
|
} // ssim
|
|
|
} // ssim
|
|
|
|
|
|
|
|
|