From 1e48f868c1800347412556c53f79ea6f3dcaae52 Mon Sep 17 00:00:00 2001 From: cflin Date: Thu, 15 Jun 2023 11:26:58 +0800 Subject: [PATCH] before debug --- .../CMakeLists.txt | 2 +- .../top-mech-regular-model-test/config.json | 14 +- examples/top-mech-regular-model-test/main.cpp | 112 ++--- .../config.json | 19 +- .../top-thermoelastic-regular-model/main.cpp | 11 +- src/FEA/FEA.h | 20 +- src/FEA/MechanicalLinearFEA.cpp | 58 +++ src/FEA/MechanicalLinearFEA.h | 10 + src/FEA/ThermalLinearFEA.cpp | 108 ++++- src/FEA/ThermalLinearFEA.h | 10 + src/Mesh/IrregularMesh.cpp | 2 +- src/Mesh/Mesh.h | 3 + src/ThermoelasticTop3d.h | 417 +++++++++++++++++- src/Top3d.cpp | 2 +- src/Top3d.h | 29 +- 15 files changed, 700 insertions(+), 117 deletions(-) diff --git a/examples/top-mech-regular-model-test/CMakeLists.txt b/examples/top-mech-regular-model-test/CMakeLists.txt index 3ea4c6f..48b97b7 100644 --- a/examples/top-mech-regular-model-test/CMakeLists.txt +++ b/examples/top-mech-regular-model-test/CMakeLists.txt @@ -1,4 +1,4 @@ -set(SUB_PROJECT_NAME top-mesh-regular-model) +set(SUB_PROJECT_NAME top-mesh-regular-model-test) #file(GLOB CPP_LIST ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp) #add_executable(${SUB_PROJECT_NAME} ${CPP_LIST}) add_executable(${SUB_PROJECT_NAME} main.cpp) diff --git a/examples/top-mech-regular-model-test/config.json b/examples/top-mech-regular-model-test/config.json index e01ec0d..9da75b5 100644 --- a/examples/top-mech-regular-model-test/config.json +++ b/examples/top-mech-regular-model-test/config.json @@ -4,23 +4,23 @@ "poisson_ratio": 0.3 }, "topology": { - "max_loop": 5, + "max_loop": 100, "volfrac": 0.5, - "penal": 1.0, + "penal": 3.0, "r_min": 1.5 }, "model": { "regular_model": { - "lx": 30, + "lx": 2, "ly": 20, - "lz": 20 + "lz": 30 }, "visual_model":"2-box-refined.obj" }, "NBC": [ { - "min": [1, 0, 0], - "max": [1, 1, 1], + "min": [0, 1, 0], + "max": [1, 1, 0], "val": [0.0, 0.0, -1.0] } @@ -28,7 +28,7 @@ "DBC": [ { "min": [0, 0, 0], - "max": [0, 1, 1], + "max": [1, 0, 1], "dir": [1, 1, 1] } ] diff --git a/examples/top-mech-regular-model-test/main.cpp b/examples/top-mech-regular-model-test/main.cpp index 46cf555..d3a63e9 100644 --- a/examples/top-mech-regular-model-test/main.cpp +++ b/examples/top-mech-regular-model-test/main.cpp @@ -4,7 +4,7 @@ #include #include #include "Boundary.h" -#include "IrregularMesh.h" +#include "Mesh/IrregularMesh.h" #include "Top3d.h" #include "Util.h" @@ -39,7 +39,9 @@ int main() { // set material parameters double E = j_config["material"]["E"]; double Poisson_ratio = j_config["material"]["poisson_ratio"]; - auto material = std::make_shared(E, Poisson_ratio,1.0); + auto material = std::make_shared(); + material->E=E; + material->poisson_ratio=Poisson_ratio; // set fea auto sp_fea=std::make_shared(material); @@ -91,59 +93,59 @@ int main() { mesh->GetOrigin() + mesh->GetLenBox()); WriteTensorToVtk(output_dir /"vtk" / "top_mach_regular_field_matrix.vtk", ten_rho, mesh); // extract stress field - auto v_tenosr = top3d.GetTensorOfStress(Eigen::Vector3d{0, 1, 2}); - write_tensor3d(output_dir / "stressX_field_matrix.txt", v_tenosr[0], - mesh->GetOrigin(), mesh->GetOrigin() + mesh->GetLenBox()); - write_tensor3d(output_dir /"stressY_field_matrix.txt", v_tenosr[1], - mesh->GetOrigin(), mesh->GetOrigin() + mesh->GetLenBox()); - write_tensor3d(output_dir /"stressZ_field_matrix.txt", v_tenosr[2], - mesh->GetOrigin(), mesh->GetOrigin() + mesh->GetLenBox()); - WriteTensorToVtk(output_dir / "stressX_field_matrix.vtk", v_tenosr[0], mesh); - WriteTensorToVtk(output_dir / "stressY_field_matrix.vtk", v_tenosr[1], mesh); - WriteTensorToVtk(output_dir / "stressZ_field_matrix.vtk", v_tenosr[2], mesh); - if (!j_config["model"].count("visual_model")) { - // visual model with tensor - // rho - { - auto [mesh_obj, vertex_propty] = top::GetMeshVertexPropty(assets_dir / j_config["model"]["visual_model"], - top3d.GetRhoFieldOneFilled(), bdr, true); - std::string vtk_path = (output_dir / "regular_with_rho.vtk").string(); - top::WriteTriVTK(vtk_path, mesh_obj.mat_coordinates, mesh_obj.mat_faces, {}, - std::vector(vertex_propty.data(), vertex_propty.data() + vertex_propty.size())); - spdlog::info("write vtk with rho to: {}", vtk_path); - top::WriteVectorXd((output_dir / "regular_rho.txt").string(), vertex_propty); - } - // stressX - { - auto [mesh_obj, vertex_propty] = top::GetMeshVertexPropty( - assets_dir / j_config["model"]["visual_model"], v_tenosr[0], bdr, true); - std::string vtk_path = (output_dir / "regular_with_stressX.vtk").string(); - top::WriteTriVTK(vtk_path, mesh_obj.mat_coordinates, mesh_obj.mat_faces, {}, - std::vector(vertex_propty.data(), vertex_propty.data() + vertex_propty.size())); - spdlog::info("write vtk with stressX to: {}", vtk_path); - top::WriteVectorXd((output_dir / "regular_stressX.txt").string(), vertex_propty); - } - // stressY - { - auto [mesh_obj, vertex_propty] = top::GetMeshVertexPropty( - assets_dir / j_config["model"]["visual_model"], v_tenosr[1], bdr, true); - std::string vtk_path = (output_dir / "regular_with_stressY.vtk").string(); - top::WriteTriVTK(vtk_path, mesh_obj.mat_coordinates, mesh_obj.mat_faces, {}, - std::vector(vertex_propty.data(), vertex_propty.data() + vertex_propty.size())); - spdlog::info("write vtk with stressY to: {}", vtk_path); - top::WriteVectorXd((output_dir / "regular_stressY.txt").string(), vertex_propty); - } - // stressZ - { - auto [mesh_obj, vertex_propty] = top::GetMeshVertexPropty( - assets_dir / j_config["model"]["visual_model"], v_tenosr[2], bdr, true); - std::string vtk_path = (output_dir / "regular_with_stressZ.vtk").string(); - top::WriteTriVTK(vtk_path, mesh_obj.mat_coordinates, mesh_obj.mat_faces, {}, - std::vector(vertex_propty.data(), vertex_propty.data() + vertex_propty.size())); - spdlog::info("write vtk with stressZ to: {}", vtk_path); - top::WriteVectorXd((output_dir / "regular_stressZ.txt").string(), vertex_propty); - } - } +// auto v_tenosr = top3d.GetTensorOfStress(Eigen::Vector3d{0, 1, 2}); +// write_tensor3d(output_dir / "stressX_field_matrix.txt", v_tenosr[0], +// mesh->GetOrigin(), mesh->GetOrigin() + mesh->GetLenBox()); +// write_tensor3d(output_dir /"stressY_field_matrix.txt", v_tenosr[1], +// mesh->GetOrigin(), mesh->GetOrigin() + mesh->GetLenBox()); +// write_tensor3d(output_dir /"stressZ_field_matrix.txt", v_tenosr[2], +// mesh->GetOrigin(), mesh->GetOrigin() + mesh->GetLenBox()); +// WriteTensorToVtk(output_dir / "stressX_field_matrix.vtk", v_tenosr[0], mesh); +// WriteTensorToVtk(output_dir / "stressY_field_matrix.vtk", v_tenosr[1], mesh); +// WriteTensorToVtk(output_dir / "stressZ_field_matrix.vtk", v_tenosr[2], mesh); +// if ( 0 && !j_config["model"].count("visual_model")) { +// // visual model with tensor +// // rho +// { +// auto [mesh_obj, vertex_propty] = top::GetMeshVertexPropty(assets_dir / j_config["model"]["visual_model"], +// top3d.GetRhoFieldOneFilled(), bdr, true); +// std::string vtk_path = (output_dir / "regular_with_rho.vtk").string(); +// top::WriteTriVTK(vtk_path, mesh_obj.mat_coordinates, mesh_obj.mat_faces, {}, +// std::vector(vertex_propty.data(), vertex_propty.data() + vertex_propty.size())); +// spdlog::info("write vtk with rho to: {}", vtk_path); +// top::WriteVectorXd((output_dir / "regular_rho.txt").string(), vertex_propty); +// } +// // stressX +// { +// auto [mesh_obj, vertex_propty] = top::GetMeshVertexPropty( +// assets_dir / j_config["model"]["visual_model"], v_tenosr[0], bdr, true); +// std::string vtk_path = (output_dir / "regular_with_stressX.vtk").string(); +// top::WriteTriVTK(vtk_path, mesh_obj.mat_coordinates, mesh_obj.mat_faces, {}, +// std::vector(vertex_propty.data(), vertex_propty.data() + vertex_propty.size())); +// spdlog::info("write vtk with stressX to: {}", vtk_path); +// top::WriteVectorXd((output_dir / "regular_stressX.txt").string(), vertex_propty); +// } +// // stressY +// { +// auto [mesh_obj, vertex_propty] = top::GetMeshVertexPropty( +// assets_dir / j_config["model"]["visual_model"], v_tenosr[1], bdr, true); +// std::string vtk_path = (output_dir / "regular_with_stressY.vtk").string(); +// top::WriteTriVTK(vtk_path, mesh_obj.mat_coordinates, mesh_obj.mat_faces, {}, +// std::vector(vertex_propty.data(), vertex_propty.data() + vertex_propty.size())); +// spdlog::info("write vtk with stressY to: {}", vtk_path); +// top::WriteVectorXd((output_dir / "regular_stressY.txt").string(), vertex_propty); +// } +// // stressZ +// { +// auto [mesh_obj, vertex_propty] = top::GetMeshVertexPropty( +// assets_dir / j_config["model"]["visual_model"], v_tenosr[2], bdr, true); +// std::string vtk_path = (output_dir / "regular_with_stressZ.vtk").string(); +// top::WriteTriVTK(vtk_path, mesh_obj.mat_coordinates, mesh_obj.mat_faces, {}, +// std::vector(vertex_propty.data(), vertex_propty.data() + vertex_propty.size())); +// spdlog::info("write vtk with stressZ to: {}", vtk_path); +// top::WriteVectorXd((output_dir / "regular_stressZ.txt").string(), vertex_propty); +// } +// } } diff --git a/examples/top-thermoelastic-regular-model/config.json b/examples/top-thermoelastic-regular-model/config.json index 1a2e741..e381535 100644 --- a/examples/top-thermoelastic-regular-model/config.json +++ b/examples/top-thermoelastic-regular-model/config.json @@ -2,17 +2,22 @@ "material": { "E": 1.0, "poisson_ratio": 0.3, - "thermal_conductivity":0.1 + "thermal_conductivity":1e11, "unit": "W/K", + "thermal_expansion_coefficient": 1e-15, "unit": "1/K" }, "topology": { "max_loop": 100, "volfrac": 0.5, - "penal": 2.0, - "r_min": 1.5 + "r_min": 1.5, + "T_ref": 295, + "T_limit": 3250, + "R_E": 8, + "R_lambda": 0, + "R_beta": 0 }, "model": { "regular_model": { - "lx": 2, + "lx": 1, "ly": 40, "lz": 30 }, @@ -23,7 +28,7 @@ { "min": [0, 1, 0], "max": [1, 1, 0], - "val": [0.0, 0.0, -1.0] + "val": [0.0, 0.0, -500.0] } ], @@ -40,12 +45,12 @@ { "min": [0, 0.25, 0.5], "max": [1, 0.25, 0.5], - "heat_flux": 0.06, "unit": "W/m^2" + "heat_flux": 0.06, "unit": "W" }, { "min": [0, 0.75, 0.5], "max": [1, 0.75, 0.5], - "heat_flux": 0.06, "unit": "W/m^2" + "heat_flux": 0.06, "unit": "W" } ], diff --git a/examples/top-thermoelastic-regular-model/main.cpp b/examples/top-thermoelastic-regular-model/main.cpp index 0ef1a06..ae6ed05 100644 --- a/examples/top-thermoelastic-regular-model/main.cpp +++ b/examples/top-thermoelastic-regular-model/main.cpp @@ -33,14 +33,19 @@ int main() { auto para = std::make_shared(); para->max_loop = j_config["topology"]["max_loop"]; para->volfrac = j_config["topology"]["volfrac"]; - para->penal = j_config["topology"]["penal"]; para->r_min = j_config["topology"]["r_min"]; + 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"]; - auto material = std::make_shared(E, Poisson_ratio, 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 @@ -128,6 +133,6 @@ int main() { 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); +// WriteUToVtk(output_dir / "vtk" / "top_thermoelastic_regular_U.vtk",ther_top3d.GetU(),sp_mech_mesh); } diff --git a/src/FEA/FEA.h b/src/FEA/FEA.h index a691cd4..c584959 100644 --- a/src/FEA/FEA.h +++ b/src/FEA/FEA.h @@ -7,6 +7,7 @@ #include #include "Eigen/Eigen" + namespace da { namespace sha { namespace top { @@ -14,11 +15,15 @@ namespace da { double E = 1.0; // Young's modulus double poisson_ratio = 0.3; double thermal_conductivity = 1; - - Material(double E, double poisson_ratio, double thermal_conductivity) : E(E), - poisson_ratio(poisson_ratio), - thermal_conductivity( - thermal_conductivity) {} + double thermal_expansion_coefficient = 1e-5; + Material()=default; + Material(double E, double poisson_ratio, double thermal_conductivity, + double thermal_expansion_coefficient) : E(E), + poisson_ratio(poisson_ratio), + thermal_conductivity( + thermal_conductivity), + thermal_expansion_coefficient( + thermal_expansion_coefficient) {} }; class FEA { @@ -27,9 +32,14 @@ namespace da { virtual Eigen::MatrixXd computeKe() = 0; + virtual Eigen::MatrixXd computeKe(double stiffness_coef)=0; + + virtual Eigen::MatrixXd computeBe() = 0; + virtual Eigen::MatrixXd computeBe(const Eigen::Vector3d &P) = 0; virtual Eigen::MatrixXd computeN(const Eigen::Vector3d &P) = 0; + virtual Eigen::MatrixXd computeD(double stiffness_coef) = 0; virtual Eigen::MatrixXd computeD() = 0; diff --git a/src/FEA/MechanicalLinearFEA.cpp b/src/FEA/MechanicalLinearFEA.cpp index d6012e7..d46bb0f 100644 --- a/src/FEA/MechanicalLinearFEA.cpp +++ b/src/FEA/MechanicalLinearFEA.cpp @@ -94,6 +94,64 @@ namespace da { return Ke; } + Eigen::MatrixXd MechanicalLinearFEA::computeBe() { + Eigen::MatrixXd B_inted = Eigen::MatrixXd::Zero(6, 24); + Eigen::Vector2d GP(-1.0 / sqrt(3.0), 1.0 / sqrt(3.0)); + Eigen::Vector2d GW(1.0, 1.0); + Eigen::Matrix tmp{ + {-a, -b, -c}, + {a, -b, -c}, + {a, b, -c}, + {-a, b, -c}, + {-a, -b, c}, + {a, -b, c}, + {a, b, c}, + {-a, b, c} + }; + + for (int i = 0; i < 2; ++i) { + for (int j = 0; j < 2; ++j) { + for (int k = 0; k < 2; ++k) { + double x = GP(i), y = GP(j), z = GP(k); + Eigen::RowVector dNx{-(1.0 - y) * (1.0 - z) / 8.0, + (1.0 - y) * (1.0 - z) / 8.0, + (1.0 + y) * (1.0 - z) / 8.0, + -(1.0 + y) * (1.0 - z) / 8.0, + -(1.0 - y) * (1.0 + z) / 8.0, + (1.0 - y) * (1.0 + z) / 8.0, + (1.0 + y) * (1.0 + z) / 8.0, + -(1.0 + y) * (1.0 + z) / 8.0}; + Eigen::RowVector dNy{-(1.0 - x) * (1.0 - z) / 8.0, + -(1.0 + x) * (1.0 - z) / 8.0, + (1.0 + x) * (1.0 - z) / 8.0, + (1.0 - x) * (1.0 - z) / 8.0, + -(1.0 - x) * (1.0 + z) / 8.0, + -(1.0 + x) * (1.0 + z) / 8.0, + (1.0 + x) * (1.0 + z) / 8.0, + (1.0 - x) * (1.0 + z) / 8.0}; + Eigen::RowVector dNz{-(1.0 - x) * (1.0 - y) / 8.0, + -(1.0 + x) * (1.0 - y) / 8.0, + -(1.0 + x) * (1.0 + y) / 8.0, + -(1.0 - x) * (1.0 + y) / 8.0, + (1.0 - x) * (1.0 - y) / 8.0, + (1.0 + x) * (1.0 - y) / 8.0, + (1.0 + x) * (1.0 + y) / 8.0, + (1.0 - x) * (1.0 + y) / 8.0}; + + Eigen::Matrix tmp1; + tmp1(0, Eigen::all) = dNx; + tmp1(1, Eigen::all) = dNy; + tmp1(2, Eigen::all) = dNz; + Eigen::Matrix3d J = tmp1 * tmp; + Eigen::Matrix3d JInv = J.inverse(); + double JDet = J.determinant(); + B_inted = B_inted + GW(i) * GW(j) * GW(k) * JDet * computeBe({x,y,z}); + } + } + } + return B_inted; + } + Eigen::MatrixXd MechanicalLinearFEA::computeBe(const Eigen::Vector3d &P) { Eigen::MatrixXd B = Eigen::MatrixXd::Zero(6, 24); double x = P(0), y = P(1), z = P(2); diff --git a/src/FEA/MechanicalLinearFEA.h b/src/FEA/MechanicalLinearFEA.h index 872a850..cbb9776 100644 --- a/src/FEA/MechanicalLinearFEA.h +++ b/src/FEA/MechanicalLinearFEA.h @@ -15,12 +15,22 @@ namespace da { public: MechanicalLinearFEA(std::shared_ptr sp_material) : FEA(sp_material) {} + Eigen::MatrixXd computeKe(double stiffness_coef) { + return computeKe() / sp_material_->E * stiffness_coef; + } + Eigen::MatrixXd computeKe(); + Eigen::MatrixXd computeBe(); + Eigen::MatrixXd computeBe(const Eigen::Vector3d &P); Eigen::MatrixXd computeN(const Eigen::Vector3d &P); + Eigen::MatrixXd computeD(double stiffness_coef) { + return computeD() / sp_material_->E * stiffness_coef; + }; + Eigen::MatrixXd computeD(); }; diff --git a/src/FEA/ThermalLinearFEA.cpp b/src/FEA/ThermalLinearFEA.cpp index abde62b..2e1cdf1 100644 --- a/src/FEA/ThermalLinearFEA.cpp +++ b/src/FEA/ThermalLinearFEA.cpp @@ -69,34 +69,94 @@ namespace da { return heatKe; } + Eigen::MatrixXd ThermalLinearFEA::computeBe() { + Eigen::Matrix Be_inted; + Be_inted.setZero(); + Eigen::Vector2d GP(-1.0 / sqrt(3.0), 1.0 / sqrt(3.0)); + Eigen::Vector2d GW(1.0, 1.0); + + Eigen::Matrix tmp{ + {-a, -b, -c}, + {a, -b, -c}, + {a, b, -c}, + {-a, b, -c}, + {-a, -b, c}, + {a, -b, c}, + {a, b, c}, + {-a, b, c} + }; + + for (int i = 0; i < 2; ++i) { + for (int j = 0; j < 2; ++j) { + for (int k = 0; k < 2; ++k) { + double x = GP(i), y = GP(j), z = GP(k); + Eigen::RowVector dNx{-(1.0 - y) * (1.0 - z) / 8.0, + (1.0 - y) * (1.0 - z) / 8.0, + (1.0 + y) * (1.0 - z) / 8.0, + -(1.0 + y) * (1.0 - z) / 8.0, + -(1.0 - y) * (1.0 + z) / 8.0, + (1.0 - y) * (1.0 + z) / 8.0, + (1.0 + y) * (1.0 + z) / 8.0, + -(1.0 + y) * (1.0 + z) / 8.0}; + Eigen::RowVector dNy{-(1.0 - x) * (1.0 - z) / 8.0, + -(1.0 + x) * (1.0 - z) / 8.0, + (1.0 + x) * (1.0 - z) / 8.0, + (1.0 - x) * (1.0 - z) / 8.0, + -(1.0 - x) * (1.0 + z) / 8.0, + -(1.0 + x) * (1.0 + z) / 8.0, + (1.0 + x) * (1.0 + z) / 8.0, + (1.0 - x) * (1.0 + z) / 8.0}; + Eigen::RowVector dNz{-(1.0 - x) * (1.0 - y) / 8.0, + -(1.0 + x) * (1.0 - y) / 8.0, + -(1.0 + x) * (1.0 + y) / 8.0, + -(1.0 - x) * (1.0 + y) / 8.0, + (1.0 - x) * (1.0 - y) / 8.0, + (1.0 + x) * (1.0 - y) / 8.0, + (1.0 + x) * (1.0 + y) / 8.0, + (1.0 - x) * (1.0 + y) / 8.0}; + + Eigen::Matrix tmp1; + tmp1(0, Eigen::all) = dNx; + tmp1(1, Eigen::all) = dNy; + tmp1(2, Eigen::all) = dNz; + Eigen::Matrix3d J = tmp1 * tmp; + Eigen::Matrix3d JInv = J.inverse(); + double JDet = J.determinant(); + Be_inted = Be_inted + GW(i) * GW(j) * GW(k) * JDet * computeBe({x, y, z}); + } + } + } + return Be_inted; + } + Eigen::MatrixXd ThermalLinearFEA::computeBe(const Eigen::Vector3d &P) { double x = P(0), y = P(1), z = P(2); assert(x >= -1.0 && x <= 1.0 && y >= -1.0 && y <= 1.0 && z >= -1.0 && z <= 1.0); Eigen::Matrix dNi; - dNi(0,0) = 0.125 * - 1 * (1.0 - y) * (1.0 - z); - dNi(0,1) = 0.125 * 1 * (1.0 - y) * (1.0 - z); - dNi(0,2) = 0.125 * 1 * (1.0 - y) * (1.0 - z); - dNi(0,3) = 0.125 * - 1 * (1.0 + y) * (1.0 - z); - dNi(0,4) = 0.125 * - 1 * (1.0 + y) * (1.0 - z); - dNi(0,5) = 0.125 * 1 * (1.0 - y) * (1.0 + z); - dNi(0,6) = 0.125 * 1 * (1.0 - y) * (1.0 + z); - dNi(0,7) = 0.125 * - 1 * (1.0 + y) * (1.0 + z); - dNi(1,0) = 0.125 * (1.0 - x) * - 1 * (1.0 - z); - dNi(1,1) = 0.125 * (1.0 + x) * - 1 * (1.0 - z); - dNi(1,2) = 0.125 * (1.0 + x) * - 1 * (1.0 - z); - dNi(1,3) = 0.125 * (1.0 - x) * + 1 * (1.0 - z); - dNi(1,4) = 0.125 * (1.0 - x) * + 1 * (1.0 - z); - dNi(1,5) = 0.125 * (1.0 + x) * - 1 * (1.0 + z); - dNi(1,6) = 0.125 * (1.0 + x) * - 1 * (1.0 + z); - dNi(1,7) = 0.125 * (1.0 - x) * + 1 * (1.0 + z); - dNi(2,0) = 0.125 * (1.0 - x) * (1.0 - y) * - 1; - dNi(2,1) = 0.125 * (1.0 + x) * (1.0 - y) * - 1; - dNi(2,2) = 0.125 * (1.0 + x) * (1.0 + y) * - 1; - dNi(2,3) = 0.125 * (1.0 - x) * (1.0 + y) * - 1; - dNi(2,4) = 0.125 * (1.0 - x) * (1.0 - y) * + 1; - dNi(2,5) = 0.125 * (1.0 + x) * (1.0 - y) * + 1; - dNi(2,6) = 0.125 * (1.0 + x) * (1.0 + y) * + 1; - dNi(2,7) = 0.125 * (1.0 - x) * (1.0 + y) * + 1; + dNi(0, 0) = 0.125 * -1 * (1.0 - y) * (1.0 - z); + dNi(0, 1) = 0.125 * 1 * (1.0 - y) * (1.0 - z); + dNi(0, 2) = 0.125 * 1 * (1.0 - y) * (1.0 - z); + dNi(0, 3) = 0.125 * -1 * (1.0 + y) * (1.0 - z); + dNi(0, 4) = 0.125 * -1 * (1.0 + y) * (1.0 - z); + dNi(0, 5) = 0.125 * 1 * (1.0 - y) * (1.0 + z); + dNi(0, 6) = 0.125 * 1 * (1.0 - y) * (1.0 + z); + dNi(0, 7) = 0.125 * -1 * (1.0 + y) * (1.0 + z); + dNi(1, 0) = 0.125 * (1.0 - x) * -1 * (1.0 - z); + dNi(1, 1) = 0.125 * (1.0 + x) * -1 * (1.0 - z); + dNi(1, 2) = 0.125 * (1.0 + x) * -1 * (1.0 - z); + dNi(1, 3) = 0.125 * (1.0 - x) * +1 * (1.0 - z); + dNi(1, 4) = 0.125 * (1.0 - x) * +1 * (1.0 - z); + dNi(1, 5) = 0.125 * (1.0 + x) * -1 * (1.0 + z); + dNi(1, 6) = 0.125 * (1.0 + x) * -1 * (1.0 + z); + dNi(1, 7) = 0.125 * (1.0 - x) * +1 * (1.0 + z); + dNi(2, 0) = 0.125 * (1.0 - x) * (1.0 - y) * -1; + dNi(2, 1) = 0.125 * (1.0 + x) * (1.0 - y) * -1; + dNi(2, 2) = 0.125 * (1.0 + x) * (1.0 + y) * -1; + dNi(2, 3) = 0.125 * (1.0 - x) * (1.0 + y) * -1; + dNi(2, 4) = 0.125 * (1.0 - x) * (1.0 - y) * +1; + dNi(2, 5) = 0.125 * (1.0 + x) * (1.0 - y) * +1; + dNi(2, 6) = 0.125 * (1.0 + x) * (1.0 + y) * +1; + dNi(2, 7) = 0.125 * (1.0 - x) * (1.0 + y) * +1; return dNi; } diff --git a/src/FEA/ThermalLinearFEA.h b/src/FEA/ThermalLinearFEA.h index 5991296..0dfe8a6 100644 --- a/src/FEA/ThermalLinearFEA.h +++ b/src/FEA/ThermalLinearFEA.h @@ -15,12 +15,22 @@ namespace da { public: ThermalLinearFEA(std::shared_ptr sp_material) : FEA(sp_material) {} + Eigen::MatrixXd computeKe(double stiffness_coef) { + return computeKe() / sp_material_->thermal_conductivity * stiffness_coef; + } + Eigen::MatrixXd computeKe(); + Eigen::MatrixXd computeBe(); + Eigen::MatrixXd computeBe(const Eigen::Vector3d &P); Eigen::MatrixXd computeN(const Eigen::Vector3d &P); + Eigen::MatrixXd computeD(double stiffness_coef) { + return computeD() / sp_material_->thermal_conductivity * stiffness_coef; + }; + Eigen::MatrixXd computeD() { return (Eigen::MatrixXd(1, 1) << sp_material_->thermal_conductivity).finished(); } diff --git a/src/Mesh/IrregularMesh.cpp b/src/Mesh/IrregularMesh.cpp index 6661e20..f82caf1 100644 --- a/src/Mesh/IrregularMesh.cpp +++ b/src/Mesh/IrregularMesh.cpp @@ -57,7 +57,7 @@ namespace da::sha { 1, 1, 1, 0, 1, 1 ).finished(); - // precompute coords needed to query + // Precompute coords needed to query const int SAMPLE_POINTS_EACH_DIM=3; double space=1.0/SAMPLE_POINTS_EACH_DIM; diff --git a/src/Mesh/Mesh.h b/src/Mesh/Mesh.h index 2aaa69d..27f3e1c 100644 --- a/src/Mesh/Mesh.h +++ b/src/Mesh/Mesh.h @@ -80,6 +80,9 @@ namespace da::sha { int Get_DOFS_EACH_NODE()const{ return DOFS_EACH_NODE; } + int Get_DOFS_EACH_ELE()const{ + return NUM_NODES_EACH_ELE*DOFS_EACH_NODE; + } Eigen::VectorXd GetInitEleRho() const{ return init_ele_rho_; } diff --git a/src/ThermoelasticTop3d.h b/src/ThermoelasticTop3d.h index 1df5cc4..5bee205 100644 --- a/src/ThermoelasticTop3d.h +++ b/src/ThermoelasticTop3d.h @@ -12,19 +12,380 @@ namespace da::sha::top { public: ThermoelasticTop3d(std::shared_ptr sp_mech_top3d, std::shared_ptr sp_thermal_top3d) : sp_mech_top3d_(sp_mech_top3d), sp_thermal_top3d_(sp_thermal_top3d) { + Precompute(); } - void AddThermalDBC(const Eigen::MatrixXi &thermal_DBC_coords, double temperature) { - sp_thermal_top3d_->AddDBC(thermal_DBC_coords, temperature); - } + Tensor3d TopOptMainLoop() { + auto &sp_mesh_ = sp_mech_top3d_->sp_mesh_; + auto &sp_para_ = sp_mech_top3d_->sp_para_; + Eigen::VectorXd xPhys_col(sp_mesh_->GetNumEles()); + xPhys_col.setConstant(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); +// } - void AddThermalNBC(const Eigen::MatrixXi &thermal_NBC_coords, double heat_flux) { - sp_thermal_top3d_->AddNBC(thermal_NBC_coords, Eigen::Vector3d(heat_flux, 0, 0)); - } + int loop = 0; + double change = 1.0; + double E0_m = sp_mech_top3d_->sp_fea_->sp_material_->E; + double lambda0 = sp_mech_top3d_->sp_fea_->sp_material_->thermal_conductivity; + double lambda_min = lambda0 * sp_mech_top3d_->sp_para_->E_factor; + double alpha0 = sp_mech_top3d_->sp_fea_->sp_material_->thermal_expansion_coefficient; + // Precompute - Tensor3d TopOptMainLoop() { - return sp_thermal_top3d_->TopOptMainLoop(); -// return sp_mech_top3d_->TopOptMainLoop(); + Eigen::VectorXd dv(sp_mesh_->GetNumEles()); + dv.setOnes(); + dv = sp_mech_top3d_->H_ * (dv.array() / sp_mech_top3d_->Hs_.array()).matrix().eval(); + + Eigen::VectorXd ele_to_write = + Eigen::VectorXd::Zero(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz()); + Eigen::VectorXi pixel_idx = sp_mesh_->GetPixelIdx(); + + // dofs of limited T + struct LimitedDof { + int dof; + 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) {} + }; + std::map> map_ele2Limit; + std::vector v_dof(sp_thermal_top3d_->set_dofs_to_load.begin(), + sp_thermal_top3d_->set_dofs_to_load.end()); + { + Eigen::MatrixXi ele2dof_map = sp_thermal_top3d_->sp_mesh_->GetEleId2DofsMap(); + // loop ele2dof_map + for (int i = 0; i < ele2dof_map.rows(); ++i) { + for (int j = 0; j < ele2dof_map.cols(); ++j) { + for (int k = 0; k < v_dof.size(); ++k) { + if (ele2dof_map(i, j) == v_dof[k]) { + if (map_ele2Limit.find(i) == map_ele2Limit.end()) { + map_ele2Limit[i] = {LimitedDof(v_dof[k], k, j)}; + } else { + map_ele2Limit[i].push_back(LimitedDof(v_dof[k], k, j)); + } + } + } + } + } + } + spdlog::info("end Precompute"); + +#ifdef USE_SUITESPARSE + spdlog::info("using suitesparse solver"); +#else + spdlog::warn("using Eigen built-in direct solver!"); +#endif + // start iteration + while (change > sp_para_->tol_x*1 && loop < sp_para_->max_loop) { + ++loop; + // filter + xPhys_col = sp_mech_top3d_->H_ * (xPhys_col.array() / sp_mech_top3d_->Hs_.array()).matrix().eval(); + auto CalR = [](double rho, double R) { + return rho / (1.0 + R * (1.0 - rho)); + }; + auto CalR_Vec=[](const Eigen::VectorXd& vec_rho,double R)->Eigen::VectorXd { + return vec_rho.array()/(1.0+R*(1.0-vec_rho.array())); + }; + auto CalDRDrho = [](double rho, double R) { + double down = 1 + R * (1 - rho); + return (1 + R) / down * down; + }; + auto CalDRDrho_Vec=[](const Eigen::VectorXd&vec_rho,double R)->Eigen::VectorXd { + auto down=1+R*(1-vec_rho.array()); + return (1+R)/down.pow(2); + }; + auto CalE_Vec = [&](const Eigen::VectorXd& vec_rho)->Eigen::VectorXd { + return CalR_Vec(vec_rho, sp_para_->R_E) * E0_m; + }; + auto CalDEDrho_Vec=[&](const Eigen::VectorXd& vec_rho)->Eigen::VectorXd { + return CalDRDrho_Vec(vec_rho,sp_para_->R_E) * E0_m; + }; + 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); + }; + auto CalDlambdaDrho = [&](double rho) { + return CalDRDrho(rho, sp_para_->R_lambda) * (lambda0 - lambda_min); + }; + auto CalDlambdaDrho_Vec=[&](const Eigen::VectorXd& vec_rho)->Eigen::VectorXd { + return CalDRDrho_Vec(vec_rho,sp_para_->R_lambda) *(lambda0-lambda_min); + }; + + auto CalBeta = [&](double rho) { + return CalR(rho, sp_para_->R_beta) * E0_m * alpha0; + }; + auto CalDBetaDrho = [&](double rho) { + return CalDRDrho(rho, sp_para_->R_beta) * E0_m * alpha0; + }; + // solve T + + Eigen::VectorXd sK_th = + (sp_thermal_top3d_->sKe_ * (lambda_min + + xPhys_col.array() / + (1.0 + sp_para_->R_lambda * (1.0 - xPhys_col.array())) * + (lambda0 - lambda_min)).matrix().transpose()) + .reshaped(); + auto v_tri_th = Vec2Triplet(sp_thermal_top3d_->iK_, sp_thermal_top3d_->jK_, sK_th); + sp_thermal_top3d_->K_.setFromTriplets(v_tri_th.begin(), v_tri_th.end()); + sp_thermal_top3d_->IntroduceFixedDofs(sp_thermal_top3d_->K_, sp_thermal_top3d_->F_); +#ifdef USE_SUITESPARSE + Eigen::CholmodSupernodalLLT> solver_th; +#else + Eigen::SimplicialLLT> solver_th; +#endif + solver_th.compute(sp_thermal_top3d_->K_); + sp_thermal_top3d_->U_ = solver_th.solve(sp_thermal_top3d_->F_); + + + // solve U + + Eigen::VectorXd sK_m = + (sp_mech_top3d_->sKe_ * + ( + xPhys_col.array() / + (1.0 + sp_para_->R_E * (1.0 - xPhys_col.array())) * + E0_m).matrix() + .transpose()) + .reshaped(); + auto v_tri_m = Vec2Triplet(sp_mech_top3d_->iK_, sp_mech_top3d_->jK_, sK_m); + sp_mech_top3d_->K_.setFromTriplets(v_tri_m.begin(), v_tri_m.end()); + // for each element + Eigen::VectorXd &T = sp_thermal_top3d_->U_; + Eigen::VectorXd F_th = Eigen::VectorXd::Zero(sp_mech_top3d_->sp_mesh_->GetNumDofs()); + + for (int i = 0; i < sp_thermal_top3d_->sp_mesh_->GetNumEles(); ++i) { + Eigen::VectorXi dofs_th = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs(i); + Eigen::VectorXi dofs_m = sp_mech_top3d_->sp_mesh_->MapEleId2Dofs(i); + double Te = T(dofs_th).mean(); + double beta_rho = CalBeta(xPhys_col(i)); + F_th(dofs_m) += beta_rho * (Te - sp_mech_top3d_->sp_para_->T_ref) * Inted_; + } + Eigen::VectorXd F = Eigen::VectorXd(sp_mech_top3d_->F_) + F_th; + sp_mech_top3d_->IntroduceFixedDofs(sp_mech_top3d_->K_, F); +#ifdef USE_SUITESPARSE + Eigen::CholmodSupernodalLLT> solver; +#else + Eigen::SimplicialLLT> solver; +#endif + solver.compute(sp_mech_top3d_->K_); + sp_mech_top3d_->U_ = solver.solve(F); + + + // compliance + Eigen::VectorXd ce(sp_mesh_->GetNumEles()); + for (int i = 0; i < sp_mesh_->GetNumEles(); ++i) { + Eigen::VectorXi dofs_in_ele_i = sp_mesh_->MapEleId2Dofs(i); + Eigen::VectorXd Ue = sp_mech_top3d_->U_(dofs_in_ele_i); + ce(i) = Ue.transpose() * sp_mech_top3d_->Ke_ * Ue; + } + double c = + ce.transpose() * (xPhys_col.array() / + (1.0 + sp_para_->R_E * (1.0 - xPhys_col.array())) * + E0_m).matrix(); + double v = flg_chosen ? xPhys_col(chosen_ele_id).sum() : xPhys_col.sum(); + + // sensitivity + // lambda_m + Eigen::VectorXd lambda_m = -sp_mech_top3d_->U_; + // dFth_drho + Eigen::SparseMatrix dFth_drho(sp_mech_top3d_->sp_mesh_->GetNumEles(), + sp_mech_top3d_->sp_mesh_->GetNumDofs()); + Eigen::VectorXd v_dFth_drho(i_dFth_drho_.size()); + for (int i = 0; i < sp_thermal_top3d_->sp_mesh_->GetNumEles(); ++i) { + Eigen::VectorXi dofs_th = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs(i); + Eigen::VectorXi dofs_m = sp_mech_top3d_->sp_mesh_->MapEleId2Dofs(i); + double Te = T(dofs_th).mean(); +// 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 + assert(ele_dFth_drho.size() == 24); + v_dFth_drho(Eigen::seqN(i * ele_dFth_drho.rows(), 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()); + + // dFth_dT + Eigen::SparseMatrix dFth_dT(sp_thermal_top3d_->sp_mesh_->GetNumDofs(), + sp_mech_top3d_->sp_mesh_->GetNumDofs()); + Eigen::VectorXd v_dFth_dT(i_dFth_dT_.size()); + for (int i = 0; i < sp_thermal_top3d_->sp_mesh_->GetNumEles(); ++i) { + Eigen::VectorXi dofs_th = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs(i); + Eigen::VectorXi dofs_m = sp_mech_top3d_->sp_mesh_->MapEleId2Dofs(i); + double beta_rho = CalBeta(xPhys_col(i)); +// F_th(dofs_m) += beta_rho * (Te - sp_mech_top3d_->sp_para_->T_ref) * Inted_; + Eigen::MatrixXd ele_dFth_dT = + Eigen::VectorXd::Ones(dofs_th.size()) * 1.0 / 8.0 * beta_rho * Inted_.transpose(); + assert(ele_dFth_dT.rows() == 8 && ele_dFth_dT.cols() == 24); + v_dFth_dT(Eigen::seqN(i * ele_dFth_dT.rows(), 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()); + + Eigen::VectorXd rhs = dFth_dT * lambda_m; + for (auto dof_value: sp_thermal_top3d_->v_dofs_to_set) { + auto [dof, value] = dof_value; + rhs(dof) = sp_thermal_top3d_->K_.coeffRef(dof, dof) * value; + } + // lambda_t + Eigen::VectorXd lambda_t = solver_th.solve(rhs); + // dF_drho + Eigen::SparseMatrix &dF_drho = dFth_drho; + // lambda_m_Mul_dKm_drho_Mul_U + Eigen::VectorXd lambda_m_Mul_dKm_drho_Mul_U = + -((1 + sp_para_->R_E) / (1.0 + sp_para_->R_E * (1.0 - xPhys_col.array())).pow(2)) * ce.array(); + // lambda_t_Mul_dKt_drho_Mul_T + Eigen::VectorXd ce_th(sp_thermal_top3d_->sp_mesh_->GetNumEles()); + for (int i = 0; i < sp_thermal_top3d_->sp_mesh_->GetNumEles(); ++i) { + Eigen::VectorXi dofs_in_ele_i = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs(i); + Eigen::VectorXd Te = sp_thermal_top3d_->U_(dofs_in_ele_i); + 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 = ((1 + (lambda0 - lambda_min)) / + (1.0 + + (lambda0 - lambda_min) * (1.0 - xPhys_col.array())).pow( + 2)) * + ce_th.array(); + // dc_drho +// Eigen::VectorXd dc_drho = lambda_t_Mul_dKt_drho_Mul_T + +// lambda_m_Mul_dKm_drho_Mul_U + +// 2 * Eigen::VectorXd(dF_drho * sp_mech_top3d_->U_); + Eigen::VectorXd dc_drho = + lambda_m_Mul_dKm_drho_Mul_U + + 2 * Eigen::VectorXd(dF_drho * sp_mech_top3d_->U_); + // dT_drho + Eigen::MatrixXd dT_drho = Eigen::MatrixXd::Zero(sp_thermal_top3d_->sp_mesh_->GetNumEles(), + sp_thermal_top3d_->set_dofs_to_load.size()); + for(auto it=map_ele2Limit.begin();it!=map_ele2Limit.end();++it){ + auto [ele_id,v_limited]=*it; + Eigen::VectorXi dofs_in_ele_i = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs(ele_id); + Eigen::VectorXd dKe_th_Mul_T = + CalDlambdaDrho(xPhys_col(ele_id)) * sp_thermal_top3d_->Ke_ * T(dofs_in_ele_i); + + Eigen::MatrixXd Ke_th(sp_thermal_top3d_->sp_mesh_->Get_DOFS_EACH_ELE(), + sp_thermal_top3d_->sp_mesh_->Get_DOFS_EACH_ELE()); + for (int j1 = 0; j1 < Ke_th.cols(); ++j1) { + for (int i1 = 0; i1 < Ke_th.rows(); ++i1) { + Ke_th(i1, j1) = sp_thermal_top3d_->K_.coeffRef(i1, j1); + } + } + Eigen::VectorXd ele_dT_drho = Ke_th.llt().solve(-dKe_th_Mul_T); + for(auto &limited:v_limited){ + dT_drho(ele_id, limited.idx_of_load_dof) = ele_dT_drho(limited.idx_in_ele); + } + } +// for (int i = 0; i < sp_thermal_top3d_->sp_mesh_->GetNumEles(); ++i) { +// Eigen::VectorXi dofs_in_ele_i = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs(i); +// Eigen::VectorXd dKe_th_Mul_T = +// CalDlambdaDrho(xPhys_col(i)) * sp_thermal_top3d_->Ke_ * T(dofs_in_ele_i); +// Eigen::MatrixXd Ke_th(sp_thermal_top3d_->sp_mesh_->Get_DOFS_EACH_ELE(), +// sp_thermal_top3d_->sp_mesh_->Get_DOFS_EACH_ELE()); +// for (int j1 = 0; j1 < Ke_th.cols(); ++j1) { +// for (int i1 = 0; i1 < Ke_th.rows(); ++i1) { +// Ke_th(i1, j1) = sp_thermal_top3d_->K_.coeffRef(i1, j1); +// } +// } +// Eigen::VectorXd ele_dT_drho = Ke_th.llt().solve(-dKe_th_Mul_T); +// dT_drho(i, dofs_in_ele_i) = ele_dT_drho.transpose(); +// } +// for (auto dof_value: sp_thermal_top3d_->v_dofs_to_set) { +// auto [dof, value] = dof_value; +// dT_drho.col(dof).setZero(); +// } + + // dc_dx + Eigen::VectorXd dc_dx = drho_dx_ * dc_drho; + // dT_dx + Eigen::MatrixXd dT_dx = drho_dx_ * dT_drho; + + // mma solver + size_t num_constraints = + 1 + dT_dx.cols();// volume and temperature constraints + size_t num_variables = flg_chosen ? chosen_ele_id.size() : sp_mesh_->GetNumEles(); + auto mma = std::make_shared(num_variables, num_constraints); + 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 / dc_dx.cwiseAbs().maxCoeff(); +// double fval = v - num_variables * sp_para_->volfrac; + Eigen::VectorXd fval = (Eigen::VectorXd(num_constraints) << (v / num_variables - sp_para_->volfrac), + T(v_dof).array() / sp_para_->T_limit - 1).finished(); +// Eigen::VectorXd dfdx = flg_chosen ? dv(chosen_ele_id) : dv; + Eigen::MatrixXd dfdx = (Eigen::MatrixXd(num_variables, num_constraints) + << 1.0 / num_variables * dv, 1.0 / sp_para_->T_limit * dT_dx).finished().transpose(); + static Eigen::VectorXd low_bounds = Eigen::VectorXd::Zero(num_variables); + static Eigen::VectorXd up_bounds = Eigen::VectorXd::Ones(num_variables); + + // spdlog::info("mma update"); + mma->Update(variables_tmp.data(), df0dx.data(), fval.data(), dfdx.data(), low_bounds.data(), + up_bounds.data()); + if (flg_chosen) { + change = (variables_tmp - xPhys_col(chosen_ele_id)).cwiseAbs().maxCoeff(); + xPhys_col(chosen_ele_id) = variables_tmp; + } else { + change = (variables_tmp - xPhys_col).cwiseAbs().maxCoeff(); + xPhys_col = variables_tmp; + } + + spdlog::critical("Iter: {:3d}, Comp: {:.3e}, Vol: {:.2f}, Change: {:f}", loop, c, v, change); + std::cout<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); + } + ten_xPhys_to_write = ten_xPhys_to_write.reshape(Eigen::array{ + sp_mesh_->GetLx(), sp_mesh_->GetLy(), sp_mesh_->GetLz()}); + top::WriteTensorToVtk( + da::WorkingResultDirectoryPath() / ("field_matrix" + std::to_string(loop) + ".vtk"), + ten_xPhys_to_write, sp_mesh_); +#endif + } + // result + sp_mech_top3d_->rho_ = xPhys_col; + // 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(); + { + 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); + for (int i = 0; i < ele_to_write.size(); ++i) { + ten_xPhys_to_write(i, 0, 0) = ele_to_write(i); + } + ten_xPhys_to_write = ten_xPhys_to_write.reshape(Eigen::array{ + sp_mesh_->GetLx(), sp_mesh_->GetLy(), sp_mesh_->GetLz()}); + sp_mech_top3d_->rho_field_zero_filled_ = ten_xPhys_to_write; + } + + { + 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); + for (int i = 0; i < ele_to_write.size(); ++i) { + ten_xPhys_to_write(i, 0, 0) = ele_to_write(i); + } + ten_xPhys_to_write = ten_xPhys_to_write.reshape(Eigen::array{ + sp_mesh_->GetLx(), sp_mesh_->GetLy(), sp_mesh_->GetLz()}); + sp_mech_top3d_->rho_field_one_filled_ = ten_xPhys_to_write; + } + + return sp_mech_top3d_->rho_field_zero_filled_; } Eigen::VectorXd GetU() const { @@ -32,8 +393,46 @@ namespace da::sha::top { // return sp_mech_top3d_->GetU(); } + private: + void Precompute() { + i_dFth_dT_ = Eigen::KroneckerProduct( + sp_thermal_top3d_->sp_mesh_->GetEleId2DofsMap(), + Eigen::VectorXi::Ones(sp_mech_top3d_->sp_mesh_->Get_DOFS_EACH_ELE())) + .transpose() + .reshaped(); + j_dFth_dT_ = Eigen::KroneckerProduct(sp_mech_top3d_->sp_mesh_->GetEleId2DofsMap(), + Eigen::RowVectorXi::Ones( + sp_thermal_top3d_->sp_mesh_->Get_DOFS_EACH_ELE())) + .transpose() + .reshaped(); + + i_dFth_drho_ = (Eigen::VectorXi::LinSpaced(sp_mech_top3d_->sp_mesh_->GetNumEles(), 0, + sp_mech_top3d_->sp_mesh_->GetNumEles()) * + Eigen::RowVectorXi::Ones(sp_mech_top3d_->sp_mesh_->Get_DOFS_EACH_ELE())).reshaped(); + j_dFth_drho_ = sp_mech_top3d_->sp_mesh_->GetEleId2DofsMap().transpose().reshaped(); + + D0_ = sp_mech_top3d_->sp_fea_->computeD(1.0); + const Eigen::MatrixXd Be = sp_mech_top3d_->sp_fea_->computeBe(); + Inted_ = Be.transpose() * D0_ * (Eigen::VectorXd(6) << 1, 1, 1, 0, 0, 0).finished(); + Eigen::VectorXi i_Hs = Eigen::VectorXi::LinSpaced(sp_mech_top3d_->Hs_.size(), 0, + sp_mech_top3d_->Hs_.size()); + Eigen::SparseMatrix sp_Hs(i_Hs.size(), i_Hs.size()); + auto v_tri = Vec2Triplet(i_Hs, i_Hs, sp_mech_top3d_->Hs_); + sp_Hs.setFromTriplets(v_tri.begin(), v_tri.end()); + drho_dx_ = sp_Hs * sp_mech_top3d_->H_; + + + } + private: std::shared_ptr sp_mech_top3d_, sp_thermal_top3d_; + private: + Eigen::VectorXi i_dFth_dT_, j_dFth_dT_; + Eigen::VectorXi i_dFth_drho_, j_dFth_drho_; + Eigen::VectorXd Inted_; + Eigen::MatrixXd D0_; + Eigen::SparseMatrix drho_dx_; + }; } diff --git a/src/Top3d.cpp b/src/Top3d.cpp index d74b3d4..581d56b 100644 --- a/src/Top3d.cpp +++ b/src/Top3d.cpp @@ -190,7 +190,7 @@ void Top3d::Precompute() { jK_ = Eigen::KroneckerProduct(mat_ele2dofs, Eigen::RowVectorXi::Ones(dofs_each_ele)) .transpose() .reshaped(); - Ke_=sp_fea_->computeKe(); + Ke_=sp_fea_->computeKe(1.0); sKe_ = Ke_.reshaped(); // precompute filter diff --git a/src/Top3d.h b/src/Top3d.h index caf55b5..30e5e86 100644 --- a/src/Top3d.h +++ b/src/Top3d.h @@ -27,12 +27,19 @@ namespace da::sha { double penal = 1.0; int max_loop = 100; double r_min = 2.0; - + double T_ref = 295; + double T_limit=325; double tol_x = 0.01; double E_factor = 1e-9; + double R_E = 8; + double R_lambda = 0; + double R_beta = 0; }; + class ThermoelasticTop3d; + class Top3d { + friend ThermoelasticTop3d; public: Top3d(std::shared_ptr sp_para, std::shared_ptr sp_fea, std::shared_ptr sp_mesh) : sp_para_(sp_para @@ -40,7 +47,7 @@ namespace da::sha { F_ = SpMat(sp_mesh_->GetNumDofs(), 1); K_ = SpMat(sp_mesh_->GetNumDofs(), sp_mesh_->GetNumDofs()); spdlog::info("DOF: {}", K_.rows()); - spdlog::info("start to precompute..."); + spdlog::info("start to Precompute..."); Precompute(); } @@ -100,7 +107,10 @@ namespace da::sha { // sp_material_->computeN(Eigen::RowVector3d(local_coords(i,0),local_coords(i,1),local_coords(i,2)), N);//TODO: fixme Eigen::VectorXd node_forces = N.transpose() * forces; for (int dofi = 0; dofi < dof_of_the_ele.size(); ++dofi) + { + set_dofs_to_load.insert(dof_of_the_ele(i)); F_.coeffRef(dof_of_the_ele(i), 0) += node_forces(i); + } } } @@ -112,7 +122,10 @@ namespace da::sha { for (int i = 0; i < node_id_to_load.size(); ++i) { for (int dir = 0; dir < sp_mesh_->Get_DOFS_EACH_NODE(); ++dir) if (forces[dir]) + { + set_dofs_to_load.insert(GetDof(node_id_to_load[i], dir)); F_.coeffRef(GetDof(node_id_to_load[i], dir), 0) += forces(dir); + } } } @@ -141,7 +154,15 @@ namespace da::sha { for (auto dof_value: v_dofs_to_set) { auto [dof, value] = dof_value; K.coeffRef(dof, dof) *= 1e10; - F.coeffRef(dof, 0) = K.coeffRef(dof, dof)* value; + F.coeffRef(dof, 0) = K.coeffRef(dof, dof) * value; + } + } + + void IntroduceFixedDofs(Eigen::SparseMatrix &K, Eigen::VectorXd &F) { + for (auto dof_value: v_dofs_to_set) { + auto [dof, value] = dof_value; + K.coeffRef(dof, dof) *= 1e10; + F(dof) = K.coeffRef(dof, dof) * value; } } @@ -158,13 +179,13 @@ namespace da::sha { private: SpMat F_; std::vector> v_dofs_to_set; + std::set set_dofs_to_load; Eigen::VectorXi iK_, jK_; Eigen::VectorXd sKe_; Eigen::MatrixXd Ke_; SpMat K_; -// std::set dofs_to_fixed_; SpMat H_; Eigen::VectorXd Hs_; // result