diff --git a/main.cpp b/main.cpp
index db92772..beb98b3 100644
--- a/main.cpp
+++ b/main.cpp
@@ -15,19 +15,22 @@ void print_tensor(const Tensor3f &t3) {
}
int main() {
-
auto para = std::make_shared<top::CtrlPara>();
- para->penal=3;
- para->max_loop=4;
+ para->max_loop=200;
+ para->volfrac=0.2;
+ para->penal=3.0;
+ para->r_min=1.5;
double E = 1.0;
double Poisson_ratio = 0.3;
auto material = std::make_shared<SIM::Material>(E, Poisson_ratio);
- auto mesh = std::make_shared<top::Mesh>(10, 200,10);
+ auto mesh = std::make_shared<top::Mesh>(20, 120,20);
top::Top3d top3d(para, material, mesh);
top::Boundary bdr(mesh);
- auto l_bd=bdr.GetLeftBoundary();
- top3d.AddDBC(l_bd,{1,1,1});
- top3d.AddNBC(bdr.GetRightBottomMidPoint(),{0,0,-1});
+ top3d.AddDBC(bdr.GetCornerPoint(0),{1,1,1});
+ top3d.AddDBC(bdr.GetCornerPoint(1),{1,1,1});
+ top3d.AddDBC(bdr.GetCornerPoint(2),{1,1,1});
+ top3d.AddDBC(bdr.GetCornerPoint(3),{1,1,1});
+ top3d.AddNBC(bdr.GetTopSurfaceCenter(),{0,0,-1});
top3d.TopOptMainLoop();
diff --git a/ref/top3d.pdf b/ref/top3d.pdf
new file mode 100644
index 0000000..bc86fa9
Binary files /dev/null and b/ref/top3d.pdf differ
diff --git a/src/Boundary.h b/src/Boundary.h
index 68abc75..4261151 100644
--- a/src/Boundary.h
+++ b/src/Boundary.h
@@ -16,26 +16,61 @@ namespace top {
}
+ Eigen::MatrixXi GetTopBoundary() {
+ Eigen::VectorXi x_sequence = Eigen::VectorXi::LinSpaced(sp_mesh_->GetLx() + 1, 0, sp_mesh_->GetLx());
+ Eigen::VectorXi y_sequence = Eigen::VectorXi::LinSpaced(sp_mesh_->GetLy() + 1, 0, sp_mesh_->GetLy());
+ Eigen::VectorXi z_sequence = Eigen::Vector<int, 1>(sp_mesh_->GetLz());
+ return GetChosenCoords(x_sequence, y_sequence, z_sequence);
+ }
+
+ Eigen::MatrixXi GetBottomBoundary() {
+ Eigen::VectorXi x_sequence = Eigen::VectorXi::LinSpaced(sp_mesh_->GetLx() + 1, 0, sp_mesh_->GetLx());
+ Eigen::VectorXi y_sequence = Eigen::VectorXi::LinSpaced(sp_mesh_->GetLy() + 1, 0, sp_mesh_->GetLy());
+ Eigen::VectorXi z_sequence = Eigen::Vector<int, 1>(0);
+ return GetChosenCoords(x_sequence, y_sequence, z_sequence);
+ }
+
Eigen::MatrixXi GetLeftBoundary() {
Eigen::VectorXi x_sequence = Eigen::VectorXi::LinSpaced(sp_mesh_->GetLx() + 1, 0, sp_mesh_->GetLx());
+ Eigen::VectorXi y_sequence = Eigen::Vector<int, 1>(0);
Eigen::VectorXi z_sequence = Eigen::VectorXi::LinSpaced(sp_mesh_->GetLz() + 1, 0, sp_mesh_->GetLz());
- Eigen::MatrixXi bound_coords((sp_mesh_->GetLx() + 1) * (sp_mesh_->GetLz() + 1), 3);
- bound_coords.col(1).setZero();
- bound_coords.col(0) = (x_sequence * Eigen::RowVectorXi::Ones(z_sequence.size())).reshaped();
- bound_coords.col(2) = (Eigen::VectorXi::Ones(x_sequence.size()) * z_sequence.transpose()).reshaped();
- return bound_coords;
+ return GetChosenCoords(x_sequence, y_sequence, z_sequence);
}
Eigen::MatrixXi GetRightBoundary() {
Eigen::VectorXi x_sequence = Eigen::VectorXi::LinSpaced(sp_mesh_->GetLx() + 1, 0, sp_mesh_->GetLx());
+ Eigen::VectorXi y_sequence = Eigen::Vector<int, 1>(sp_mesh_->GetLy());
Eigen::VectorXi z_sequence = Eigen::VectorXi::LinSpaced(sp_mesh_->GetLz() + 1, 0, sp_mesh_->GetLz());
- Eigen::MatrixXi bound_coords((sp_mesh_->GetLx() + 1) * (sp_mesh_->GetLz() + 1), 3);
- bound_coords.col(1).setConstant(sp_mesh_->GetLy());
- bound_coords.col(0) = (x_sequence * Eigen::RowVectorXi::Ones(z_sequence.size())).reshaped();
- bound_coords.col(2) = (Eigen::VectorXi::Ones(x_sequence.size()) * z_sequence.transpose()).reshaped();
- return bound_coords;
+ return GetChosenCoords(x_sequence, y_sequence, z_sequence);
}
+ Eigen::MatrixXi GetTopSurfaceCenter(){
+ Eigen::MatrixXi coord(1, 3);
+ coord << sp_mesh_->GetLx() / 2, sp_mesh_->GetLy() / 2,sp_mesh_->GetLz();
+ return coord;
+ }
+ Eigen::MatrixXi GetCornerPoint(int i){
+ bool is_top=i%4;
+ i%=4;
+ Eigen::MatrixXi coord(1,3);
+ int lx=sp_mesh_->GetLx(),ly=sp_mesh_->GetLy(),lz=sp_mesh_->GetLz();
+ switch (i) {
+ case 0:coord<<0,0,0;
+ break;
+ case 1:coord<<lx,0,0;
+ break;
+ case 2:coord<<lx,ly,0;
+ break;
+ case 3:coord<<0,ly,0;
+ break;
+ default:
+ spdlog::debug("never arrive!");
+ }
+ if(is_top){
+ coord.col(2).setConstant(lz);
+ }
+ return coord;
+ }
Eigen::MatrixXi GetRightBottomMidPoint() {
int lx = sp_mesh_->GetLx();
int ly = sp_mesh_->GetLy();
@@ -50,7 +85,19 @@ namespace top {
lx / 2 + 1, ly, 0;
return bound_coords;
+ }
+
+// private:
+ Eigen::MatrixXi GetChosenCoords(const Eigen::VectorXi &x_sequence, const Eigen::VectorXi &y_sequence,
+ const Eigen::VectorXi &z_sequence) {
+ Eigen::MatrixXi chosen_coords(x_sequence.size() * y_sequence.size() * z_sequence.size(), 3);
+ chosen_coords.col(0) = x_sequence.replicate(y_sequence.size() * z_sequence.size(), 1);
+ chosen_coords.col(1) = y_sequence.transpose().replicate(x_sequence.size(), 1).reshaped().replicate(
+ z_sequence.size(), 1);
+ chosen_coords.col(2) = z_sequence.transpose().replicate(x_sequence.size() * y_sequence.size(),
+ 1).reshaped();
+ return chosen_coords;
}
private:
diff --git a/src/IrregularMesh.cpp b/src/IrregularMesh.cpp
new file mode 100644
index 0000000..c04d07d
--- /dev/null
+++ b/src/IrregularMesh.cpp
@@ -0,0 +1,8 @@
+//
+// Created by cflin on 4/24/23.
+//
+
+#include "IrregularMesh.h"
+
+namespace top {
+} // top
\ No newline at end of file
diff --git a/src/IrregularMesh.h b/src/IrregularMesh.h
new file mode 100644
index 0000000..1068c5b
--- /dev/null
+++ b/src/IrregularMesh.h
@@ -0,0 +1,137 @@
+//
+// Created by cflin on 4/24/23.
+//
+
+#ifndef TOP3D_IRREGULARMESH_H
+#define TOP3D_IRREGULARMESH_H
+
+#include "Mesh.h"
+
+namespace top {
+ struct ModelMesh {
+ Eigen::MatrixX3d points;
+ Eigen::MatrixXi surfaces;
+ };
+
+ class IrregularMesh : public Mesh {
+ public:
+ IrregularMesh(const ModelMesh &arbitrary_mesh) : Mesh() {
+ // get num_pixels
+ min_point_box_ = arbitrary_mesh.points.colwise().minCoeff();
+ Eigen::Vector3d box_max_point = arbitrary_mesh.points.colwise().maxCoeff();
+ Eigen::Vector3d box_len = box_max_point - min_point_box_;
+ const double percentage_of_min_len = 0.20;
+ len_pixel_ = box_len.minCoeff() * percentage_of_min_len;
+ Eigen::Vector3d d_num_pixels = box_len / len_pixel_;
+ Eigen::Vector3i num_pixels(std::ceil(d_num_pixels(0)), std::ceil(d_num_pixels(1)),
+ std::ceil(d_num_pixels(2)));
+
+ lx_ = num_pixels(0);
+ ly_ = num_pixels(1);
+ lz_ = num_pixels(2);
+ num_node_ = (lx_ + 1) * (ly_ + 1) * (lz_ + 1);
+ num_ele_ = lx_ * ly_ * lz_;
+
+ ten_ele_coord2ele_id = Tensor3i(lx_, ly_, lz_);
+ ten_node_coord2node_id.setConstant(-1);
+ ten_node_coord2node_id = Tensor3i(lx_ + 1, ly_ + 1, lz_ + 1);
+ ten_node_coord2node_id.setConstant(-1);
+
+ auto EvaluateSDF = [](const Eigen::Vector3d point) {
+ Eigen::Vector3d center(1, 1, 1);
+ double r = 5;
+ return r - (point - center).norm();
+ };
+ auto GetWorldCoordByElementIdx = [&](const Eigen::Vector3i &ele_idx) {
+ Eigen::Vector3d local_coord = ele_idx.cast<double>();
+ Eigen::Vector3d world_coord = ((local_coord.array() + 0.5) * len_pixel_).matrix() + min_point_box_;
+ return world_coord;
+ };
+ auto Is_in_model = [&](const Eigen::Vector3i &ele_idx) {
+ return EvaluateSDF(GetWorldCoordByElementIdx(ele_idx)) >= 0;
+ };
+
+ static const Eigen::MatrixXi delta_coord = (Eigen::MatrixXi(8, 3) <<
+ 0, 0, 0,
+ 1, 0, 0,
+ 1, 1, 0,
+ 0, 1, 0,
+ 0, 0, 1,
+ 1, 0, 1,
+ 1, 1, 1,
+ 0, 1, 1
+ ).finished();
+ // get num_pixel_ && fill pixel id
+ int cnt_pixel = 0;
+ for (int k = 0; k < lz_; ++k) {
+ for (int j = 0; j < ly_; ++j) {
+ for (int i = 0; i < lx_; ++i) {
+ if (Is_in_model({i, j, k})) {
+ ten_ele_coord2ele_id(i, j, k) = cnt_pixel++;
+ for (int di = 0; di < delta_coord.rows(); ++di) {
+ Eigen::Vector3i cur_delta_coord = delta_coord.row(di);
+ ten_node_coord2node_id(i + cur_delta_coord(0), j + cur_delta_coord(1),
+ k + cur_delta_coord(2)) = 1;
+ }
+ }
+ }
+ }
+ }
+ num_pixel_ = cnt_pixel;
+ // get_num_node_pixel && fill node_id
+ int cnt_node_pixel = 0;
+ for (int k = 0; k < lz_ + 1; ++k) {
+ for (int j = 0; j < ly_ + 1; ++j) {
+ for (int i = 0; i < lx_ + 1; ++i) {
+ if (ten_node_coord2node_id(i, j, k) == 1) {
+ ten_node_coord2node_id(i, j, k) = cnt_node_pixel++;
+ }
+ }
+ }
+ }
+ num_node_pixel_ = cnt_node_pixel;
+ // fill mat_ele_id2dofs_
+ mat_ele_id2dofs_.resize(num_pixel_, NUM_NODES_EACH_ELE * 3);
+ for (int k = 0; k < lz_; ++k) {
+ for (int j = 0; j < ly_; ++j) {
+ for (int i = 0; i < lx_; ++i) {
+ int cur_ele_id = ten_ele_coord2ele_id(i, j, k);
+ if (cur_ele_id == -1) {
+ continue;
+ }
+ Eigen::MatrixXi world_node_coords = delta_coord.rowwise() + Eigen::RowVector3i(i, j, k);
+ assert(world_node_coords.rows() == 8 && world_node_coords.cols() == 3);
+ Eigen::Vector<int, 8> node_ids = ten_node_coord2node_id(world_node_coords);
+ for (int nodi = 0; nodi < 8; ++nodi) {
+ mat_ele_id2dofs_(cur_ele_id, {3 * nodi, 3 * nodi + 1, 3 * nodi + 2}) = Eigen::Vector3i(
+ 3 * node_ids(nodi), 3 * node_ids(nodi) + 1, 3 * node_ids(nodi) + 2);
+ }
+ }
+ }
+ }
+
+ }
+
+ int GetNumDofs() const override {
+ return num_node_pixel_ * 3;
+ }
+
+ int GetNumEles() const override {
+ return num_pixel_;
+ }
+
+ int GetNumNodes() const override {
+ return num_node_pixel_;
+ }
+
+ private:
+ double len_pixel_;
+ double num_pixel_;
+ double num_node_pixel_;
+ Eigen::Vector3d min_point_box_;
+
+ };
+
+} // top
+
+#endif //TOP3D_IRREGULARMESH_H
diff --git a/src/Mesh.h b/src/Mesh.h
index b20498c..fa8c012 100644
--- a/src/Mesh.h
+++ b/src/Mesh.h
@@ -10,6 +10,7 @@
namespace top {
class Mesh {
public:
+ Mesh()=default;
Mesh(int len_x, int len_y, int len_z) : lx_(len_x), ly_(len_y), lz_(len_z),
num_node_((lx_ + 1) * (ly_ + 1) * (lz_ + 1)),
num_ele_(lx_ * ly_ * lz_) {
@@ -64,13 +65,13 @@ namespace top {
return lz_;
}
- int GetNumDofs() const{
+ virtual int GetNumDofs() const{
return num_node_*3;
}
- int GetNumEles()const{
+ virtual int GetNumEles()const{
return num_ele_;
}
- int GetNumNodes()const{
+ virtual int GetNumNodes()const{
return num_node_;
}
Eigen::MatrixXi GetEleId2DofsMap()const{
@@ -90,7 +91,7 @@ namespace top {
return mat_ele_id2dofs_(ele_id,all);
}
- private:
+ protected:
static const int NUM_NODES_EACH_ELE = 8;
int lx_;
int ly_;
diff --git a/src/Top3d.h b/src/Top3d.h
index 83c540c..73304a6 100644
--- a/src/Top3d.h
+++ b/src/Top3d.h
@@ -26,9 +26,9 @@ namespace top {
double volfrac = 0.5;
double penal = 1.0;
int max_loop = 100;
- int tol_x = 0.01;
double r_min = 2.0;
+ double tol_x = 0.01;
double E_factor = 1e-9;
};
@@ -39,7 +39,8 @@ namespace top {
), sp_material_(sp_material), sp_mesh_(sp_mesh) {
F_ = SpMat(sp_mesh_->GetNumDofs(), 1);
K_ = SpMat(sp_mesh_->GetNumDofs(), sp_mesh_->GetNumDofs());
-
+ spdlog::info("sizeof K: {}",K_.rows());
+ spdlog::info("start to precompute...");
Precompute();
}
@@ -50,6 +51,7 @@ namespace top {
if (directions[dir])
dofs_to_fixed_.insert(node_id_to_fix[i] * 3 + dir);
}
+ auto ttt=dofs_to_fixed_;
}
void AddNBC(const Eigen::MatrixXi &NBC_coords, const Eigen::Vector3d &forces) {
@@ -86,7 +88,6 @@ namespace top {
#ifdef USE_SUITESPARSE
Eigen::CholmodSupernodalLLT<Eigen::SparseMatrix<double>> solver;
-// Eigen::CholmodSimplicialLDLT<Eigen::SparseMatrix<double>> solver;
solver.compute(K_);
Eigen::VectorXd U = solver.solve(F_);
diff --git a/src/Util.h b/src/Util.h
index 3bce920..8c8796f 100644
--- a/src/Util.h
+++ b/src/Util.h
@@ -17,7 +17,7 @@
#include <vtkFloatArray.h>
#include <vtkStructuredGridWriter.h>
#include <Eigen/Dense>
-
+#include <spdlog/spdlog.h>
#include <iostream>
namespace top {
using Tensor3d=TensorWrapper<double,3>;