diff --git a/.gitignore b/.gitignore index 438227d..556be6f 100644 --- a/.gitignore +++ b/.gitignore @@ -2,5 +2,4 @@ cmake-*/ build/ .idea/ -output/*/* -output/* \ No newline at end of file +output/**/* diff --git a/CMakeLists.txt b/CMakeLists.txt index 9f69899..e6294d3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,12 +3,9 @@ project(top3d LANGUAGES CUDA CXX) option(PROJECT_WITH_SIMD "Enable SIMD" ON) option(VERBOSE "Show more infos" ON) # choose linear solver -option(ENABLE_AMGCL "Use AMGCL" OFF) -if (ENABLE_AMGCL) - option(ENABLE_AMGCL_CUDA "use Cuda to speed up AMGCL" OFF) -else () - option(ENABLE_SUITESPARSE "Use SuiteSparse" ON) -endif () +option(ENABLE_AMGCL "Use AMGCL" ON) +option(ENABLE_AMGCL_CUDA "use Cuda to speed up AMGCL" ON) +option(ENABLE_SUITESPARSE "Use SuiteSparse" OFF) list(PREPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake ${CMAKE_CURRENT_SOURCE_DIR}/cmake/find) @@ -34,7 +31,9 @@ if (ENABLE_AMGCL AND ENABLE_AMGCL_CUDA) else () File(GLOB CPP_FILES src/*.cpp - src/*/*.cpp + src/LinearSolver/Amgcl.cpp + src/FEA/*.cpp + src/Mesh/*.cpp ) endif () include_directories(${CMAKE_SOURCE_DIR}/src) @@ -90,6 +89,7 @@ if (ENABLE_AMGCL) ) target_link_libraries(${PROJECT_NAME}_cuda_lib PUBLIC ${PROJECT_NAME}_lib) target_compile_options(${PROJECT_NAME}_cuda_lib PRIVATE -Xcompiler -fopenmp) + target_compile_definitions(${PROJECT_NAME}_cuda_lib PUBLIC USE_AMGCL_CUDA) endif () endif () @@ -98,4 +98,11 @@ target_compile_definitions(${PROJECT_NAME}_lib PUBLIC CMAKE_SOURCE_DIR="${CMAKE_ target_compile_definitions(${PROJECT_NAME}_lib PUBLIC WRITE_TENSOR) target_compile_definitions(${PROJECT_NAME}_lib PUBLIC DEBUG) +# if we want to use the temperature limit +#target_compile_definitions(${PROJECT_NAME}_lib PUBLIC WITH_T_LIMIT) +target_compile_definitions(${PROJECT_NAME}_lib PUBLIC MECH_ONLY) + +target_compile_definitions(${PROJECT_NAME}_lib PUBLIC OUTPUT_DIR="${CMAKE_SOURCE_DIR}/output") +target_compile_definitions(${PROJECT_NAME}_lib PUBLIC ASSETS_DIR="${CMAKE_SOURCE_DIR}/assets") + add_subdirectory(${CMAKE_SOURCE_DIR}/examples) \ No newline at end of file diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index cd7d281..f693283 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -5,6 +5,7 @@ foreach (EXAMPLE ${EXAMPLE_LIST}) if ((ENABLE_AMGCL AND ENABLE_AMGCL_CUDA AND ${EXAMPLE} MATCHES "cuda$") OR (NOT ENABLE_AMGCL_CUDA AND NOT ${EXAMPLE} MATCHES "cuda$")) add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/${EXAMPLE}) + message(STATUS "Added example: ${EXAMPLE}") endif () endif () endif () diff --git a/examples/top-thermoelastic-BiclampedStructure-cuda/CMakeLists.txt b/examples/top-thermoelastic-BiclampedStructure-cuda/CMakeLists.txt index 72d876c..5192e2c 100644 --- a/examples/top-thermoelastic-BiclampedStructure-cuda/CMakeLists.txt +++ b/examples/top-thermoelastic-BiclampedStructure-cuda/CMakeLists.txt @@ -4,5 +4,3 @@ add_executable(${SUB_PROJECT_NAME} main.cu) target_compile_options(${SUB_PROJECT_NAME} PRIVATE -Xcompiler -fopenmp) target_link_libraries(${SUB_PROJECT_NAME} PUBLIC ${PROJECT_NAME}_cuda_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") diff --git a/examples/top-thermoelastic-BiclampedStructure/CMakeLists.txt b/examples/top-thermoelastic-BiclampedStructure/CMakeLists.txt index 50d9b5e..1b5cd82 100644 --- a/examples/top-thermoelastic-BiclampedStructure/CMakeLists.txt +++ b/examples/top-thermoelastic-BiclampedStructure/CMakeLists.txt @@ -1,6 +1,4 @@ set(SUB_PROJECT_NAME top-thermoelastic-BiclampedStructure) 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 CONFIG_FILE="${CMAKE_CURRENT_SOURCE_DIR}/config.json") \ No newline at end of file diff --git a/examples/top-thermoelastic-Lshape-condition/CMakeLists.txt b/examples/top-thermoelastic-Lshape-condition/CMakeLists.txt index 53c2766..00db91c 100644 --- a/examples/top-thermoelastic-Lshape-condition/CMakeLists.txt +++ b/examples/top-thermoelastic-Lshape-condition/CMakeLists.txt @@ -2,5 +2,3 @@ set(SUB_PROJECT_NAME top-thermoelastic-Lshape-condition) 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") diff --git a/examples/top-thermoelastic-Lshape-condition/main.cpp b/examples/top-thermoelastic-Lshape-condition/main.cpp index 324daf7..d737103 100644 --- a/examples/top-thermoelastic-Lshape-condition/main.cpp +++ b/examples/top-thermoelastic-Lshape-condition/main.cpp @@ -72,7 +72,8 @@ int main() { 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) { + int the_len=len_x; + if (j > the_len && k > the_len) { L_shape_model(i, j, k) = 0; } } diff --git a/examples/top-thermoelastic-beam-cuda/CMakeLists.txt b/examples/top-thermoelastic-beam-cuda/CMakeLists.txt new file mode 100644 index 0000000..21cc516 --- /dev/null +++ b/examples/top-thermoelastic-beam-cuda/CMakeLists.txt @@ -0,0 +1,6 @@ +set(SUB_PROJECT_NAME top-thermoelastic-beam-cu) +add_executable(${SUB_PROJECT_NAME} main.cu) +# OpenMP Required! +target_compile_options(${SUB_PROJECT_NAME} PRIVATE -Xcompiler -fopenmp) +target_link_libraries(${SUB_PROJECT_NAME} PUBLIC ${PROJECT_NAME}_cuda_lib) +target_compile_definitions(${SUB_PROJECT_NAME} PUBLIC CONFIG_FILE="${CMAKE_CURRENT_SOURCE_DIR}/config.json") diff --git a/examples/top-thermoelastic-beam-cuda/config.json b/examples/top-thermoelastic-beam-cuda/config.json new file mode 100644 index 0000000..cab9237 --- /dev/null +++ b/examples/top-thermoelastic-beam-cuda/config.json @@ -0,0 +1,67 @@ +{ + "TopologyOptimizationExample": "beam","//": "user defined project name", + "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.15, + "penal": 3.0, + "r_min": 1.2, + "T_ref": 312, + "T_limit": 3996.52, + "R_E":28, + "R_lambda":28, + "R_beta":0 + }, + "model": { + "regular_model": { + "lx": 2, + "ly": 30, + "lz": 10 + } + }, + "mechanical_boundary_condition":{ + "DBC": [ + { + "min": [0, 0, 0], + "max": [1, 0, 1], + "dir": [1, 1, 1] + } + + ], + "NBC": [ + { + "min": [0.5, 1, 0], + "max": [0.5, 1, 0], + "val": [0.0,0.0 , -1e7] + } + + ] + }, + "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-beam-cuda/main.cu b/examples/top-thermoelastic-beam-cuda/main.cu new file mode 100644 index 0000000..496ff1d --- /dev/null +++ b/examples/top-thermoelastic-beam-cuda/main.cu @@ -0,0 +1,210 @@ +// +// 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("assets 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); + // read title + std::string condition = j_config["TopologyOptimizationExample"]; + spdlog::critical("TopologyOptimizationExample: {}", condition); + + // 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; + if (false/* only for Lshape condition*/) { + // L-shape condition + spdlog::info("Using Lshape condition!"); + top::Tensor3d L_shape_model(len_x, len_y, len_z); + L_shape_model.setConstant(1); + // set the initial voxel model + for (int k = 0; k < len_z; ++k) { + for (int j = 0; j < len_y; ++j) { + for (int i = 0; i < len_x; ++i) { + int the_len=len_x; + if (j > the_len && k > the_len) { + 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(); + top::Tensor3d ten_rho =sp_mech_top3d->TopOptMainLoop(); + + // extract rho (txt or vtk) + fs_path rho_txt_path = output_dir / "txt" / (condition + "_rho.txt"); + write_tensor3d(rho_txt_path, ten_rho, sp_mech_mesh->GetOrigin(), + sp_mech_mesh->GetOrigin() + sp_mech_mesh->GetLenBox()); + spdlog::info("write density txt to: {}", rho_txt_path.c_str()); + + fs_path rho_vtk_path = output_dir / "vtk" / (condition + "_rho.vtk"); + fs_path rho_vtk_path_threshold = output_dir / "vtk" / (condition + "_rho_01.vtk"); + WriteTensorToVtk(rho_vtk_path, ten_rho, sp_mech_mesh); + WriteTensorToVtk(rho_vtk_path_threshold, ten_rho, sp_mech_mesh, 0.5); + spdlog::info("write density vtk to: {}", rho_vtk_path.c_str()); + spdlog::info("write density(0 or 1) vtk to: {}", rho_vtk_path_threshold.c_str()); + + // extract temperature(vtk) + fs_path U_vtk_path = output_dir / "vtk" / (condition + "_U.vtk"); + WriteUToVtk(U_vtk_path, ther_top3d.GetTemperature(), sp_mech_mesh); + spdlog::info("write displacement norm vtk to: {}", U_vtk_path.c_str()); + + // extract stress field(txt or vtk) + auto v_tenosr = ther_top3d.GetTensorOfStress(Eigen::Vector3d{0, 1, 2}); + fs_path stressX_txt_pah = output_dir / "txt" / (condition + "_stressX.txt"); + write_tensor3d(stressX_txt_pah, v_tenosr[0], sp_mech_mesh->GetOrigin(), + sp_mech_mesh->GetOrigin() + sp_mech_mesh->GetLenBox()); + spdlog::info("write stressX txt to: {}", stressX_txt_pah.c_str()); + fs_path stressX_vtk_path = output_dir / "vtk" / (condition + "_stressX.vtk"); + WriteTensorToVtk(stressX_vtk_path, v_tenosr[0], sp_mech_mesh); + spdlog::info("write stressX vtk to: {}", stressX_vtk_path.c_str()); + + fs_path stressY_txt_pah = output_dir / "txt" / (condition + "_stressY.txt"); + write_tensor3d(stressY_txt_pah, v_tenosr[1], sp_mech_mesh->GetOrigin(), + sp_mech_mesh->GetOrigin() + sp_mech_mesh->GetLenBox()); + spdlog::info("write stressY txt to: {}", stressY_txt_pah.c_str()); + fs_path stressY_vtk_path = output_dir / "vtk" / (condition + "_stressY.vtk"); + WriteTensorToVtk(stressY_vtk_path, v_tenosr[1], sp_mech_mesh); + spdlog::info("write stressY vtk to: {}", stressY_vtk_path.c_str()); + + fs_path stressZ_txt_pah = output_dir / "txt" / (condition + "_stressZ.txt"); + write_tensor3d(stressZ_txt_pah, v_tenosr[2], sp_mech_mesh->GetOrigin(), + sp_mech_mesh->GetOrigin() + sp_mech_mesh->GetLenBox()); + spdlog::info("write stressZ txt to: {}", stressZ_txt_pah.c_str()); + fs_path stressZ_vtk_path = output_dir / "vtk" / (condition + "_stressZ.vtk"); + WriteTensorToVtk(stressZ_vtk_path, v_tenosr[2], sp_mech_mesh); + spdlog::info("write stressZ vtk to: {}", stressZ_vtk_path.c_str()); + +} diff --git a/examples/top-thermoelastic-beam/CMakeLists.txt b/examples/top-thermoelastic-beam/CMakeLists.txt new file mode 100644 index 0000000..31b66d4 --- /dev/null +++ b/examples/top-thermoelastic-beam/CMakeLists.txt @@ -0,0 +1,4 @@ +set(SUB_PROJECT_NAME top-thermoelastic-beam) +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") diff --git a/examples/top-thermoelastic-beam/config.json b/examples/top-thermoelastic-beam/config.json new file mode 100644 index 0000000..0bc95ef --- /dev/null +++ b/examples/top-thermoelastic-beam/config.json @@ -0,0 +1,67 @@ +{ + "TopologyOptimizationExample": "beam","//": "user defined project name", + "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": 100, + "volfrac": 0.3, + "penal": 3.0, + "r_min": 1.5, + "T_ref": 312, + "T_limit": 3996.52, + "R_E": 28, + "R_lambda": 28, + "R_beta":0 + }, + "model": { + "regular_model": { + "lx": 2, + "ly": 30, + "lz": 10 + } + }, + "mechanical_boundary_condition":{ + "DBC": [ + { + "min": [0, 0, 0], + "max": [1, 0, 1], + "dir": [1, 1, 1] + } + + ], + "NBC": [ + { + "min": [0, 1, 0], + "max": [1, 1, 0], + "val": [0.0,0.0 , -1e11] + } + + ] + }, + "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-beam/main.cpp b/examples/top-thermoelastic-beam/main.cpp new file mode 100644 index 0000000..053effd --- /dev/null +++ b/examples/top-thermoelastic-beam/main.cpp @@ -0,0 +1,206 @@ +// +// 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("assets 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); + // read title + std::string condition = j_config["TopologyOptimizationExample"]; + spdlog::critical("TopologyOptimizationExample: {}", condition); + + // 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; + if (false/* only for Lshape condition*/) { + // L-shape condition + spdlog::info("Using Lshape condition!"); + top::Tensor3d L_shape_model(len_x, len_y, len_z); + L_shape_model.setConstant(1); + // set the initial voxel model + for (int k = 0; k < len_z; ++k) { + for (int j = 0; j < len_y; ++j) { + for (int i = 0; i < len_x; ++i) { + int the_len=len_x; + if (j > the_len && k > the_len) { + 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 rho (txt or vtk) + fs_path rho_txt_path = output_dir / "txt" / (condition + "_rho.txt"); + write_tensor3d(rho_txt_path, ten_rho, sp_mech_mesh->GetOrigin(), + sp_mech_mesh->GetOrigin() + sp_mech_mesh->GetLenBox()); + spdlog::info("write density txt to: {}", rho_txt_path.c_str()); + + fs_path rho_vtk_path = output_dir / "vtk" / (condition + "_rho.vtk"); + WriteTensorToVtk(rho_vtk_path, ten_rho, sp_mech_mesh); + spdlog::info("write density vtk to: {}", rho_vtk_path.c_str()); + + // extract temperature(vtk) + fs_path U_vtk_path = output_dir / "vtk" / (condition + "_U.vtk"); + WriteUToVtk(U_vtk_path, ther_top3d.GetTemperature(), sp_mech_mesh); + spdlog::info("write displacement norm vtk to: {}", U_vtk_path.c_str()); + + // extract stress field(txt or vtk) + auto v_tenosr = ther_top3d.GetTensorOfStress(Eigen::Vector3d{0, 1, 2}); + fs_path stressX_txt_pah = output_dir / "txt" / (condition + "_stressX.txt"); + write_tensor3d(stressX_txt_pah, v_tenosr[0], sp_mech_mesh->GetOrigin(), + sp_mech_mesh->GetOrigin() + sp_mech_mesh->GetLenBox()); + spdlog::info("write stressX txt to: {}", stressX_txt_pah.c_str()); + fs_path stressX_vtk_path = output_dir / "vtk" / (condition + "_stressX.vtk"); + WriteTensorToVtk(stressX_vtk_path, v_tenosr[0], sp_mech_mesh); + spdlog::info("write stressX vtk to: {}", stressX_vtk_path.c_str()); + + fs_path stressY_txt_pah = output_dir / "txt" / (condition + "_stressY.txt"); + write_tensor3d(stressY_txt_pah, v_tenosr[1], sp_mech_mesh->GetOrigin(), + sp_mech_mesh->GetOrigin() + sp_mech_mesh->GetLenBox()); + spdlog::info("write stressY txt to: {}", stressY_txt_pah.c_str()); + fs_path stressY_vtk_path = output_dir / "vtk" / (condition + "_stressY.vtk"); + WriteTensorToVtk(stressY_vtk_path, v_tenosr[1], sp_mech_mesh); + spdlog::info("write stressY vtk to: {}", stressY_vtk_path.c_str()); + + fs_path stressZ_txt_pah = output_dir / "txt" / (condition + "_stressZ.txt"); + write_tensor3d(stressZ_txt_pah, v_tenosr[2], sp_mech_mesh->GetOrigin(), + sp_mech_mesh->GetOrigin() + sp_mech_mesh->GetLenBox()); + spdlog::info("write stressZ txt to: {}", stressZ_txt_pah.c_str()); + fs_path stressZ_vtk_path = output_dir / "vtk" / (condition + "_stressZ.vtk"); + WriteTensorToVtk(stressZ_vtk_path, v_tenosr[2], sp_mech_mesh); + spdlog::info("write stressZ vtk to: {}", stressZ_vtk_path.c_str()); + +} diff --git a/src/LinearSolver/AmgclCuda.h b/src/LinearSolver/AmgclCuda.h index 0a104a2..5ddf425 100644 --- a/src/LinearSolver/AmgclCuda.h +++ b/src/LinearSolver/AmgclCuda.h @@ -41,8 +41,8 @@ namespace da { static typename Solver::params &GetPrm() { static typename Solver::params prm; - prm.solver.tol = 1e-10; - prm.solver.maxiter = 1500; + prm.solver.tol = 1e-6; + prm.solver.maxiter = 6000; return prm; } @@ -69,7 +69,7 @@ namespace da { thrust::device_vector d_rhs(rhs); thrust::device_vector d_x(n_, 0.0); auto [iter, error] = solver_(d_rhs, d_x); - if ( verbose && (error > 1e-4 || iter > 500)) { + if ( verbose && (error > 1e-6)) { std::cout << "tol: " << solver_.solver().prm.tol << '\t' << "solve_iter: " << iter << diff --git a/src/ThermoelasticTop3d.cpp b/src/ThermoelasticTop3d.cpp index 7e03c65..4d425c9 100644 --- a/src/ThermoelasticTop3d.cpp +++ b/src/ThermoelasticTop3d.cpp @@ -4,7 +4,10 @@ #include "ThermoelasticTop3d.h" #include "LinearSolver/Amgcl.h" + +#ifdef USE_AMGCL_CUDA #include "LinearSolver/AmgclCuda.h" +#endif da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { auto &sp_mesh_ = sp_mech_top3d_->sp_mesh_; @@ -38,7 +41,8 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { dv = (sp_mech_top3d_->H_ * dv).array() / sp_mech_top3d_->Hs_.array(); Eigen::VectorXd ele_to_write = - Eigen::VectorXd::Zero(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz()); + Eigen::VectorXd::Zero( + sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz()); Eigen::VectorXi global_idx_of_ele_in_use = sp_mesh_->GetGlobalIdxOfEleInUse(); // dofs of limited T @@ -65,7 +69,8 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { 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)); + map_ele2Limit[i].push_back( + LimitedDof(v_dof[k], k, j)); } } } @@ -80,36 +85,45 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { while (change > sp_para_->tol_x && loop < sp_para_->max_loop) { ++loop; // filter - xPhys_col = (sp_mech_top3d_->H_ * xPhys_col).array() / sp_mech_top3d_->Hs_.array(); + xPhys_col = (sp_mech_top3d_->H_ * xPhys_col).array() / + sp_mech_top3d_->Hs_.array(); 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 { + 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 CalDRDrho_Vec = [](const Eigen::VectorXd &vec_rho, + double R) -> Eigen::VectorXd { auto down = 1 + R * (1 - vec_rho.array()); return (1 + R) / (down * down); }; auto CalE_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { - return E_min + CalR_Vec(vec_rho, sp_para_->R_E).array() * (E0_m - E_min); + return E_min + + CalR_Vec(vec_rho, sp_para_->R_E).array() * (E0_m - E_min); }; - auto CalDEDrho_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { + auto CalDEDrho_Vec = [&]( + const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { return CalDRDrho_Vec(vec_rho, sp_para_->R_E) * (E0_m - E_min); }; - auto CalLambda_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { + 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); + 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 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) { @@ -125,10 +139,12 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { Eigen::VectorXd sK_th = (sp_thermal_top3d_->sKe_ * CalLambda_Vec(xPhys_col).transpose()) .reshaped(); - auto v_tri_th = Vec2Triplet(sp_thermal_top3d_->iK_, sp_thermal_top3d_->jK_, sK_th); + 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_); - INIT_SOLVER(solver_th,sp_thermal_top3d_->K_); + sp_thermal_top3d_->IntroduceFixedDofs(sp_thermal_top3d_->K_, + sp_thermal_top3d_->F_); + INIT_SOLVER(solver_th, sp_thermal_top3d_->K_); sp_thermal_top3d_->U_ = solver_th.solve(sp_thermal_top3d_->F_); #endif @@ -137,22 +153,27 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { Eigen::VectorXd sK_m = (sp_mech_top3d_->sKe_ * CalE_Vec(xPhys_col).transpose()) .reshaped(); - auto v_tri_m = Vec2Triplet(sp_mech_top3d_->iK_, sp_mech_top3d_->jK_, sK_m); + 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()); #ifndef MECH_ONLY // F_th Eigen::VectorXd &T = sp_thermal_top3d_->U_; - Eigen::VectorXd F_th = Eigen::VectorXd::Zero(sp_mech_top3d_->sp_mesh_->GetNumDofs()); + 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_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_; + F_th(dofs_m) += + beta_rho * (Te - sp_mech_top3d_->sp_para_->T_ref) * Inted_; } if (loop == 1 || loop == sp_para_->max_loop) - spdlog::info("||Fth|| / ||Fm||: {}", F_th.norm() / sp_mech_top3d_->F_.norm()); + spdlog::info("||Fth|| / ||Fm||: {}", + F_th.norm() / sp_mech_top3d_->F_.norm()); #endif Eigen::VectorXd F = Eigen::VectorXd(sp_mech_top3d_->F_) #ifndef MECH_ONLY @@ -163,7 +184,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { sp_mech_top3d_->IntroduceFixedDofs(sp_mech_top3d_->K_, F); // Timer::tic(); - INIT_SOLVER(solver,sp_mech_top3d_->K_); + INIT_SOLVER(solver, sp_mech_top3d_->K_); sp_mech_top3d_->U_ = solver.solve(F); // Timer::toc(); @@ -171,43 +192,52 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { // compliance Eigen::VectorXd ce(sp_mech_top3d_->sp_mesh_->GetNumEles()); for (int i = 0; i < sp_mech_top3d_->sp_mesh_->GetNumEles(); ++i) { - Eigen::VectorXi dofs_in_ele_i = sp_mech_top3d_->sp_mesh_->MapEleId2Dofs(i); + Eigen::VectorXi dofs_in_ele_i = sp_mech_top3d_->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() * CalE_Vec(xPhys_col); - double v = flg_chosen ? xPhys_col(chosen_ele_id).sum() : xPhys_col.sum(); + double v = flg_chosen ? xPhys_col(chosen_ele_id).sum() + : xPhys_col.sum(); // sensitivity // lambda_m Eigen::VectorXd lambda_m = -sp_mech_top3d_->U_; #ifndef MECH_ONLY // dFth_drho - Eigen::SparseMatrix dFth_drho(sp_mech_top3d_->sp_mesh_->GetNumEles(), - sp_mech_top3d_->sp_mesh_->GetNumDofs()); + 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_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) * + 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; + 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()); + 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::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); @@ -215,8 +245,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::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 * + 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(); assert(ele_dFth_dT.rows() == 8 && ele_dFth_dT.cols() == 24); v_dFth_dT(Eigen::seqN(i * ele_dFth_dT.size(), @@ -240,7 +273,8 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { // 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::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; @@ -261,13 +295,18 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { 2 * Eigen::VectorXd(dFth_drho * sp_mech_top3d_->U_); #endif #ifndef MECH_ONLY +#ifdef WITH_T_LIMIT // dT_drho - auto CalDKth_Mul_T__For_DTi_Drhoj = [&](int rho_i) -> Eigen::SparseVector { - Eigen::VectorXi dofs_in_ele = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs(rho_i); + auto CalDKth_Mul_T__For_DTi_Drhoj = [&]( + int rho_i) -> Eigen::SparseVector { + Eigen::VectorXi dofs_in_ele = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs( + rho_i); Eigen::VectorXd ele_T = sp_thermal_top3d_->U_(dofs_in_ele); - Eigen::SparseVector sp_dKth_Mul_T(sp_thermal_top3d_->sp_mesh_->GetNumDofs()); + Eigen::SparseVector sp_dKth_Mul_T( + sp_thermal_top3d_->sp_mesh_->GetNumDofs()); Eigen::VectorXd dKe_th_Mul_T = - CalDlambdaDrho(xPhys_col(rho_i)) * sp_thermal_top3d_->Ke_ * ele_T; + CalDlambdaDrho(xPhys_col(rho_i)) * sp_thermal_top3d_->Ke_ * + ele_T; for (int i = 0; i < dofs_in_ele.size(); ++i) { sp_dKth_Mul_T.coeffRef(dofs_in_ele(i)) = dKe_th_Mul_T(i); } @@ -284,18 +323,21 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { Eigen::VectorXd lambda_i = solver_th.solve(Li); return lambda_i.transpose() * sp_dKth_Mul_T; }; - Eigen::MatrixXd dT_drho = Eigen::MatrixXd::Zero(sp_thermal_top3d_->sp_mesh_->GetNumEles(), - sp_thermal_top3d_->set_dofs_to_load.size()); + 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; auto sp_dKth_Mul_T = CalDKth_Mul_T__For_DTi_Drhoj(ele_id); for (auto &limited: v_limited) { - dT_drho(ele_id, limited.idx_of_load_dof) = CalDTi_Drhoj(limited.dof, sp_dKth_Mul_T); + dT_drho(ele_id, limited.idx_of_load_dof) = CalDTi_Drhoj( + limited.dof, sp_dKth_Mul_T); } } // dT_dx Eigen::MatrixXd dT_dx = sp_mech_top3d_->drho_dx_ * dT_drho; +#endif #endif // dc_dx Eigen::VectorXd dc_dx = sp_mech_top3d_->drho_dx_ * dc_drho; @@ -305,68 +347,68 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { size_t num_constraints = 1 #ifndef MECH_ONLY + #ifdef WITH_T_LIMIT + dT_dx.cols()// volume and temperature constraints +#endif #endif ; - size_t num_variables = flg_chosen ? chosen_ele_id.size() : sp_mesh_->GetNumEles(); + 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; + 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() * SENSITIVITY_SCALE_COEF - : dc_dx / dc_dx.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(num_constraints) + << (v / num_variables - sp_para_->volfrac) #ifndef MECH_ONLY +#ifdef WITH_T_LIMIT , T(v_dof).array() / sp_para_->T_limit - 1 +#endif #endif ).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 #ifndef MECH_ONLY +#ifdef WITH_T_LIMIT , 1.0 / sp_para_->T_limit * dT_dx +#endif #endif ).finished().transpose() * SENSITIVITY_SCALE_COEF; - static Eigen::VectorXd low_bounds = Eigen::VectorXd::Zero(num_variables); + 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(), + 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(); + 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, + 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(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); - } - 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; @@ -374,31 +416,38 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { assert(xPhys_col.size()); Eigen::VectorXi continue_idx = Eigen::VectorXi::LinSpaced(xPhys_col.size(), 0, xPhys_col.size()); - Eigen::VectorXi unchosen_idx = flg_chosen ? SetDifference(continue_idx, chosen_ele_id) + Eigen::VectorXi unchosen_idx = flg_chosen ? SetDifference(continue_idx, + chosen_ele_id) : Eigen::VectorXi(); { xPhys_col(unchosen_idx).setZero(); 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); + 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()}); + 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(global_idx_of_ele_in_use) = xPhys_col; - Tensor3d ten_xPhys_to_write(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz(), 1, - 1); + 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()}); + 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; } diff --git a/src/ThermoelasticTop3d.cu b/src/ThermoelasticTop3d.cu index 7e03c65..4d425c9 100644 --- a/src/ThermoelasticTop3d.cu +++ b/src/ThermoelasticTop3d.cu @@ -4,7 +4,10 @@ #include "ThermoelasticTop3d.h" #include "LinearSolver/Amgcl.h" + +#ifdef USE_AMGCL_CUDA #include "LinearSolver/AmgclCuda.h" +#endif da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { auto &sp_mesh_ = sp_mech_top3d_->sp_mesh_; @@ -38,7 +41,8 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { dv = (sp_mech_top3d_->H_ * dv).array() / sp_mech_top3d_->Hs_.array(); Eigen::VectorXd ele_to_write = - Eigen::VectorXd::Zero(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz()); + Eigen::VectorXd::Zero( + sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz()); Eigen::VectorXi global_idx_of_ele_in_use = sp_mesh_->GetGlobalIdxOfEleInUse(); // dofs of limited T @@ -65,7 +69,8 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { 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)); + map_ele2Limit[i].push_back( + LimitedDof(v_dof[k], k, j)); } } } @@ -80,36 +85,45 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { while (change > sp_para_->tol_x && loop < sp_para_->max_loop) { ++loop; // filter - xPhys_col = (sp_mech_top3d_->H_ * xPhys_col).array() / sp_mech_top3d_->Hs_.array(); + xPhys_col = (sp_mech_top3d_->H_ * xPhys_col).array() / + sp_mech_top3d_->Hs_.array(); 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 { + 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 CalDRDrho_Vec = [](const Eigen::VectorXd &vec_rho, + double R) -> Eigen::VectorXd { auto down = 1 + R * (1 - vec_rho.array()); return (1 + R) / (down * down); }; auto CalE_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { - return E_min + CalR_Vec(vec_rho, sp_para_->R_E).array() * (E0_m - E_min); + return E_min + + CalR_Vec(vec_rho, sp_para_->R_E).array() * (E0_m - E_min); }; - auto CalDEDrho_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { + auto CalDEDrho_Vec = [&]( + const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { return CalDRDrho_Vec(vec_rho, sp_para_->R_E) * (E0_m - E_min); }; - auto CalLambda_Vec = [&](const Eigen::VectorXd &vec_rho) -> Eigen::VectorXd { + 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); + 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 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) { @@ -125,10 +139,12 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { Eigen::VectorXd sK_th = (sp_thermal_top3d_->sKe_ * CalLambda_Vec(xPhys_col).transpose()) .reshaped(); - auto v_tri_th = Vec2Triplet(sp_thermal_top3d_->iK_, sp_thermal_top3d_->jK_, sK_th); + 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_); - INIT_SOLVER(solver_th,sp_thermal_top3d_->K_); + sp_thermal_top3d_->IntroduceFixedDofs(sp_thermal_top3d_->K_, + sp_thermal_top3d_->F_); + INIT_SOLVER(solver_th, sp_thermal_top3d_->K_); sp_thermal_top3d_->U_ = solver_th.solve(sp_thermal_top3d_->F_); #endif @@ -137,22 +153,27 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { Eigen::VectorXd sK_m = (sp_mech_top3d_->sKe_ * CalE_Vec(xPhys_col).transpose()) .reshaped(); - auto v_tri_m = Vec2Triplet(sp_mech_top3d_->iK_, sp_mech_top3d_->jK_, sK_m); + 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()); #ifndef MECH_ONLY // F_th Eigen::VectorXd &T = sp_thermal_top3d_->U_; - Eigen::VectorXd F_th = Eigen::VectorXd::Zero(sp_mech_top3d_->sp_mesh_->GetNumDofs()); + 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_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_; + F_th(dofs_m) += + beta_rho * (Te - sp_mech_top3d_->sp_para_->T_ref) * Inted_; } if (loop == 1 || loop == sp_para_->max_loop) - spdlog::info("||Fth|| / ||Fm||: {}", F_th.norm() / sp_mech_top3d_->F_.norm()); + spdlog::info("||Fth|| / ||Fm||: {}", + F_th.norm() / sp_mech_top3d_->F_.norm()); #endif Eigen::VectorXd F = Eigen::VectorXd(sp_mech_top3d_->F_) #ifndef MECH_ONLY @@ -163,7 +184,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { sp_mech_top3d_->IntroduceFixedDofs(sp_mech_top3d_->K_, F); // Timer::tic(); - INIT_SOLVER(solver,sp_mech_top3d_->K_); + INIT_SOLVER(solver, sp_mech_top3d_->K_); sp_mech_top3d_->U_ = solver.solve(F); // Timer::toc(); @@ -171,43 +192,52 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { // compliance Eigen::VectorXd ce(sp_mech_top3d_->sp_mesh_->GetNumEles()); for (int i = 0; i < sp_mech_top3d_->sp_mesh_->GetNumEles(); ++i) { - Eigen::VectorXi dofs_in_ele_i = sp_mech_top3d_->sp_mesh_->MapEleId2Dofs(i); + Eigen::VectorXi dofs_in_ele_i = sp_mech_top3d_->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() * CalE_Vec(xPhys_col); - double v = flg_chosen ? xPhys_col(chosen_ele_id).sum() : xPhys_col.sum(); + double v = flg_chosen ? xPhys_col(chosen_ele_id).sum() + : xPhys_col.sum(); // sensitivity // lambda_m Eigen::VectorXd lambda_m = -sp_mech_top3d_->U_; #ifndef MECH_ONLY // dFth_drho - Eigen::SparseMatrix dFth_drho(sp_mech_top3d_->sp_mesh_->GetNumEles(), - sp_mech_top3d_->sp_mesh_->GetNumDofs()); + 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_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) * + 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; + 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()); + 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::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); @@ -215,8 +245,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::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 * + 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(); assert(ele_dFth_dT.rows() == 8 && ele_dFth_dT.cols() == 24); v_dFth_dT(Eigen::seqN(i * ele_dFth_dT.size(), @@ -240,7 +273,8 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { // 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::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; @@ -261,13 +295,18 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { 2 * Eigen::VectorXd(dFth_drho * sp_mech_top3d_->U_); #endif #ifndef MECH_ONLY +#ifdef WITH_T_LIMIT // dT_drho - auto CalDKth_Mul_T__For_DTi_Drhoj = [&](int rho_i) -> Eigen::SparseVector { - Eigen::VectorXi dofs_in_ele = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs(rho_i); + auto CalDKth_Mul_T__For_DTi_Drhoj = [&]( + int rho_i) -> Eigen::SparseVector { + Eigen::VectorXi dofs_in_ele = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs( + rho_i); Eigen::VectorXd ele_T = sp_thermal_top3d_->U_(dofs_in_ele); - Eigen::SparseVector sp_dKth_Mul_T(sp_thermal_top3d_->sp_mesh_->GetNumDofs()); + Eigen::SparseVector sp_dKth_Mul_T( + sp_thermal_top3d_->sp_mesh_->GetNumDofs()); Eigen::VectorXd dKe_th_Mul_T = - CalDlambdaDrho(xPhys_col(rho_i)) * sp_thermal_top3d_->Ke_ * ele_T; + CalDlambdaDrho(xPhys_col(rho_i)) * sp_thermal_top3d_->Ke_ * + ele_T; for (int i = 0; i < dofs_in_ele.size(); ++i) { sp_dKth_Mul_T.coeffRef(dofs_in_ele(i)) = dKe_th_Mul_T(i); } @@ -284,18 +323,21 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { Eigen::VectorXd lambda_i = solver_th.solve(Li); return lambda_i.transpose() * sp_dKth_Mul_T; }; - Eigen::MatrixXd dT_drho = Eigen::MatrixXd::Zero(sp_thermal_top3d_->sp_mesh_->GetNumEles(), - sp_thermal_top3d_->set_dofs_to_load.size()); + 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; auto sp_dKth_Mul_T = CalDKth_Mul_T__For_DTi_Drhoj(ele_id); for (auto &limited: v_limited) { - dT_drho(ele_id, limited.idx_of_load_dof) = CalDTi_Drhoj(limited.dof, sp_dKth_Mul_T); + dT_drho(ele_id, limited.idx_of_load_dof) = CalDTi_Drhoj( + limited.dof, sp_dKth_Mul_T); } } // dT_dx Eigen::MatrixXd dT_dx = sp_mech_top3d_->drho_dx_ * dT_drho; +#endif #endif // dc_dx Eigen::VectorXd dc_dx = sp_mech_top3d_->drho_dx_ * dc_drho; @@ -305,68 +347,68 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { size_t num_constraints = 1 #ifndef MECH_ONLY + #ifdef WITH_T_LIMIT + dT_dx.cols()// volume and temperature constraints +#endif #endif ; - size_t num_variables = flg_chosen ? chosen_ele_id.size() : sp_mesh_->GetNumEles(); + 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; + 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() * SENSITIVITY_SCALE_COEF - : dc_dx / dc_dx.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(num_constraints) + << (v / num_variables - sp_para_->volfrac) #ifndef MECH_ONLY +#ifdef WITH_T_LIMIT , T(v_dof).array() / sp_para_->T_limit - 1 +#endif #endif ).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 #ifndef MECH_ONLY +#ifdef WITH_T_LIMIT , 1.0 / sp_para_->T_limit * dT_dx +#endif #endif ).finished().transpose() * SENSITIVITY_SCALE_COEF; - static Eigen::VectorXd low_bounds = Eigen::VectorXd::Zero(num_variables); + 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(), + 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(); + 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, + 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(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); - } - 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; @@ -374,31 +416,38 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { assert(xPhys_col.size()); Eigen::VectorXi continue_idx = Eigen::VectorXi::LinSpaced(xPhys_col.size(), 0, xPhys_col.size()); - Eigen::VectorXi unchosen_idx = flg_chosen ? SetDifference(continue_idx, chosen_ele_id) + Eigen::VectorXi unchosen_idx = flg_chosen ? SetDifference(continue_idx, + chosen_ele_id) : Eigen::VectorXi(); { xPhys_col(unchosen_idx).setZero(); 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); + 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()}); + 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(global_idx_of_ele_in_use) = xPhys_col; - Tensor3d ten_xPhys_to_write(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz(), 1, - 1); + 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()}); + 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; } diff --git a/src/Top3d.cpp b/src/Top3d.cpp index 5952615..bb6daa1 100644 --- a/src/Top3d.cpp +++ b/src/Top3d.cpp @@ -8,7 +8,9 @@ #include "Eigen/src/Core/Matrix.h" #include "Util.h" #include "LinearSolver/Amgcl.h" +#ifdef USE_AMGCL_CUDA #include "LinearSolver/AmgclCuda.h" +#endif namespace da::sha { namespace top { Tensor3d Top3d::TopOptMainLoop() { @@ -105,19 +107,7 @@ namespace da::sha { spdlog::critical("Iter: {:3d}, Comp: {:.3e}, Vol: {:.2f}, Change: {:f}", loop, c, v, change); -#ifdef WRITE_TENSOR_IN_LOOP - // extract vtk - 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()}); - top::WriteTensorToVtk( - da::WorkingResultDirectoryPath() / ("field_matrix" + std::to_string(loop) + ".vtk"), - ten_xPhys_to_write, sp_mesh_); -#endif + } // result rho_ = xPhys_col; diff --git a/src/Top3d.cu b/src/Top3d.cu index 5952615..bb6daa1 100644 --- a/src/Top3d.cu +++ b/src/Top3d.cu @@ -8,7 +8,9 @@ #include "Eigen/src/Core/Matrix.h" #include "Util.h" #include "LinearSolver/Amgcl.h" +#ifdef USE_AMGCL_CUDA #include "LinearSolver/AmgclCuda.h" +#endif namespace da::sha { namespace top { Tensor3d Top3d::TopOptMainLoop() { @@ -105,19 +107,7 @@ namespace da::sha { spdlog::critical("Iter: {:3d}, Comp: {:.3e}, Vol: {:.2f}, Change: {:f}", loop, c, v, change); -#ifdef WRITE_TENSOR_IN_LOOP - // extract vtk - 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()}); - top::WriteTensorToVtk( - da::WorkingResultDirectoryPath() / ("field_matrix" + std::to_string(loop) + ".vtk"), - ten_xPhys_to_write, sp_mesh_); -#endif + } // result rho_ = xPhys_col; diff --git a/src/Util.cpp b/src/Util.cpp index 0e72635..e059941 100644 --- a/src/Util.cpp +++ b/src/Util.cpp @@ -115,7 +115,7 @@ namespace da::sha { } void - WriteTensorToVtk(const boost::filesystem::path &file_path, const Tensor3d &t3, std::shared_ptr sp_mesh) { + WriteTensorToVtk(const boost::filesystem::path &file_path, const Tensor3d &t3, std::shared_ptr sp_mesh,double threshold) { HexahedralMatMesh matmesh; Eigen::Vector3d origin = sp_mesh->GetOrigin(); double pixel_len = sp_mesh->GetPixelLen(); @@ -129,6 +129,18 @@ namespace da::sha { v_ele_data.push_back(*(t3.data() + pixel_idx[i])); } assert(matmesh.IsHexahedral()); + if(threshold!=-1){ + if(threshold<0 || threshold>1){ + throw std::runtime_error("threshold should be in [0,1]"); + } + for(auto &d:v_ele_data){ + if(d>=threshold){ + d=1; + }else{ + d=0; + } + } + } WriteHexahedralMatmeshToVtk(file_path, matmesh, std::vector(), v_ele_data); } diff --git a/src/Util.h b/src/Util.h index 54e30bd..ffc7a43 100644 --- a/src/Util.h +++ b/src/Util.h @@ -183,7 +183,7 @@ namespace da::sha { std::shared_ptr sp_mesh); void WriteTensorToVtk(const fs_path &file_path, const Tensor3d &t3, - std::shared_ptr sp_mesh); + std::shared_ptr sp_mesh,double threshold=-1); inline void clear_dir(const std::filesystem::path &path) { for (const auto &entry: std::filesystem::recursive_directory_iterator(path)) {