cflin
/
Thermo-elasticTopologyOptimization
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

todo: fix Emin

multiple_top
cflin 2 years ago
parent
1e48f868c1
commit
6e84519dc7
6 changed files with 623 additions and 608 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 2
      examples/top-mech-regular-model-test/config.json
  2. 18
      examples/top-thermoelastic-regular-model/config.json
  3. 378
      src/ThermoelasticTop3d.cpp
  4. 381
      src/ThermoelasticTop3d.h
  5. 451
      src/Top3d.cpp
  6. 1
      src/Top3d.h

2
examples/top-mech-regular-model-test/config.json
View File

@ -6,7 +6,7 @@
"topology": { "topology": {
"max_loop": 100, "max_loop": 100,
"volfrac": 0.5, "volfrac": 0.5,
"penal": 3.0, "penal": 1.0,
"r_min": 1.5 "r_min": 1.5
}, },
"model": { "model": {

18
examples/top-thermoelastic-regular-model/config.json
View File

@ -2,23 +2,23 @@
"material": { "material": {
"E": 1.0, "E": 1.0,
"poisson_ratio": 0.3, "poisson_ratio": 0.3,
"thermal_conductivity":1e11, "unit": "W/K", "thermal_conductivity":0.004, "unit": "W/K",
"thermal_expansion_coefficient": 1e-15, "unit": "1/K" "thermal_expansion_coefficient": 1e-6, "unit": "1/K"
}, },
"topology": { "topology": {
"max_loop": 100, "max_loop": 100,
"volfrac": 0.5, "volfrac": 0.5,
"r_min": 1.5, "r_min": 1.5,
"T_ref": 295, "T_ref": 295,
"T_limit": 3250, "T_limit": 325,
"R_E": 8, "R_E": 0,
"R_lambda": 0, "R_lambda": 0,
"R_beta": 0 "R_beta": 0
}, },
"model": { "model": {
"regular_model": { "regular_model": {
"lx": 1, "lx": 2,
"ly": 40, "ly": 20,
"lz": 30 "lz": 30
}, },
"visual_model":"2-box-refined.obj" "visual_model":"2-box-refined.obj"
@ -28,7 +28,7 @@
{ {
"min": [0, 1, 0], "min": [0, 1, 0],
"max": [1, 1, 0], "max": [1, 1, 0],
"val": [0.0, 0.0, -500.0] "val": [0.0, 0.0, -1.0]
} }
], ],
@ -45,12 +45,12 @@
{ {
"min": [0, 0.25, 0.5], "min": [0, 0.25, 0.5],
"max": [1, 0.25, 0.5], "max": [1, 0.25, 0.5],
"heat_flux": 0.06, "unit": "W" "heat_flux": 1e-15, "unit": "W"
}, },
{ {
"min": [0, 0.75, 0.5], "min": [0, 0.75, 0.5],
"max": [1, 0.75, 0.5], "max": [1, 0.75, 0.5],
"heat_flux": 0.06, "unit": "W" "heat_flux": 1e-15, "unit": "W"
} }
], ],

378
src/ThermoelasticTop3d.cpp
View File

@ -3,3 +3,381 @@
// //
#include "ThermoelasticTop3d.h" #include "ThermoelasticTop3d.h"
da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::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);
// }
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
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<int, std::vector<LimitedDof>> map_ele2Limit;
std::vector<int> 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 && 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<Eigen::SparseMatrix<double>> solver_th;
//#else
// Eigen::SimplicialLLT<Eigen::SparseMatrix<double>> 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_ *
(1e-6+
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<Eigen::SparseMatrix<double>> solver;
#else
Eigen::SimplicialLLT<Eigen::SparseMatrix<double>> 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<double> 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<double> 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<double> &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 = CalDlambdaDrho_Vec(xPhys_col).array() * ce_th.array();
// dc_drho
// Eigen::VectorXd dc_drho = lambda_t_Mul_dKt_drho_Mul_T +
// lambda_m_Mul_dKm_drho_Mul_U +
// 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 = sp_mech_top3d_->drho_dx_ * dc_drho;
// // dT_dx
// Eigen::MatrixXd dT_dx = sp_mech_top3d_->drho_dx_ * dT_drho;
// mma solver
size_t num_constraints =
1 ;// volume and temperature constraints
size_t num_variables = flg_chosen ? chosen_ele_id.size() : sp_mesh_->GetNumEles();
auto mma = std::make_shared<MMASolver>(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 fval = (Eigen::VectorXd(num_constraints) << (v / num_variables - sp_para_->volfrac)).finished();
Eigen::VectorXd fval = (Eigen::VectorXd(num_constraints) << (v - num_variables *sp_para_->volfrac)).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).finished().transpose();
Eigen::MatrixXd dfdx = (Eigen::MatrixXd(num_variables, num_constraints)
<< dv).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 << fval.transpose() << std::endl;
#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<Eigen::DenseIndex, 3>{
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<Eigen::DenseIndex, 3>{
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<Eigen::DenseIndex, 3>{
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_;
}

381
src/ThermoelasticTop3d.h
View File

@ -15,378 +15,7 @@ namespace da::sha::top {
Precompute(); Precompute();
} }
Tensor3d TopOptMainLoop() { 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);
// }
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
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<int, std::vector<LimitedDof>> map_ele2Limit;
std::vector<int> 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<Eigen::SparseMatrix<double>> solver_th;
#else
Eigen::SimplicialLLT<Eigen::SparseMatrix<double>> 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<Eigen::SparseMatrix<double>> solver;
#else
Eigen::SimplicialLLT<Eigen::SparseMatrix<double>> 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<double> 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<double> 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<double> &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<MMASolver>(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<<fval.transpose()<<std::endl;
#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<Eigen::DenseIndex, 3>{
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<Eigen::DenseIndex, 3>{
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<Eigen::DenseIndex, 3>{
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 { Eigen::VectorXd GetU() const {
return sp_thermal_top3d_->GetU(); return sp_thermal_top3d_->GetU();
@ -414,13 +43,6 @@ namespace da::sha::top {
D0_ = sp_mech_top3d_->sp_fea_->computeD(1.0); D0_ = sp_mech_top3d_->sp_fea_->computeD(1.0);
const Eigen::MatrixXd Be = sp_mech_top3d_->sp_fea_->computeBe(); const Eigen::MatrixXd Be = sp_mech_top3d_->sp_fea_->computeBe();
Inted_ = Be.transpose() * D0_ * (Eigen::VectorXd(6) << 1, 1, 1, 0, 0, 0).finished(); 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<double> 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_;
} }
@ -431,7 +53,6 @@ namespace da::sha::top {
Eigen::VectorXi i_dFth_drho_, j_dFth_drho_; Eigen::VectorXi i_dFth_drho_, j_dFth_drho_;
Eigen::VectorXd Inted_; Eigen::VectorXd Inted_;
Eigen::MatrixXd D0_; Eigen::MatrixXd D0_;
Eigen::SparseMatrix<double> drho_dx_;
}; };
} }

451
src/Top3d.cpp
View File

@ -7,235 +7,250 @@
#include <cassert> #include <cassert>
#include "Eigen/src/Core/Matrix.h" #include "Eigen/src/Core/Matrix.h"
#include "Util.h" #include "Util.h"
namespace da::sha { namespace da::sha {
namespace top { namespace top {
Tensor3d Top3d::TopOptMainLoop() { Tensor3d Top3d::TopOptMainLoop() {
Eigen::VectorXd xPhys_col(sp_mesh_->GetNumEles()); Eigen::VectorXd xPhys_col(sp_mesh_->GetNumEles());
Eigen::VectorXi chosen_ele_id(sp_mesh_->GetChosenEleIdx()); Eigen::VectorXi chosen_ele_id(sp_mesh_->GetChosenEleIdx());
bool flg_chosen = chosen_ele_id.size() != 0; bool flg_chosen = chosen_ele_id.size() != 0;
if (!flg_chosen) { if (!flg_chosen) {
// no part chosen // no part chosen
xPhys_col.setConstant(sp_para_->volfrac); xPhys_col.setConstant(sp_para_->volfrac);
} else { } else {
// pick chosen part to sim // pick chosen part to sim
xPhys_col = sp_mesh_->GetInitEleRho(); xPhys_col = sp_mesh_->GetInitEleRho();
xPhys_col(chosen_ele_id).setConstant(sp_para_->volfrac); xPhys_col(chosen_ele_id).setConstant(sp_para_->volfrac);
} }
int loop = 0; int loop = 0;
double change = 1.0; double change = 1.0;
double E0 = sp_fea_->sp_material_->E; double E0 = sp_fea_->sp_material_->E;
double Emin = E0 * sp_para_->E_factor; double Emin = E0 * sp_para_->E_factor;
// precompute // precompute
Eigen::VectorXd dv(sp_mesh_->GetNumEles()); Eigen::VectorXd dv(sp_mesh_->GetNumEles());
dv.setOnes(); dv.setOnes();
dv = H_ * (dv.array() / Hs_.array()).matrix().eval(); dv = H_ * (dv.array() / Hs_.array()).matrix().eval();
Eigen::VectorXd ele_to_write = 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 pixel_idx = sp_mesh_->GetPixelIdx(); Eigen::VectorXi pixel_idx = sp_mesh_->GetPixelIdx();
spdlog::info("end precompute"); spdlog::info("end precompute");
// // clear output dir // // clear output dir
// clear_dir(CMAKE_SOURCE_DIR "/output"); // clear_dir(CMAKE_SOURCE_DIR "/output");
#ifdef USE_SUITESPARSE #ifdef USE_SUITESPARSE
spdlog::info("using suitesparse solver"); spdlog::info("using suitesparse solver");
#else #else
spdlog::warn("using Eigen built-in direct solver!"); spdlog::warn("using Eigen built-in direct solver!");
#endif #endif
// start iteration // start iteration
while (change > sp_para_->tol_x && loop < sp_para_->max_loop) { while (change > sp_para_->tol_x && loop < sp_para_->max_loop) {
++loop; ++loop;
Eigen::VectorXd sK = // filter
(sKe_ * (Emin + xPhys_col.array().pow(sp_para_->penal) * (E0 - Emin)).matrix().transpose()) xPhys_col = H_ * (xPhys_col.array() / Hs_.array()).matrix().eval();
.reshaped(); Eigen::VectorXd sK =
auto v_tri = Vec2Triplet(iK_, jK_, sK); (sKe_ * (Emin + xPhys_col.array().pow(sp_para_->penal) * (E0 - Emin)).matrix().transpose())
K_.setFromTriplets(v_tri.begin(), v_tri.end()); .reshaped();
IntroduceFixedDofs(K_, F_); auto v_tri = Vec2Triplet(iK_, jK_, sK);
K_.setFromTriplets(v_tri.begin(), v_tri.end());
IntroduceFixedDofs(K_, F_);
#ifdef USE_SUITESPARSE #ifdef USE_SUITESPARSE
Eigen::CholmodSupernodalLLT<Eigen::SparseMatrix<double>> solver; Eigen::CholmodSupernodalLLT<Eigen::SparseMatrix<double>> solver;
#else #else
Eigen::SimplicialLLT<Eigen::SparseMatrix<double>> solver; Eigen::SimplicialLLT<Eigen::SparseMatrix<double>> solver;
#endif #endif
solver.compute(K_); solver.compute(K_);
U_ = solver.solve(F_); U_ = solver.solve(F_);
// compliance // compliance
Eigen::VectorXd ce(sp_mesh_->GetNumEles()); Eigen::VectorXd ce(sp_mesh_->GetNumEles());
for (int i = 0; i < sp_mesh_->GetNumEles(); ++i) { for (int i = 0; i < sp_mesh_->GetNumEles(); ++i) {
Eigen::VectorXi dofs_in_ele_i = sp_mesh_->MapEleId2Dofs(i); Eigen::VectorXi dofs_in_ele_i = sp_mesh_->MapEleId2Dofs(i);
Eigen::VectorXd Ue = U_(dofs_in_ele_i); Eigen::VectorXd Ue = U_(dofs_in_ele_i);
ce(i) = Ue.transpose() * Ke_ * Ue; ce(i) = Ue.transpose() * Ke_ * Ue;
} }
double c = double c =
ce.transpose() * (Emin + xPhys_col.array().pow(sp_para_->penal) * (E0 - Emin)).matrix(); ce.transpose() * (Emin + xPhys_col.array().pow(sp_para_->penal) * (E0 - Emin)).matrix();
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();
Eigen::VectorXd dc = // sensitivity
-sp_para_->penal * (E0 - Emin) * xPhys_col.array().pow(sp_para_->penal - 1.0) * ce.array(); Eigen::VectorXd dc_drho =
-sp_para_->penal * (E0 - Emin) * xPhys_col.array().pow(sp_para_->penal - 1.0) * ce.array();
// mma solver Eigen::VectorXd dc_dx = drho_dx_ * dc_drho;
size_t num_constrants = 1;
size_t num_variables = flg_chosen ? chosen_ele_id.size() : sp_mesh_->GetNumEles(); // mma solver
auto mma = std::make_shared<MMASolver>(num_variables, num_constrants); size_t num_constrants = 1;
Eigen::VectorXd variables_tmp = flg_chosen ? xPhys_col(chosen_ele_id) : xPhys_col; size_t num_variables = flg_chosen ? chosen_ele_id.size() : sp_mesh_->GetNumEles();
double f0val = c; auto mma = std::make_shared<MMASolver>(num_variables, num_constrants);
Eigen::VectorXd df0dx = flg_chosen Eigen::VectorXd variables_tmp = flg_chosen ? xPhys_col(chosen_ele_id) : xPhys_col;
? dc(chosen_ele_id).eval() / dc(chosen_ele_id).cwiseAbs().maxCoeff() double f0val = c;
: dc / dc.cwiseAbs().maxCoeff(); Eigen::VectorXd df0dx = flg_chosen
double fval = v - num_variables * sp_para_->volfrac; ? dc_dx(chosen_ele_id).eval() / dc_dx(chosen_ele_id).cwiseAbs().maxCoeff()
Eigen::VectorXd dfdx = flg_chosen ? dv(chosen_ele_id) : dv; : dc_dx / dc_dx.cwiseAbs().maxCoeff();
static Eigen::VectorXd low_bounds = Eigen::VectorXd::Zero(num_variables); double fval = v - num_variables * sp_para_->volfrac;
static Eigen::VectorXd up_bounds = Eigen::VectorXd::Ones(num_variables); Eigen::VectorXd dfdx = flg_chosen ? dv(chosen_ele_id) : dv;
static Eigen::VectorXd low_bounds = Eigen::VectorXd::Zero(num_variables);
// spdlog::info("mma update"); static Eigen::VectorXd up_bounds = Eigen::VectorXd::Ones(num_variables);
mma->Update(variables_tmp.data(), df0dx.data(), &fval, dfdx.data(), low_bounds.data(),
up_bounds.data()); // spdlog::info("mma update");
if (flg_chosen) { mma->Update(variables_tmp.data(), df0dx.data(), &fval, dfdx.data(), low_bounds.data(),
change = (variables_tmp - xPhys_col(chosen_ele_id)).cwiseAbs().maxCoeff(); up_bounds.data());
xPhys_col(chosen_ele_id) = variables_tmp; if (flg_chosen) {
} else { change = (variables_tmp - xPhys_col(chosen_ele_id)).cwiseAbs().maxCoeff();
change = (variables_tmp - xPhys_col).cwiseAbs().maxCoeff(); xPhys_col(chosen_ele_id) = variables_tmp;
xPhys_col = 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); }
spdlog::critical("Iter: {:3d}, Comp: {:.3e}, Vol: {:.2f}, Change: {:f}", loop, c, v, change);
#ifdef WRITE_TENSOR_IN_LOOP #ifdef WRITE_TENSOR_IN_LOOP
// extract vtk // extract vtk
ele_to_write(pixel_idx) = xPhys_col; ele_to_write(pixel_idx) = 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) { 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(i, 0, 0) = ele_to_write(i);
} }
ten_xPhys_to_write = ten_xPhys_to_write.reshape(Eigen::array<Eigen::DenseIndex, 3>{ ten_xPhys_to_write = ten_xPhys_to_write.reshape(Eigen::array<Eigen::DenseIndex, 3>{
sp_mesh_->GetLx(), sp_mesh_->GetLy(), sp_mesh_->GetLz()}); sp_mesh_->GetLx(), sp_mesh_->GetLy(), sp_mesh_->GetLz()});
top::WriteTensorToVtk( top::WriteTensorToVtk(
da::WorkingResultDirectoryPath() / ("field_matrix" + std::to_string(loop) + ".vtk"), da::WorkingResultDirectoryPath() / ("field_matrix" + std::to_string(loop) + ".vtk"),
ten_xPhys_to_write, sp_mesh_); ten_xPhys_to_write, sp_mesh_);
#endif #endif
} }
// result // result
rho_ = xPhys_col; rho_ = xPhys_col;
// set 0 to rho of unchosen part // set 0 to rho of unchosen part
assert(xPhys_col.size()); assert(xPhys_col.size());
Eigen::VectorXi continue_idx = Eigen::VectorXi continue_idx =
Eigen::VectorXi::LinSpaced(xPhys_col.size(), 0, xPhys_col.size() - 1); Eigen::VectorXi::LinSpaced(xPhys_col.size(), 0, xPhys_col.size() - 1);
Eigen::VectorXi unchosen_idx =flg_chosen? SetDifference(continue_idx, chosen_ele_id): Eigen::VectorXi(); Eigen::VectorXi unchosen_idx = flg_chosen ? SetDifference(continue_idx, chosen_ele_id) : Eigen::VectorXi();
{ {
xPhys_col(unchosen_idx).setZero(); xPhys_col(unchosen_idx).setZero();
ele_to_write(pixel_idx) = xPhys_col; ele_to_write(pixel_idx) = 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) { 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(i, 0, 0) = ele_to_write(i);
} }
ten_xPhys_to_write = ten_xPhys_to_write.reshape(Eigen::array<Eigen::DenseIndex, 3>{ ten_xPhys_to_write = ten_xPhys_to_write.reshape(Eigen::array<Eigen::DenseIndex, 3>{
sp_mesh_->GetLx(), sp_mesh_->GetLy(), sp_mesh_->GetLz()}); sp_mesh_->GetLx(), sp_mesh_->GetLy(), sp_mesh_->GetLz()});
rho_field_zero_filled_ = ten_xPhys_to_write; rho_field_zero_filled_ = ten_xPhys_to_write;
} }
{ {
xPhys_col(unchosen_idx).setOnes(); xPhys_col(unchosen_idx).setOnes();
ele_to_write(pixel_idx) = xPhys_col; ele_to_write(pixel_idx) = 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) { 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(i, 0, 0) = ele_to_write(i);
} }
ten_xPhys_to_write = ten_xPhys_to_write.reshape(Eigen::array<Eigen::DenseIndex, 3>{ ten_xPhys_to_write = ten_xPhys_to_write.reshape(Eigen::array<Eigen::DenseIndex, 3>{
sp_mesh_->GetLx(), sp_mesh_->GetLy(), sp_mesh_->GetLz()}); sp_mesh_->GetLx(), sp_mesh_->GetLy(), sp_mesh_->GetLz()});
rho_field_one_filled_ = ten_xPhys_to_write; rho_field_one_filled_ = ten_xPhys_to_write;
} }
return rho_field_zero_filled_; return rho_field_zero_filled_;
}
std::vector<Tensor3d> Top3d::GetTensorOfStress(const Eigen::VectorXd &which_col_of_stress) {
Eigen::VectorXd ele_to_write =
Eigen::VectorXd::Zero(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz());
Eigen::VectorXi pixel_idx = sp_mesh_->GetPixelIdx();
// stress
Eigen::MatrixXd mat_stress(sp_mesh_->GetNumEles(), 6);
Eigen::MatrixXd B=sp_fea_->computeBe({0, 0, 0});
for (int i = 0; i < sp_mesh_->GetNumEles(); ++i) {
Eigen::VectorXi dofs_in_ele_i = sp_mesh_->MapEleId2Dofs(i);
Eigen::VectorXd Ue = U_(dofs_in_ele_i);
mat_stress.row(i) = rho_(i) * sp_fea_->computeD() * B * Ue;
}
// fill
std::vector<Tensor3d> vt;
for (int i = 0; i < which_col_of_stress.size(); ++i) {
ele_to_write(pixel_idx) = mat_stress.col(which_col_of_stress(i));
vt.push_back(GetTensorFromCol(ele_to_write));
}
return vt;
}
Tensor3d Top3d::GetTensorFromCol(const Eigen::VectorXd &proprty_col) {
Tensor3d ten_prop_to_write(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz(), 1, 1);
assert(proprty_col.size() == ten_prop_to_write.size());
for (int i = 0; i < proprty_col.size(); ++i) {
ten_prop_to_write(i, 0, 0) = proprty_col(i);
}
ten_prop_to_write = ten_prop_to_write.reshape(
Eigen::array<Eigen::DenseIndex, 3>{sp_mesh_->GetLx(), sp_mesh_->GetLy(), sp_mesh_->GetLz()});
return ten_prop_to_write;
}
void Top3d::Precompute() {
Eigen::MatrixXi mat_ele2dofs = sp_mesh_->GetEleId2DofsMap();
int dofs_each_ele = sp_mesh_->Get_NUM_NODES_EACH_ELE() *
sp_mesh_->Get_DOFS_EACH_NODE(); // 24 for mathe; 8 for heat
iK_ = Eigen::KroneckerProduct(mat_ele2dofs, Eigen::VectorXi::Ones(dofs_each_ele))
.transpose()
.reshaped();
jK_ = Eigen::KroneckerProduct(mat_ele2dofs, Eigen::RowVectorXi::Ones(dofs_each_ele))
.transpose()
.reshaped();
Ke_=sp_fea_->computeKe(1.0);
sKe_ = Ke_.reshaped();
// precompute filter
Eigen::VectorXi iH = Eigen::VectorXi::Ones(
sp_mesh_->GetNumEles() * std::pow(2.0 * (std::ceil(sp_para_->r_min) - 1.0) + 1, 3));
Eigen::VectorXi jH = iH;
Eigen::VectorXd sH(iH.size());
sH.setZero();
int cnt = 0;
Hs_ = Eigen::VectorXd(sp_mesh_->GetNumEles());
int delta = std::ceil(sp_para_->r_min) - 1;
for (int k = 0; k < sp_mesh_->GetLz(); ++k) {
for (int j = 0; j < sp_mesh_->GetLy(); ++j) {
for (int i = 0; i < sp_mesh_->GetLx(); ++i) {
int ele_id0 = sp_mesh_->MapEleCoord2Id((Eigen::MatrixXi(1, 3) << i, j, k).finished())(0);
if (ele_id0 == -1) {
continue;
} }
for (int k2 = std::max(k - delta, 0); k2 <= std::min(k + delta, sp_mesh_->GetLz() - 1);
++k2) { std::vector<Tensor3d> Top3d::GetTensorOfStress(const Eigen::VectorXd &which_col_of_stress) {
for (int j2 = std::max(j - delta, 0); j2 <= std::min(j + delta, sp_mesh_->GetLy() - 1); Eigen::VectorXd ele_to_write =
++j2) { Eigen::VectorXd::Zero(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz());
for (int i2 = std::max(i - delta, 0); i2 <= std::min(i + delta, sp_mesh_->GetLx() - 1); Eigen::VectorXi pixel_idx = sp_mesh_->GetPixelIdx();
++i2) { // stress
int ele_id1 = Eigen::MatrixXd mat_stress(sp_mesh_->GetNumEles(), 6);
sp_mesh_->MapEleCoord2Id((Eigen::MatrixXi(1, 3) << i2, j2, k2).finished())(0); Eigen::MatrixXd B = sp_fea_->computeBe({0, 0, 0});
if (ele_id1 == -1) { for (int i = 0; i < sp_mesh_->GetNumEles(); ++i) {
continue; Eigen::VectorXi dofs_in_ele_i = sp_mesh_->MapEleId2Dofs(i);
} Eigen::VectorXd Ue = U_(dofs_in_ele_i);
iH(cnt) = ele_id0; mat_stress.row(i) = rho_(i) * sp_fea_->computeD() * B * Ue;
jH(cnt) = ele_id1;
sH(cnt) =
std::max(0.0, sp_para_->r_min - Eigen::Vector3d(i - i2, j - j2, k - k2).norm());
Hs_(ele_id0) += sH(cnt);
++cnt;
} }
} // fill
std::vector<Tensor3d> vt;
for (int i = 0; i < which_col_of_stress.size(); ++i) {
ele_to_write(pixel_idx) = mat_stress.col(which_col_of_stress(i));
vt.push_back(GetTensorFromCol(ele_to_write));
}
return vt;
}
Tensor3d Top3d::GetTensorFromCol(const Eigen::VectorXd &proprty_col) {
Tensor3d ten_prop_to_write(sp_mesh_->GetLx() * sp_mesh_->GetLy() * sp_mesh_->GetLz(), 1, 1);
assert(proprty_col.size() == ten_prop_to_write.size());
for (int i = 0; i < proprty_col.size(); ++i) {
ten_prop_to_write(i, 0, 0) = proprty_col(i);
}
ten_prop_to_write = ten_prop_to_write.reshape(
Eigen::array<Eigen::DenseIndex, 3>{sp_mesh_->GetLx(), sp_mesh_->GetLy(), sp_mesh_->GetLz()});
return ten_prop_to_write;
}
void Top3d::Precompute() {
Eigen::MatrixXi mat_ele2dofs = sp_mesh_->GetEleId2DofsMap();
int dofs_each_ele = sp_mesh_->Get_NUM_NODES_EACH_ELE() *
sp_mesh_->Get_DOFS_EACH_NODE(); // 24 for mathe; 8 for heat
iK_ = Eigen::KroneckerProduct(mat_ele2dofs, Eigen::VectorXi::Ones(dofs_each_ele))
.transpose()
.reshaped();
jK_ = Eigen::KroneckerProduct(mat_ele2dofs, Eigen::RowVectorXi::Ones(dofs_each_ele))
.transpose()
.reshaped();
Ke_ = sp_fea_->computeKe(1.0);
sKe_ = Ke_.reshaped();
// precompute filter
Eigen::VectorXi iH = Eigen::VectorXi::Ones(
sp_mesh_->GetNumEles() * std::pow(2.0 * (std::ceil(sp_para_->r_min) - 1.0) + 1, 3));
Eigen::VectorXi jH = iH;
Eigen::VectorXd sH(iH.size());
sH.setZero();
int cnt = 0;
Hs_ = Eigen::VectorXd(sp_mesh_->GetNumEles());
int delta = std::ceil(sp_para_->r_min) - 1;
for (int k = 0; k < sp_mesh_->GetLz(); ++k) {
for (int j = 0; j < sp_mesh_->GetLy(); ++j) {
for (int i = 0; i < sp_mesh_->GetLx(); ++i) {
int ele_id0 = sp_mesh_->MapEleCoord2Id((Eigen::MatrixXi(1, 3) << i, j, k).finished())(0);
if (ele_id0 == -1) {
continue;
}
for (int k2 = std::max(k - delta, 0); k2 <= std::min(k + delta, sp_mesh_->GetLz() - 1);
++k2) {
for (int j2 = std::max(j - delta, 0); j2 <= std::min(j + delta, sp_mesh_->GetLy() - 1);
++j2) {
for (int i2 = std::max(i - delta, 0); i2 <= std::min(i + delta, sp_mesh_->GetLx() - 1);
++i2) {
int ele_id1 =
sp_mesh_->MapEleCoord2Id((Eigen::MatrixXi(1, 3) << i2, j2, k2).finished())(
0);
if (ele_id1 == -1) {
continue;
}
iH(cnt) = ele_id0;
jH(cnt) = ele_id1;
sH(cnt) =
std::max(0.0,
sp_para_->r_min - Eigen::Vector3d(i - i2, j - j2, k - k2).norm());
Hs_(ele_id0) += sH(cnt);
++cnt;
}
}
}
}
}
}
std::vector<Eigen::Triplet<double>> v_tri = Vec2Triplet(iH, jH, sH);
H_ = SpMat(sp_mesh_->GetNumEles(), sp_mesh_->GetNumEles());
H_.setFromTriplets(v_tri.begin(), v_tri.end());
Eigen::VectorXi i_Hs = Eigen::VectorXi::LinSpaced(Hs_.size(), 0,
Hs_.size());
Eigen::SparseMatrix<double> sp_inv_Hs(i_Hs.size(), i_Hs.size());
auto v_inv_Hs_tri = Vec2Triplet(i_Hs, i_Hs, Eigen::VectorXd(1.0 / Hs_.array()));
sp_inv_Hs.setFromTriplets(v_inv_Hs_tri.begin(), v_inv_Hs_tri.end());
drho_dx_ = sp_inv_Hs * H_;
} }
} } // namespace top
}
}
std::vector<Eigen::Triplet<double>> v_tri = Vec2Triplet(iH, jH, sH);
H_ = SpMat(sp_mesh_->GetNumEles(), sp_mesh_->GetNumEles());
H_.setFromTriplets(v_tri.begin(), v_tri.end());
}
} // namespace top
} // namespace da::sha } // namespace da::sha

1
src/Top3d.h
View File

@ -188,6 +188,7 @@ namespace da::sha {
SpMat H_; SpMat H_;
Eigen::VectorXd Hs_; Eigen::VectorXd Hs_;
Eigen::SparseMatrix<double> drho_dx_;
// result // result
Eigen::VectorXd U_; Eigen::VectorXd U_;
Eigen::VectorXd rho_; Eigen::VectorXd rho_;

Write Preview
Loading…
Cancel
Save