From 23d3793112c863f6067076e2f35efbce6ef6f94e Mon Sep 17 00:00:00 2001 From: cflin Date: Thu, 6 Jul 2023 17:24:33 +0800 Subject: [PATCH] before clear --- da-sha/sha-topology-optimization-3d/Mesh.h | 2 +- .../CMakeLists.txt | 9 + .../config.json | 67 +++++++ .../main.cpp | 171 +++++++++++++++++ .../config.json | 54 +++--- .../config_case2.json | 65 +++++++ .../top-thermoelastic-regular-model/main.cpp | 72 +++++--- src/Mesh/HeatMesh.h | 7 +- src/Mesh/IrregularMesh.cpp | 6 +- src/Mesh/Mesh.cpp | 173 +++++++++++++++--- src/Mesh/Mesh.h | 28 ++- src/TensorWrapper.h | 14 ++ src/ThermoelasticTop3d.cpp | 82 +++++---- src/ThermoelasticTop3d.h | 1 - src/Top3d.cpp | 6 +- src/Util.cpp | 2 +- 16 files changed, 635 insertions(+), 124 deletions(-) create mode 100644 examples/top-thermoelastic-Lshape-condition/CMakeLists.txt create mode 100644 examples/top-thermoelastic-Lshape-condition/config.json create mode 100644 examples/top-thermoelastic-Lshape-condition/main.cpp create mode 100644 examples/top-thermoelastic-regular-model/config_case2.json diff --git a/da-sha/sha-topology-optimization-3d/Mesh.h b/da-sha/sha-topology-optimization-3d/Mesh.h index 2aaa69d..a9115a1 100644 --- a/da-sha/sha-topology-optimization-3d/Mesh.h +++ b/da-sha/sha-topology-optimization-3d/Mesh.h @@ -30,7 +30,7 @@ namespace da::sha { return lz_; } - Eigen::VectorXi GetPixelIdx() const { + Eigen::VectorXi GetGlobalIdxOfEleInUse() const { return valid_ele_idx_; } diff --git a/examples/top-thermoelastic-Lshape-condition/CMakeLists.txt b/examples/top-thermoelastic-Lshape-condition/CMakeLists.txt new file mode 100644 index 0000000..b68a990 --- /dev/null +++ b/examples/top-thermoelastic-Lshape-condition/CMakeLists.txt @@ -0,0 +1,9 @@ +set(SUB_PROJECT_NAME top-thermoelastic-Lshape-condition) +#file(GLOB CPP_LIST ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp) +#add_executable(${SUB_PROJECT_NAME} ${CPP_LIST}) +add_executable(${SUB_PROJECT_NAME} main.cpp) +target_link_libraries(${SUB_PROJECT_NAME} ${PROJECT_NAME}_lib) +target_compile_definitions(${SUB_PROJECT_NAME} PUBLIC CONFIG_FILE="${CMAKE_CURRENT_SOURCE_DIR}/config.json") +target_compile_definitions(${SUB_PROJECT_NAME} PUBLIC OUTPUT_DIR="${CMAKE_SOURCE_DIR}/output") +target_compile_definitions(${SUB_PROJECT_NAME} PUBLIC ASSETS_DIR="${CMAKE_SOURCE_DIR}/assets") +target_compile_definitions(${SUB_PROJECT_NAME} PUBLIC DA_CMD) diff --git a/examples/top-thermoelastic-Lshape-condition/config.json b/examples/top-thermoelastic-Lshape-condition/config.json new file mode 100644 index 0000000..fcda67a --- /dev/null +++ b/examples/top-thermoelastic-Lshape-condition/config.json @@ -0,0 +1,67 @@ +{ + "//Topology optimization for L-shape condition": "Lshape Condition", + "material": { + "E": 2.1e11, + "poisson_ratio": 0.3, + "thermal_conductivity":43, "//unit": "W/(m*K)", + "thermal_expansion_coefficient": 1.21e-5, "//unit": "1/K" + }, + "topology": { + "max_loop": 150, + "volfrac": 0.5, + "penal": 3.0, + "r_min": 1.05, + "T_ref": 312, + "T_limit": 3996.52, + "R_E": 28, + "R_lambda": 28, + "R_beta":0 + }, + "model": { + "regular_model": { + "lx":8, + "ly": 30, + "lz": 20 + } + }, + "mechanical_boundary_condition":{ + "NBC": [ + { + "min": [0, 1, 0], + "max": [1, 1, 0], + "val": [0.0,0.0 , -1e8] + } + + ], + "DBC": [ + { + "min": [0, 0, 1], + "max": [1, 0.5, 1], + "dir": [1, 1, 1] + } + + ] + }, + "thermal_boundary_condition": { + "NBC": [ + { + "min": [0.5, 0.2, 0.2], + "max": [0.5, 0.2, 0.2], + "heat_flux": 3.5, "//unit": "W" + }, + { + "min": [0.5, 0.8, 0.2], + "max": [0.5, 0.8, 0.2], + "heat_flux": 3.5, "//unit": "W" + } + + ], + "DBC": [ + { + "min": [0, 0, 1], + "max": [1, 0.5, 1], + "temperature": 312, "//unit": "K" + } + ] + } +} \ No newline at end of file diff --git a/examples/top-thermoelastic-Lshape-condition/main.cpp b/examples/top-thermoelastic-Lshape-condition/main.cpp new file mode 100644 index 0000000..64d7fe9 --- /dev/null +++ b/examples/top-thermoelastic-Lshape-condition/main.cpp @@ -0,0 +1,171 @@ +// +// Created by cflin on 6/9/23. +// +#include +#include +#include "Boundary.h" +#include "Mesh/HeatMesh.h" +#include "Util.h" +#include "ThermoelasticTop3d.h" +#include "FEA/MechanicalLinearFEA.h" +#include "FEA/ThermalLinearFEA.h" +#include "TensorWrapper.h" + +int main() { + using namespace da::sha; + using top::fs_path; + using std::string; + top::fs_path output_dir(OUTPUT_DIR); + top::fs_path config_file(CONFIG_FILE); + top::fs_path assets_dir(ASSETS_DIR); + spdlog::info("Algo read from '{}'", config_file.string()); + spdlog::info("Algo output to '{}'", output_dir.string()); + spdlog::info("asserts dir: '{}'", assets_dir.string()); + + // read json + std::ifstream f(config_file.c_str()); + if (!f) { + spdlog::critical("f open fail!"); + exit(-7); + } + nlohmann::json j_config = nlohmann::json::parse(f); + + // set topology parameters + auto para = std::make_shared(); + para->max_loop = j_config["topology"]["max_loop"]; + para->volfrac = j_config["topology"]["volfrac"]; + para->r_min = j_config["topology"]["r_min"]; + para->penal = j_config["topology"]["penal"]; + para->T_ref = j_config["topology"]["T_ref"]; + para->T_limit = j_config["topology"]["T_limit"]; + para->R_E = j_config["topology"]["R_E"]; + para->R_lambda = j_config["topology"]["R_lambda"]; + para->R_beta = j_config["topology"]["R_beta"]; + + // set material parameters + double E = j_config["material"]["E"]; + double Poisson_ratio = j_config["material"]["poisson_ratio"]; + double thermal_conductivity = j_config["material"]["thermal_conductivity"]; + double thermal_expansion_coefficient = j_config["material"]["thermal_expansion_coefficient"]; + auto material = std::make_shared(E, Poisson_ratio, thermal_conductivity, + thermal_expansion_coefficient); + // set fea + auto sp_mech_fea = std::make_shared(material);// for mechanical + auto sp_thermal_fea = std::make_shared(material);// for thermal + + // set mesh(regular) + int len_x = j_config["model"]["regular_model"]["lx"]; + int len_y = j_config["model"]["regular_model"]["ly"]; + int len_z = j_config["model"]["regular_model"]["lz"]; + // NOTE: USER DEFINE GRID HERE!!! //TODO + std::shared_ptr sp_mech_mesh; + std::shared_ptr sp_thermal_mesh; + std::string condition = j_config["model"]["predefined_condition"]; + if (condition == "L_shape") { + // L-shape condition + spdlog::info("Using {}!", condition); + top::Tensor3d L_shape_model(len_x, len_y, len_z); + L_shape_model.setConstant(1); + for (int k = 0; k < len_z; ++k) { + for (int j = 0; j < len_y; ++j) { + for (int i = 0; i < len_x; ++i) { + if (j > len_y * 0.5 && k > len_z * 0.5) { + L_shape_model(i, j, k)=0; + } + } + } + } + auto &model = L_shape_model; + sp_mech_mesh = std::make_shared(len_x, len_y, len_z, model); + sp_thermal_mesh = std::make_shared(len_x, len_y, len_z, model); + } else { + sp_mech_mesh = std::make_shared(len_x, len_y, len_z); + sp_thermal_mesh = std::make_shared(len_x, len_y, len_z); + } + // initialize Top3d + auto sp_mech_top3d = std::make_shared(para, sp_mech_fea, sp_mech_mesh); + auto sp_thermal_top3d = std::make_shared(para, sp_thermal_fea, sp_thermal_mesh); + + // auxiliary class(Boundary) help to get boundary coordinates + // see the comments in Boundary.h for more information + top::Boundary mech_bdr(sp_mech_mesh); + top::Boundary thermal_bdr(sp_thermal_mesh); + + { + // add Dirichlet boundary, see the comments in Top3d::AddDBC for more information + assert(j_config.count("mechanical_boundary_condition")); + assert(j_config["mechanical_boundary_condition"].count("DBC")); + assert(j_config["mechanical_boundary_condition"].count("NBC")); + int DBCNum = j_config["mechanical_boundary_condition"]["DBC"].size(); + for (int _i = 0; _i < DBCNum; ++_i) { + const auto &DBCI = j_config["mechanical_boundary_condition"]["DBC"][_i]; + Eigen::Vector3d minBBox(DBCI["min"][0], DBCI["min"][1], DBCI["min"][2]); + Eigen::Vector3d maxBBox(DBCI["max"][0], DBCI["max"][1], DBCI["max"][2]); + Eigen::Vector3i dir(DBCI["dir"][0], DBCI["dir"][1], DBCI["dir"][2]); + top::Dir t_dir(minBBox, maxBBox, dir); + + sp_mech_top3d->AddDBC(mech_bdr.GetChosenCoordsByRelativeAlignedBox(t_dir.box), + t_dir.dir); + } + + // add Neumann boundary, see the comments in Top3d::AddNBC for more information + int NBCNum = j_config["mechanical_boundary_condition"]["NBC"].size(); + for (int _i = 0; _i < NBCNum; ++_i) { + const auto &NBCI = j_config["mechanical_boundary_condition"]["NBC"][_i]; + Eigen::Vector3d minBBox(NBCI["min"][0], NBCI["min"][1], NBCI["min"][2]); + Eigen::Vector3d maxBBox(NBCI["max"][0], NBCI["max"][1], NBCI["max"][2]); + Eigen::Vector3d val(NBCI["val"][0], NBCI["val"][1], NBCI["val"][2]); + top::Neu t_neu(minBBox, maxBBox, val); + + Eigen::MatrixXi coords = mech_bdr.GetChosenCoordsByRelativeAlignedBox(t_neu.box); + sp_mech_top3d->AddNBC(coords, t_neu.val / coords.rows()); + } + } + + { + // add Dirichlet boundary, see the comments in Top3d::AddDBC for more information + assert(j_config.count("thermal_boundary_condition")); + assert(j_config["thermal_boundary_condition"].count("DBC")); + assert(j_config["thermal_boundary_condition"].count("NBC")); + int DBCNum = j_config["thermal_boundary_condition"]["DBC"].size(); + for (int _i = 0; _i < DBCNum; ++_i) { + const auto &DBCI = j_config["thermal_boundary_condition"]["DBC"][_i]; + Eigen::Vector3d minBBox(DBCI["min"][0], DBCI["min"][1], DBCI["min"][2]); + Eigen::Vector3d maxBBox(DBCI["max"][0], DBCI["max"][1], DBCI["max"][2]); + top::Dir t_dir(minBBox, maxBBox, Eigen::Vector3i(1, 0, 0)); + + sp_thermal_top3d->AddDBC(mech_bdr.GetChosenCoordsByRelativeAlignedBox(t_dir.box), + t_dir.dir, DBCI["temperature"]); + } + + // add Neumann boundary, see the comments in Top3d::AddNBC for more information + int NBCNum = j_config["thermal_boundary_condition"]["NBC"].size(); + for (int _i = 0; _i < NBCNum; ++_i) { + const auto &NBCI = j_config["thermal_boundary_condition"]["NBC"][_i]; + Eigen::Vector3d minBBox(NBCI["min"][0], NBCI["min"][1], NBCI["min"][2]); + Eigen::Vector3d maxBBox(NBCI["max"][0], NBCI["max"][1], NBCI["max"][2]); + Eigen::Vector3d val(NBCI["heat_flux"], 0, 0); + top::Neu t_neu(minBBox, maxBBox, val); + + Eigen::MatrixXi coords = mech_bdr.GetChosenCoordsByRelativeAlignedBox(t_neu.box); + sp_thermal_top3d->AddNBC(coords, t_neu.val / coords.rows()); + } + } + + + // init thermoelastic top3d + top::ThermoelasticTop3d ther_top3d(sp_mech_top3d, sp_thermal_top3d); + + // process topology optimization + top::Tensor3d ten_rho = ther_top3d.TopOptMainLoop(); + + // extract txt or vtk + write_tensor3d(output_dir / "txt" / "top_mach_regular_field_matrix.txt", ten_rho, + sp_mech_mesh->GetOrigin(), + sp_mech_mesh->GetOrigin() + sp_mech_mesh->GetLenBox()); + WriteTensorToVtk(output_dir / "vtk" / "top_mach_regular_field_matrix.vtk", ten_rho, + sp_mech_mesh); + WriteUToVtk(output_dir / "vtk" / "top_thermoelastic_regular_U.vtk", ther_top3d.GetU(), + sp_mech_mesh); + +} diff --git a/examples/top-thermoelastic-regular-model/config.json b/examples/top-thermoelastic-regular-model/config.json index c8f1fa6..30e1d8b 100644 --- a/examples/top-thermoelastic-regular-model/config.json +++ b/examples/top-thermoelastic-regular-model/config.json @@ -6,59 +6,61 @@ "thermal_expansion_coefficient": 1.21e-5, "unit": "1/K" }, "topology": { - "max_loop": 200, - "volfrac": 0.4, + "max_loop": 150, + "volfrac": 0.5, "penal": 3.0, - "r_min": 1.1, - "T_ref": 0, - "T_limit": 26, + "r_min": 1.05, + "T_ref": 312, + "T_limit": 3996.52, "R_E": 28, "R_lambda": 28, "R_beta":0 }, "model": { + "predefined_condition": "L_shape", "regular_model": { - "lx": 1, - "ly": 50, - "lz": 35 - }, - "visual_model":"2-box-refined.obj" + "lx":8, + "ly": 30, + "lz": 20 + } }, "mechanical_boundary_condition":{ "NBC": [ { - "min": [0, 0.44, 0], - "max": [1, 0.56, 0], - "val": [0.0, 0.0, -8.4e8] + "min": [0, 1, 0], + "max": [1, 1, 0], + "val": [0.0,0.0 , -1e8] } ], "DBC": [ { - "min": [0, 0, 0], - "max": [1, 0, 1], - "dir": [1, 1, 1] - }, - { - "min": [0, 1, 0], - "max": [1, 1, 1], + "min": [0, 0, 1], + "max": [1, 0.5, 1], "dir": [1, 1, 1] } + ] }, "thermal_boundary_condition": { "NBC": [ { - "min": [0, 1, 0.3], - "max": [1, 1, 0.7], - "heat_flux": 5, "unit": "W" + "min": [0.5, 0.2, 0.2], + "max": [0.5, 0.2, 0.2], + "heat_flux": 3.5,"back":12.35, "unit": "W" + }, + { + "min": [0.5, 0.8, 0.2], + "max": [0.5, 0.8, 0.2], + "heat_flux": 3.5,"back":12.35, "unit": "W" } + ], "DBC": [ { - "min": [0, 0, 0], - "max": [1, 0, 1], - "temperature": 0, "unit": "K" + "min": [0, 0, 1], + "max": [1, 0.5, 1], + "temperature": 312, "unit": "K" } ] } diff --git a/examples/top-thermoelastic-regular-model/config_case2.json b/examples/top-thermoelastic-regular-model/config_case2.json new file mode 100644 index 0000000..d9e83ee --- /dev/null +++ b/examples/top-thermoelastic-regular-model/config_case2.json @@ -0,0 +1,65 @@ +{ + "material": { + "E": 2.1e11, + "poisson_ratio": 0.3, + "thermal_conductivity":43, "unit": "W/(m*K)", + "thermal_expansion_coefficient": 1.21e-5, "unit": "1/K" + }, + "topology": { + "max_loop": 200, + "volfrac": 0.4, + "penal": 3.0, + "r_min": 1.1, + "T_ref": 0, + "T_limit": 28, + "R_E": 28, + "R_lambda": 28, + "R_beta":0 + }, + "model": { + "regular_model": { + "lx": 1, + "ly": 50, + "lz": 35 + }, + "visual_model":"2-box-refined.obj" + }, + "mechanical_boundary_condition":{ + "NBC": [ + { + "min": [0, 0.44, 0], + "max": [1, 0.56, 0], + "val": [0.0, 0.0, -8.4e8] + } + + ], + "DBC": [ + { + "min": [0, 0, 0], + "max": [1, 0, 1], + "dir": [1, 1, 1] + }, + { + "min": [0, 1, 0], + "max": [1, 1, 1], + "dir": [1, 1, 1] + } + ] + }, + "thermal_boundary_condition": { + "NBC": [ + { + "min": [0, 1, 0.3], + "max": [1, 1, 0.7], + "heat_flux": 5, "unit": "W" + } + ], + "DBC": [ + { + "min": [0, 0, 0], + "max": [1, 0, 1], + "temperature": 0, "unit": "K" + } + ] + } +} \ No newline at end of file diff --git a/examples/top-thermoelastic-regular-model/main.cpp b/examples/top-thermoelastic-regular-model/main.cpp index d40e2a2..64d7fe9 100644 --- a/examples/top-thermoelastic-regular-model/main.cpp +++ b/examples/top-thermoelastic-regular-model/main.cpp @@ -9,6 +9,7 @@ #include "ThermoelasticTop3d.h" #include "FEA/MechanicalLinearFEA.h" #include "FEA/ThermalLinearFEA.h" +#include "TensorWrapper.h" int main() { using namespace da::sha; @@ -34,19 +35,20 @@ int main() { para->max_loop = j_config["topology"]["max_loop"]; para->volfrac = j_config["topology"]["volfrac"]; para->r_min = j_config["topology"]["r_min"]; - para->penal=j_config["topology"]["penal"]; - para->T_ref=j_config["topology"]["T_ref"]; - para->T_limit=j_config["topology"]["T_limit"]; - para->R_E=j_config["topology"]["R_E"]; - para->R_lambda=j_config["topology"]["R_lambda"]; - para->R_beta=j_config["topology"]["R_beta"]; + para->penal = j_config["topology"]["penal"]; + para->T_ref = j_config["topology"]["T_ref"]; + para->T_limit = j_config["topology"]["T_limit"]; + para->R_E = j_config["topology"]["R_E"]; + para->R_lambda = j_config["topology"]["R_lambda"]; + para->R_beta = j_config["topology"]["R_beta"]; // set material parameters double E = j_config["material"]["E"]; double Poisson_ratio = j_config["material"]["poisson_ratio"]; double thermal_conductivity = j_config["material"]["thermal_conductivity"]; - double thermal_expansion_coefficient=j_config["material"]["thermal_expansion_coefficient"]; - auto material = std::make_shared(E, Poisson_ratio, thermal_conductivity,thermal_expansion_coefficient); + double thermal_expansion_coefficient = j_config["material"]["thermal_expansion_coefficient"]; + auto material = std::make_shared(E, Poisson_ratio, thermal_conductivity, + thermal_expansion_coefficient); // set fea auto sp_mech_fea = std::make_shared(material);// for mechanical auto sp_thermal_fea = std::make_shared(material);// for thermal @@ -55,8 +57,31 @@ int main() { int len_x = j_config["model"]["regular_model"]["lx"]; int len_y = j_config["model"]["regular_model"]["ly"]; int len_z = j_config["model"]["regular_model"]["lz"]; - auto sp_mech_mesh = std::make_shared(len_x, len_y, len_z); - auto sp_thermal_mesh = std::make_shared(len_x, len_y, len_z); + // NOTE: USER DEFINE GRID HERE!!! //TODO + std::shared_ptr sp_mech_mesh; + std::shared_ptr sp_thermal_mesh; + std::string condition = j_config["model"]["predefined_condition"]; + if (condition == "L_shape") { + // L-shape condition + spdlog::info("Using {}!", condition); + top::Tensor3d L_shape_model(len_x, len_y, len_z); + L_shape_model.setConstant(1); + for (int k = 0; k < len_z; ++k) { + for (int j = 0; j < len_y; ++j) { + for (int i = 0; i < len_x; ++i) { + if (j > len_y * 0.5 && k > len_z * 0.5) { + L_shape_model(i, j, k)=0; + } + } + } + } + auto &model = L_shape_model; + sp_mech_mesh = std::make_shared(len_x, len_y, len_z, model); + sp_thermal_mesh = std::make_shared(len_x, len_y, len_z, model); + } else { + sp_mech_mesh = std::make_shared(len_x, len_y, len_z); + sp_thermal_mesh = std::make_shared(len_x, len_y, len_z); + } // initialize Top3d auto sp_mech_top3d = std::make_shared(para, sp_mech_fea, sp_mech_mesh); auto sp_thermal_top3d = std::make_shared(para, sp_thermal_fea, sp_thermal_mesh); @@ -79,7 +104,8 @@ int main() { Eigen::Vector3i dir(DBCI["dir"][0], DBCI["dir"][1], DBCI["dir"][2]); top::Dir t_dir(minBBox, maxBBox, dir); - sp_mech_top3d->AddDBC(mech_bdr.GetChosenCoordsByRelativeAlignedBox(t_dir.box), t_dir.dir); + sp_mech_top3d->AddDBC(mech_bdr.GetChosenCoordsByRelativeAlignedBox(t_dir.box), + t_dir.dir); } // add Neumann boundary, see the comments in Top3d::AddNBC for more information @@ -91,8 +117,8 @@ int main() { Eigen::Vector3d val(NBCI["val"][0], NBCI["val"][1], NBCI["val"][2]); top::Neu t_neu(minBBox, maxBBox, val); - Eigen::MatrixXi coords=mech_bdr.GetChosenCoordsByRelativeAlignedBox(t_neu.box); - sp_mech_top3d->AddNBC(coords, t_neu.val/coords.rows()); + Eigen::MatrixXi coords = mech_bdr.GetChosenCoordsByRelativeAlignedBox(t_neu.box); + sp_mech_top3d->AddNBC(coords, t_neu.val / coords.rows()); } } @@ -106,9 +132,10 @@ int main() { const auto &DBCI = j_config["thermal_boundary_condition"]["DBC"][_i]; Eigen::Vector3d minBBox(DBCI["min"][0], DBCI["min"][1], DBCI["min"][2]); Eigen::Vector3d maxBBox(DBCI["max"][0], DBCI["max"][1], DBCI["max"][2]); - top::Dir t_dir(minBBox, maxBBox, Eigen::Vector3i(1,0,0)); + top::Dir t_dir(minBBox, maxBBox, Eigen::Vector3i(1, 0, 0)); - sp_thermal_top3d->AddDBC(mech_bdr.GetChosenCoordsByRelativeAlignedBox(t_dir.box), t_dir.dir,DBCI["temperature"]); + sp_thermal_top3d->AddDBC(mech_bdr.GetChosenCoordsByRelativeAlignedBox(t_dir.box), + t_dir.dir, DBCI["temperature"]); } // add Neumann boundary, see the comments in Top3d::AddNBC for more information @@ -117,11 +144,11 @@ int main() { const auto &NBCI = j_config["thermal_boundary_condition"]["NBC"][_i]; Eigen::Vector3d minBBox(NBCI["min"][0], NBCI["min"][1], NBCI["min"][2]); Eigen::Vector3d maxBBox(NBCI["max"][0], NBCI["max"][1], NBCI["max"][2]); - Eigen::Vector3d val( NBCI["heat_flux"],0,0); + Eigen::Vector3d val(NBCI["heat_flux"], 0, 0); top::Neu t_neu(minBBox, maxBBox, val); - Eigen::MatrixXi coords=mech_bdr.GetChosenCoordsByRelativeAlignedBox(t_neu.box); - sp_thermal_top3d->AddNBC(coords, t_neu.val/coords.rows()); + Eigen::MatrixXi coords = mech_bdr.GetChosenCoordsByRelativeAlignedBox(t_neu.box); + sp_thermal_top3d->AddNBC(coords, t_neu.val / coords.rows()); } } @@ -133,9 +160,12 @@ int main() { top::Tensor3d ten_rho = ther_top3d.TopOptMainLoop(); // extract txt or vtk - write_tensor3d(output_dir / "txt" / "top_mach_regular_field_matrix.txt", ten_rho, sp_mech_mesh->GetOrigin(), + write_tensor3d(output_dir / "txt" / "top_mach_regular_field_matrix.txt", ten_rho, + sp_mech_mesh->GetOrigin(), sp_mech_mesh->GetOrigin() + sp_mech_mesh->GetLenBox()); - WriteTensorToVtk(output_dir / "vtk" / "top_mach_regular_field_matrix.vtk", ten_rho, sp_mech_mesh); - WriteUToVtk(output_dir / "vtk" / "top_thermoelastic_regular_U.vtk",ther_top3d.GetU(),sp_mech_mesh); + WriteTensorToVtk(output_dir / "vtk" / "top_mach_regular_field_matrix.vtk", ten_rho, + sp_mech_mesh); + WriteUToVtk(output_dir / "vtk" / "top_thermoelastic_regular_U.vtk", ther_top3d.GetU(), + sp_mech_mesh); } diff --git a/src/Mesh/HeatMesh.h b/src/Mesh/HeatMesh.h index 4647cdf..0977762 100644 --- a/src/Mesh/HeatMesh.h +++ b/src/Mesh/HeatMesh.h @@ -14,8 +14,13 @@ namespace da::sha { public: HeatMesh(int len_x, int len_y, int len_z) : Mesh(len_x, len_y, len_z, 1) { } - }; + HeatMesh(int len_x, int len_y, int len_z, const Tensor3d &user_defined_grid) + : + Mesh(len_x, len_y, + len_z, + user_defined_grid, 1) {} + }; } // top } // da::sha diff --git a/src/Mesh/IrregularMesh.cpp b/src/Mesh/IrregularMesh.cpp index f82caf1..723f65c 100644 --- a/src/Mesh/IrregularMesh.cpp +++ b/src/Mesh/IrregularMesh.cpp @@ -128,7 +128,7 @@ namespace da::sha { } } chosen_ele_id_=Eigen::Map(v_chosen_ele_id.data(),v_chosen_ele_id.size()); - init_ele_rho_=Eigen::Map(v_rho_non_empty.data(),v_rho_non_empty.size()); + init_rho_of_ele_in_use_=Eigen::Map(v_rho_non_empty.data(), v_rho_non_empty.size()); num_pixel_ = cnt_pixel; // get_num_node_pixel && fill node_id @@ -175,12 +175,12 @@ namespace da::sha { } } // fill valid_ele_idx_ - valid_ele_idx_.resize(num_pixel_); + global_idx_of_ele_in_use.resize(num_pixel_); Tensor3i ten_tmp_col = ten_ele_coord2ele_id_.reshape(Eigen::array{num_ele_, 1, 1}); cnt_pixel = 0; for (int i = 0; i < ten_tmp_col.size(); ++i) { if (ten_tmp_col(i, 0, 0) != -1) { - valid_ele_idx_(cnt_pixel) = i; + global_idx_of_ele_in_use(cnt_pixel) = i; ++cnt_pixel; } } diff --git a/src/Mesh/Mesh.cpp b/src/Mesh/Mesh.cpp index 7c8093c..5b33efd 100644 --- a/src/Mesh/Mesh.cpp +++ b/src/Mesh/Mesh.cpp @@ -1,11 +1,3 @@ -/* - * @Author: lab pc yjxkwp@foxmail.com - * @Date: 2023-04-28 22:24:13 - * @LastEditors: lab pc yjxkwp@foxmail.com - * @LastEditTime: 2023-04-29 12:17:53 - * @FilePath: /designauto/da-sha/sha-topology-optimization-3d/Mesh.cpp - * @Description: 这是默认设置,请设置`customMade`, 打开koroFileHeader查看配置 进行设置: https://github.com/OBKoro1/koro1FileHeader/wiki/%E9%85%8D%E7%BD%AE - */ // // Created by cflin on 4/20/23. // @@ -15,11 +7,141 @@ namespace da::sha { namespace top { - Mesh::Mesh(int len_x, int len_y, int len_z,int dofs_each_node) : DOFS_EACH_NODE(dofs_each_node), lx_(len_x), ly_(len_y), lz_(len_z), - num_node_((lx_ + 1) * (ly_ + 1) * (lz_ + 1)), - num_ele_(lx_ * ly_ * lz_), origin_(0,0,0), len_pixel_(1.0), len_box_(len_x,len_y,len_z){ - valid_ele_idx_ = Eigen::VectorXi::LinSpaced(num_ele_, 0, num_ele_ - 1); - init_ele_rho_=Eigen::VectorXd::Ones(num_ele_); + Mesh::Mesh(int len_x, int len_y, int len_z, const Tensor3d &user_defined_grid, int + dofs_each_node) : DOFS_EACH_NODE(dofs_each_node), lx_(len_x), ly_(len_y), lz_(len_z), + num_node_((lx_ + 1) * (ly_ + 1) * (lz_ + 1)), + num_ele_(lx_ * ly_ * lz_), + origin_(0, 0, 0), + len_pixel_(1.0), + len_box_(len_x, len_y, len_z) { + + if (!(user_defined_grid.dimension(0) == lx_ && user_defined_grid.dimension(1) == ly_ + && user_defined_grid.dimension(2) == lz_)) { + spdlog::critical("pls check dimension!"); + exit(-1); + }; + // valid_ele_idx_ + // init_rho_of_ele_in_use_ + + ten_ele_coord2ele_id_ = Tensor3i(lx_, ly_, lz_); + ten_ele_coord2ele_id_.setConstant(-1); + ten_node_coord2node_id_ = Tensor3i(lx_ + 1, ly_ + 1, lz_ + 1); + ten_node_coord2node_id_.setConstant(-1); + + 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(); + + // count ele really in use && fill ten_ele_coord2ele_id_ && set flg for ten_node_coord2node_id_ + const int NOT_EMPTY_FLG = 1; + int cnt_ele_in_use = 0; +// std::vector v_chosen_ele_id; + std::vector v_rho_non_empty; + for (int k = 0; k < lz_; ++k) { + for (int j = 0; j < ly_; ++j) { + for (int i = 0; i < lx_; ++i) { + if (user_defined_grid(i, j, k) > 0) { + ten_ele_coord2ele_id_(i, j, k) = cnt_ele_in_use; +// v_chosen_ele_id.push_back(cnt_ele_in_use); + ++cnt_ele_in_use; + v_rho_non_empty.push_back(user_defined_grid(i, j, k)); + 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)) = NOT_EMPTY_FLG; + } + } + } + } + } +// chosen_ele_id_=Eigen::VectorXi::LinSpaced(cnt_ele_in_use,0,cnt_ele_in_use); + init_rho_of_ele_in_use_ = Eigen::Map(v_rho_non_empty.data(), + v_rho_non_empty.size()); + num_ele_in_use_ = cnt_ele_in_use; + + // get_num_node_pixel && fill node_id + int cnt_node_pixel = 0; + std::vector v_node_coords; + 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) == NOT_EMPTY_FLG) { + ten_node_coord2node_id_(i, j, k) = cnt_node_pixel; + v_node_coords.push_back({i, j, k}); + ++cnt_node_pixel; + } + } + } + } + num_node_in_use_ = cnt_node_pixel; + // fill mat_node_id 2 node_coord + mat_node_id2node_coord_ = Eigen::MatrixXi(v_node_coords.size(), 3); + for (int i = 0; i < v_node_coords.size(); ++i) { + mat_node_id2node_coord_.row(i) = v_node_coords[i]; + } + + + // fill mat_ele_id2dofs_ + mat_ele_id2dofs_.resize(num_ele_in_use_, NUM_NODES_EACH_ELE * DOFS_EACH_NODE); + mat_ele_id2node_id_.resize(num_ele_in_use_, NUM_NODES_EACH_ELE); + 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 node_ids = ten_node_coord2node_id_(world_node_coords); + mat_ele_id2node_id_.row(cur_ele_id) = node_ids; + for (int nodi = 0; nodi < NUM_NODES_EACH_ELE; ++nodi) { + for (int dofi = 0; dofi < DOFS_EACH_NODE; ++dofi) { + mat_ele_id2dofs_(cur_ele_id, DOFS_EACH_NODE * nodi + dofi) = + DOFS_EACH_NODE * node_ids(nodi) + dofi; + } + } + } + } + } + + // fill valid_ele_idx_ + global_idx_of_ele_in_use.resize(num_ele_in_use_); + cnt_ele_in_use = 0; + for (int k = 0; k < lz_; ++k) { + for (int j = 0; j < ly_; ++j) { + for (int i = 0; i < lx_; ++i) { + int idx=k*ly_*lx_+ j *lx_ +i; + if (ten_ele_coord2ele_id_(i, j, k) != -1) { + global_idx_of_ele_in_use(cnt_ele_in_use++) =idx; + } + } + } + } + assert(cnt_ele_in_use == num_ele_in_use_); + + } + + Mesh::Mesh(int len_x, int len_y, int len_z, int dofs_each_node) + : + DOFS_EACH_NODE(dofs_each_node), lx_(len_x), ly_(len_y), lz_(len_z), + num_node_((lx_ + 1) * (ly_ + 1) * (lz_ + 1)), + num_ele_(lx_ * ly_ * lz_), origin_(0, 0, 0), len_pixel_(1.0), + len_box_(len_x, len_y, len_z) { + num_node_in_use_ = num_node_; + num_ele_in_use_ = num_ele_; + + global_idx_of_ele_in_use = Eigen::VectorXi::LinSpaced(num_ele_, 0, num_ele_); + init_rho_of_ele_in_use_ = Eigen::VectorXd::Ones(num_ele_); // ten_node_coord2node_id_ = Tensor3i(num_node_, 1, 1); // for (int i = 0; i < num_node_; ++i) { @@ -30,20 +152,20 @@ namespace da::sha { // get node_coord <---> node_id int cnt_node_pixel = 0; std::vector v_node_coords; - ten_node_coord2node_id_=Tensor3i(lx_ + 1, ly_ + 1, lz_ + 1); + ten_node_coord2node_id_ = Tensor3i(lx_ + 1, ly_ + 1, lz_ + 1); for (int k = 0; k < lz_ + 1; ++k) { for (int j = 0; j < ly_ + 1; ++j) { for (int i = 0; i < lx_ + 1; ++i) { ten_node_coord2node_id_(i, j, k) = cnt_node_pixel; - v_node_coords.push_back({i,j,k}); + v_node_coords.push_back({i, j, k}); ++cnt_node_pixel; } } } // fill mat_node_id 2 node_coord - mat_node_id2node_coord_=Eigen::MatrixXi(v_node_coords.size(),3); - for(int i=0;i{lx_, ly_, lz_}); + ten_ele_coord2ele_id_ = ten_ele_coord2ele_id_.reshape( + Eigen::array{lx_, ly_, lz_}); mat_ele_id2dofs_.resize(num_ele_, NUM_NODES_EACH_ELE * DOFS_EACH_NODE); - mat_ele_id2node_id_.resize(num_ele_,NUM_NODES_EACH_ELE); + mat_ele_id2node_id_.resize(num_ele_, NUM_NODES_EACH_ELE); static const Eigen::MatrixXi delta_coord = (Eigen::MatrixXi(8, 3) << 0, 0, 0, 1, 0, 0, @@ -69,13 +192,15 @@ namespace da::sha { 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); - Eigen::MatrixXi world_node_coords = delta_coord.rowwise() + Eigen::RowVector3i(i, j, k); + 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 node_ids = ten_node_coord2node_id_(world_node_coords); - mat_ele_id2node_id_.row(cur_ele_id)=node_ids; + mat_ele_id2node_id_.row(cur_ele_id) = node_ids; for (int nodi = 0; nodi < NUM_NODES_EACH_ELE; ++nodi) { - for(int dofi=0; dofi < DOFS_EACH_NODE; ++dofi){ - mat_ele_id2dofs_(cur_ele_id, DOFS_EACH_NODE * nodi + dofi)= DOFS_EACH_NODE * node_ids(nodi) + dofi; + for (int dofi = 0; dofi < DOFS_EACH_NODE; ++dofi) { + mat_ele_id2dofs_(cur_ele_id, DOFS_EACH_NODE * nodi + dofi) = + DOFS_EACH_NODE * node_ids(nodi) + dofi; } } } diff --git a/src/Mesh/Mesh.h b/src/Mesh/Mesh.h index 27f3e1c..9544814 100644 --- a/src/Mesh/Mesh.h +++ b/src/Mesh/Mesh.h @@ -15,7 +15,8 @@ namespace da::sha { class Mesh { public: Mesh():origin_(0,0,0),len_pixel_(1.0){}; - + Mesh(int len_x, int len_y, int len_z,const Tensor3d& user_defined_grid, int + dofs_each_node=3); Mesh(int len_x, int len_y, int len_z,int dofs_each_node=3); int GetLx() const { @@ -30,21 +31,28 @@ namespace da::sha { return lz_; } - Eigen::VectorXi GetPixelIdx() const { - return valid_ele_idx_; + Eigen::VectorXi GetGlobalIdxOfEleInUse() const { + return global_idx_of_ele_in_use; } virtual int GetNumDofs() const { - return num_node_ * DOFS_EACH_NODE; + return num_node_in_use_ * DOFS_EACH_NODE; } virtual int GetNumEles() const { - return num_ele_; + return num_ele_in_use_; } virtual int GetNumNodes() const { - return num_node_; + return num_node_in_use_; } +// virtual int GetNumEleInUse() const { +// return num_ele_in_use_; +// } +// +// virtual int GetNumNodesInUse() const { +// return num_node_in_use_; +// } Eigen::MatrixXi GetEleId2DofsMap() const { return mat_ele_id2dofs_; @@ -84,7 +92,7 @@ namespace da::sha { return NUM_NODES_EACH_ELE*DOFS_EACH_NODE; } Eigen::VectorXd GetInitEleRho() const{ - return init_ele_rho_; + return init_rho_of_ele_in_use_; } Eigen::VectorXi MapNodeCoord2Id(const Eigen::MatrixXi &node_coords) const { @@ -114,8 +122,10 @@ namespace da::sha { int ly_; int lz_; int num_ele_; + int num_ele_in_use_; int num_node_; - Eigen::VectorXi valid_ele_idx_; + int num_node_in_use_; + Eigen::VectorXi global_idx_of_ele_in_use; Tensor3i ten_node_coord2node_id_; Eigen::MatrixXi mat_node_id2node_coord_;// num_node x 3 Tensor3i ten_ele_coord2ele_id_; @@ -127,7 +137,7 @@ namespace da::sha { Eigen::Vector3d origin_; Eigen::Vector3d len_box_; double len_pixel_; - Eigen::VectorXd init_ele_rho_; + Eigen::VectorXd init_rho_of_ele_in_use_; }; diff --git a/src/TensorWrapper.h b/src/TensorWrapper.h index 923afa7..40c39e7 100644 --- a/src/TensorWrapper.h +++ b/src/TensorWrapper.h @@ -39,6 +39,20 @@ namespace da::sha { } }; + template + std::ostream& operator<<(std::ostream & os,const TensorWrapper& t3){ + for(int k=0;ksp_mesh_; auto &sp_para_ = sp_mech_top3d_->sp_para_; Eigen::VectorXd xPhys_col(sp_mesh_->GetNumEles()); - xPhys_col.setConstant(sp_para_->volfrac/sp_para_->volfrac); + xPhys_col.setConstant(sp_para_->volfrac / sp_para_->volfrac); bool flg_chosen = false; - Eigen::VectorXi chosen_ele_id; -// Eigen::VectorXi chosen_ele_id(sp_mesh_->GetChosenEleIdx()); -// bool flg_chosen = chosen_ele_id.size() != 0; -// if (!flg_chosen) { -// // no part chosen -// xPhys_col.setConstant(sp_para_->volfrac); -// } else { -// // pick chosen part to sim -// xPhys_col = sp_mesh_->GetInitEleRho(); -// xPhys_col(chosen_ele_id).setConstant(sp_para_->volfrac); -// } + Eigen::VectorXi chosen_ele_id(sp_mesh_->GetChosenEleIdx()); +// bool flg_chosen = chosen_ele_id.size() != 0; +// if (!flg_chosen) { +// // no part chosen +// xPhys_col.setConstant(sp_para_->volfrac); +// } else { +// // pick chosen part to sim +// xPhys_col = sp_mesh_->GetInitEleRho(); +// xPhys_col(chosen_ele_id).setConstant(sp_para_->volfrac); +// } int loop = 0; double change = 1.0; @@ -38,7 +37,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { Eigen::VectorXd ele_to_write = Eigen::VectorXd::Zero(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz()); - Eigen::VectorXi pixel_idx = sp_mesh_->GetPixelIdx(); + Eigen::VectorXi global_idx_of_ele_in_use = sp_mesh_->GetGlobalIdxOfEleInUse(); // dofs of limited T struct LimitedDof { @@ -46,8 +45,10 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { int idx_of_load_dof; int idx_in_ele; - LimitedDof(int dof, int idx_of_load_dof, int idx_in_ele) : dof(dof), idx_of_load_dof(idx_of_load_dof), - idx_in_ele(idx_in_ele) {} + LimitedDof(int dof, int idx_of_load_dof, int idx_in_ele) + : + dof(dof), idx_of_load_dof(idx_of_load_dof), + idx_in_ele(idx_in_ele) {} }; std::map> map_ele2Limit; std::vector v_dof(sp_thermal_top3d_->set_dofs_to_load.begin(), @@ -102,7 +103,8 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { return CalDRDrho_Vec(vec_rho, sp_para_->R_E) * (E0_m - E_min); }; auto CalLambda_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { - return lambda_min + CalR_Vec(vec_rho, sp_para_->R_lambda).array() * (lambda0 - lambda_min); + return lambda_min + + CalR_Vec(vec_rho, sp_para_->R_lambda).array() * (lambda0 - lambda_min); }; auto CalDlambdaDrho = [&](double rho) { return CalDRDrho(rho, sp_para_->R_lambda) * (lambda0 - lambda_min); @@ -202,9 +204,11 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { // double beta_rho = CalBeta(xPhys_col(i)); // F_th(dofs_m) += beta_rho * (Te - sp_mech_top3d_->sp_para_->T_ref) * Inted_; Eigen::VectorXd ele_dFth_drho = - CalDBetaDrho(xPhys_col(i)) * (Te - sp_mech_top3d_->sp_para_->T_ref) * Inted_;// 24x1 + CalDBetaDrho(xPhys_col(i)) * (Te - sp_mech_top3d_->sp_para_->T_ref) * + Inted_;// 24x1 assert(ele_dFth_drho.size() == 24); - v_dFth_drho(Eigen::seqN(i * ele_dFth_drho.size(), ele_dFth_drho.size())) = ele_dFth_drho; + v_dFth_drho( + Eigen::seqN(i * ele_dFth_drho.size(), ele_dFth_drho.size())) = ele_dFth_drho; } auto v_dFth_drho_tri = Vec2Triplet(i_dFth_drho_, j_dFth_drho_, v_dFth_drho); dFth_drho.setFromTriplets(v_dFth_drho_tri.begin(), v_dFth_drho_tri.end()); @@ -221,9 +225,11 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { // F_th(dofs_m) += beta_rho * (Te - sp_mech_top3d_->sp_para_->T_ref) * Inted_; Eigen::MatrixXd ele_dFth_dT = Eigen::VectorXd::Ones(sp_thermal_top3d_->sp_mesh_->Get_DOFS_EACH_ELE()) * 1.0 / - sp_thermal_top3d_->sp_mesh_->Get_DOFS_EACH_ELE() * beta_rho * Inted_.transpose(); + sp_thermal_top3d_->sp_mesh_->Get_DOFS_EACH_ELE() * beta_rho * + Inted_.transpose(); assert(ele_dFth_dT.rows() == 8 && ele_dFth_dT.cols() == 24); - v_dFth_dT(Eigen::seqN(i * ele_dFth_dT.size(), ele_dFth_dT.size())) = ele_dFth_dT.reshaped(); + v_dFth_dT(Eigen::seqN(i * ele_dFth_dT.size(), + ele_dFth_dT.size())) = ele_dFth_dT.reshaped(); } auto v_dFth_dT_tri = Vec2Triplet(i_dFth_dT_, j_dFth_dT_, v_dFth_dT); dFth_dT.setFromTriplets(v_dFth_dT_tri.begin(), v_dFth_dT_tri.end()); @@ -248,7 +254,8 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { Eigen::VectorXd lambda_t_e = lambda_t(dofs_in_ele_i); ce_th(i) = lambda_t_e.transpose() * sp_thermal_top3d_->Ke_ * Te; } - Eigen::VectorXd lambda_t_Mul_dKt_drho_Mul_T = CalDlambdaDrho_Vec(xPhys_col).array() * ce_th.array(); + Eigen::VectorXd lambda_t_Mul_dKt_drho_Mul_T = + CalDlambdaDrho_Vec(xPhys_col).array() * ce_th.array(); // dc_drho // Eigen::VectorXd dc_drho = lambda_t_Mul_dKt_drho_Mul_T + // lambda_m_Mul_dKm_drho_Mul_U + @@ -279,7 +286,8 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { } return sp_dKth_Mul_T; }; - auto CalDTi_Drhoj = [&](int Tdof, const Eigen::SparseVector &sp_dKth_Mul_T) -> double { + auto CalDTi_Drhoj = [&](int Tdof, + const Eigen::SparseVector &sp_dKth_Mul_T) -> double { Eigen::VectorXd Li = Eigen::VectorXd::Zero(sp_dKth_Mul_T.rows()); Li(Tdof) = -1; for (auto dof_value: sp_thermal_top3d_->v_dofs_to_set) { @@ -320,16 +328,18 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { Eigen::VectorXd variables_tmp = flg_chosen ? xPhys_col(chosen_ele_id) : xPhys_col; double f0val = c; Eigen::VectorXd df0dx = flg_chosen - ? dc_dx(chosen_ele_id).eval() / dc_dx(chosen_ele_id).cwiseAbs().maxCoeff() * + ? dc_dx(chosen_ele_id).eval() / + dc_dx(chosen_ele_id).cwiseAbs().maxCoeff() * SENSITIVITY_SCALE_COEF : dc_dx / dc_dx.cwiseAbs().maxCoeff() * SENSITIVITY_SCALE_COEF; // double fval = v - num_variables * sp_para_->volfrac; - Eigen::VectorXd fval = (Eigen::VectorXd(num_constraints) << (v / num_variables - sp_para_->volfrac) + Eigen::VectorXd fval = + (Eigen::VectorXd(num_constraints) << (v / num_variables - sp_para_->volfrac) #ifndef MECH_ONLY - , T(v_dof).array() / sp_para_->T_limit - 1 + , T(v_dof).array() / sp_para_->T_limit - 1 #endif - ).finished() * SENSITIVITY_SCALE_COEF; + ).finished() * SENSITIVITY_SCALE_COEF; // Eigen::VectorXd dfdx = flg_chosen ? dv(chosen_ele_id) : dv; Eigen::MatrixXd dfdx = (Eigen::MatrixXd(num_variables, num_constraints) << 1.0 / num_variables * dv @@ -354,11 +364,12 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { xPhys_col = variables_tmp; } - spdlog::critical("Iter: {:3d}, Comp: {:.3e}, Vol: {:.2f}, Change: {:f}", loop, c, v, change); + spdlog::critical("Iter: {:3d}, Comp: {:.3e}, Vol: {:.2f}, Change: {:f}", loop, c, v, + change); std::cout << fval.transpose() << std::endl; #ifdef WRITE_TENSOR_IN_LOOP // extract vtk - ele_to_write(pixel_idx) = xPhys_col; + ele_to_write(idx_of_ele_in_use) = xPhys_col; Tensor3d ten_xPhys_to_write(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz(), 1, 1); for (int i = 0; i < ele_to_write.size(); ++i) { ten_xPhys_to_write(i, 0, 0) = ele_to_write(i); @@ -375,12 +386,14 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { // set 0 to rho of unchosen part assert(xPhys_col.size()); Eigen::VectorXi continue_idx = - Eigen::VectorXi::LinSpaced(xPhys_col.size(), 0, xPhys_col.size() - 1); - Eigen::VectorXi unchosen_idx = flg_chosen ? SetDifference(continue_idx, chosen_ele_id) : Eigen::VectorXi(); + Eigen::VectorXi::LinSpaced(xPhys_col.size(), 0, xPhys_col.size()); + Eigen::VectorXi unchosen_idx = flg_chosen ? SetDifference(continue_idx, chosen_ele_id) + : Eigen::VectorXi(); { xPhys_col(unchosen_idx).setZero(); - ele_to_write(pixel_idx) = xPhys_col; - Tensor3d ten_xPhys_to_write(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz(), 1, 1); + ele_to_write(global_idx_of_ele_in_use) = xPhys_col; + Tensor3d ten_xPhys_to_write(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz(), 1, + 1); for (int i = 0; i < ele_to_write.size(); ++i) { ten_xPhys_to_write(i, 0, 0) = ele_to_write(i); } @@ -391,8 +404,9 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { { xPhys_col(unchosen_idx).setOnes(); - ele_to_write(pixel_idx) = xPhys_col; - Tensor3d ten_xPhys_to_write(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz(), 1, 1); + ele_to_write(global_idx_of_ele_in_use) = xPhys_col; + Tensor3d ten_xPhys_to_write(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz(), 1, + 1); for (int i = 0; i < ele_to_write.size(); ++i) { ten_xPhys_to_write(i, 0, 0) = ele_to_write(i); } diff --git a/src/ThermoelasticTop3d.h b/src/ThermoelasticTop3d.h index 7c51b27..bbfb6c6 100644 --- a/src/ThermoelasticTop3d.h +++ b/src/ThermoelasticTop3d.h @@ -19,7 +19,6 @@ namespace da::sha::top { Eigen::VectorXd GetU() const { return sp_thermal_top3d_->GetU(); -// return sp_mech_top3d_->GetU(); } private: diff --git a/src/Top3d.cpp b/src/Top3d.cpp index cfd3e90..6a55991 100644 --- a/src/Top3d.cpp +++ b/src/Top3d.cpp @@ -35,7 +35,7 @@ namespace da::sha { Eigen::VectorXd ele_to_write = Eigen::VectorXd::Zero(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz()); - Eigen::VectorXi pixel_idx = sp_mesh_->GetPixelIdx(); + Eigen::VectorXi pixel_idx = sp_mesh_->GetGlobalIdxOfEleInUse(); spdlog::info("end precompute"); // // clear output dir // clear_dir(CMAKE_SOURCE_DIR "/output"); @@ -156,7 +156,7 @@ namespace da::sha { std::vector Top3d::GetTensorOfStress(const Eigen::VectorXd &which_col_of_stress) { Eigen::VectorXd ele_to_write = Eigen::VectorXd::Zero(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz()); - Eigen::VectorXi pixel_idx = sp_mesh_->GetPixelIdx(); + Eigen::VectorXi pixel_idx = sp_mesh_->GetGlobalIdxOfEleInUse(); // stress Eigen::MatrixXd mat_stress(sp_mesh_->GetNumEles(), 6); Eigen::MatrixXd B = sp_fea_->computeBe({0, 0, 0}); @@ -204,7 +204,7 @@ namespace da::sha { Eigen::VectorXd sH(iH.size()); sH.setZero(); int cnt = 0; - Hs_ = Eigen::VectorXd(sp_mesh_->GetNumEles()); + Hs_ = Eigen::VectorXd::Zero(sp_mesh_->GetNumEles()); int delta = std::ceil(sp_para_->r_min) - 1; for (int k = 0; k < sp_mesh_->GetLz(); ++k) { diff --git a/src/Util.cpp b/src/Util.cpp index 8d4ea1d..0e72635 100644 --- a/src/Util.cpp +++ b/src/Util.cpp @@ -124,7 +124,7 @@ namespace da::sha { matmesh.mat_hexahedrons = sp_mesh->GetEleId2NodeIdsMap(); std::vector v_ele_data; - Eigen::VectorXi pixel_idx = sp_mesh->GetPixelIdx(); + Eigen::VectorXi pixel_idx = sp_mesh->GetGlobalIdxOfEleInUse(); for (int i = 0; i < pixel_idx.size(); ++i) { v_ele_data.push_back(*(t3.data() + pixel_idx[i])); }