add IrregularMesh

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();
 
 

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:

src/IrregularMesh.cpp

@@ -0,0 +1,8 @@
+//
+// Created by cflin on 4/24/23.
+//
+
+#include "IrregularMesh.h"
+
+namespace top {
+} // top

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

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_;

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_);
 

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>;