| 
						
						
							
								
							
						
						
					 | 
				
				 | 
				
					@ -8,7 +8,7 @@ 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); | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					    xPhys_col.setConstant(sp_para_->volfrac/sp_para_->volfrac); | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					    bool flg_chosen = false; | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					    Eigen::VectorXi chosen_ele_id; | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					//            Eigen::VectorXi chosen_ele_id(sp_mesh_->GetChosenEleIdx());
 | 
				
			
			
		
	
	
		
			
				
					| 
						
							
								
							
						
						
							
								
							
						
						
					 | 
				
				 | 
				
					@ -118,7 +118,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					            return CalDRDrho(rho, sp_para_->R_beta) * E0_m * alpha0; | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        }; | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#ifndef MECH_ONLY | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        // solve T
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        Eigen::VectorXd sK_th = | 
				
			
			
		
	
	
		
			
				
					| 
						
						
						
							
								
							
						
					 | 
				
				 | 
				
					@ -134,7 +134,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#endif | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        solver_th.compute(sp_thermal_top3d_->K_); | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        sp_thermal_top3d_->U_ = solver_th.solve(sp_thermal_top3d_->F_); | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#endif | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        // solve U
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					
 | 
				
			
			
		
	
	
		
			
				
					| 
						
						
						
							
								
							
						
					 | 
				
				 | 
				
					@ -143,7 +143,8 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                        .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
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#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()); | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					
 | 
				
			
			
		
	
	
		
			
				
					| 
						
						
						
							
								
							
						
					 | 
				
				 | 
				
					@ -156,7 +157,13 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        } | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        if (loop == 1 || loop == sp_para_->max_loop) | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					            spdlog::info("||Fth|| / ||Fm||: {}", F_th.norm() / sp_mech_top3d_->F_.norm()); | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        Eigen::VectorXd F = Eigen::VectorXd(sp_mech_top3d_->F_) + F_th; | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#endif | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        Eigen::VectorXd F = Eigen::VectorXd(sp_mech_top3d_->F_) | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                            #ifndef MECH_ONLY | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                            + F_th | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#endif | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        ; | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        sp_mech_top3d_->IntroduceFixedDofs(sp_mech_top3d_->K_, F); | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#ifdef USE_SUITESPARSE | 
				
			
			
		
	
	
		
			
				
					| 
						
							
								
							
						
						
							
								
							
						
						
					 | 
				
				 | 
				
					@ -183,7 +190,7 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        // sensitivity
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        // lambda_m
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        Eigen::VectorXd lambda_m = -sp_mech_top3d_->U_; | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#ifndef  MECH_ONLY | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        // dFth_drho
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        Eigen::SparseMatrix<double> dFth_drho(sp_mech_top3d_->sp_mesh_->GetNumEles(), | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                                              sp_mech_top3d_->sp_mesh_->GetNumDofs()); | 
				
			
			
		
	
	
		
			
				
					| 
						
							
								
							
						
						
							
								
							
						
						
					 | 
				
				 | 
				
					@ -228,10 +235,11 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        } | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        // lambda_t
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        Eigen::VectorXd lambda_t = solver_th.solve(rhs); | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#endif | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        // lambda_m_Mul_dKm_drho_Mul_U
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        Eigen::VectorXd lambda_m_Mul_dKm_drho_Mul_U = | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                -CalDEDrho_Vec(xPhys_col).array() * ce.array(); | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#ifndef MECH_ONLY | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        // 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) { | 
				
			
			
		
	
	
		
			
				
					| 
						
						
						
							
								
							
						
					 | 
				
				 | 
				
					@ -245,12 +253,20 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					//        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_);
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#endif | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#ifdef MECH_ONLY | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        Eigen::VectorXd dc_drho = | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                lambda_m_Mul_dKm_drho_Mul_U | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                ; | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#else | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        Eigen::VectorXd dc_drho = | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                lambda_t_Mul_dKt_drho_Mul_T | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                + | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                lambda_m_Mul_dKm_drho_Mul_U | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                + | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                2 * Eigen::VectorXd(dFth_drho * sp_mech_top3d_->U_); | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#endif | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#ifndef MECH_ONLY | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        // dT_drho
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        auto CalDKth_Mul_T__For_DTi_Drhoj = [&](int rho_i) -> Eigen::SparseVector<double> { | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					            Eigen::VectorXi dofs_in_ele = sp_thermal_top3d_->sp_mesh_->MapEleId2Dofs(rho_i); | 
				
			
			
		
	
	
		
			
				
					| 
						
							
								
							
						
						
							
								
							
						
						
					 | 
				
				 | 
				
					@ -283,15 +299,20 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					            } | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        } | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        // 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; | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#endif | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        // dc_dx
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        Eigen::VectorXd dc_dx = sp_mech_top3d_->drho_dx_ * dc_drho; | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        // mma solver
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#define SENSITIVITY_SCALE_COEF 200 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        size_t num_constraints = | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                1 + dT_dx.cols();// volume and temperature constraints
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                1 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                #ifndef MECH_ONLY | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                + dT_dx.cols()// volume and temperature constraints
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#endif | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        ; | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        size_t num_variables = flg_chosen ? chosen_ele_id.size() : sp_mesh_->GetNumEles(); | 
				
			
			
		
	
	
		
			
				
					| 
						
						
						
							
								
							
						
					 | 
				
				 | 
				
					@ -299,15 +320,23 @@ da::sha::top::Tensor3d da::sha::top::ThermoelasticTop3d::TopOptMainLoop() { | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        Eigen::VectorXd variables_tmp = flg_chosen ? xPhys_col(chosen_ele_id) : xPhys_col; | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        double f0val = c; | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        Eigen::VectorXd df0dx = flg_chosen | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                                ? dc_dx(chosen_ele_id).eval() / dc_dx(chosen_ele_id).cwiseAbs().maxCoeff()*SENSITIVITY_SCALE_COEF | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                                ? dc_dx(chosen_ele_id).eval() / dc_dx(chosen_ele_id).cwiseAbs().maxCoeff() * | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                                  SENSITIVITY_SCALE_COEF | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                                : dc_dx / dc_dx.cwiseAbs().maxCoeff() * SENSITIVITY_SCALE_COEF; | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					//                double fval = v - num_variables * sp_para_->volfrac;
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        Eigen::VectorXd fval = (Eigen::VectorXd(num_constraints) << (v / num_variables - sp_para_->volfrac), | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                T(v_dof).array() / sp_para_->T_limit - 1).finished() * SENSITIVITY_SCALE_COEF; | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					        Eigen::VectorXd fval = (Eigen::VectorXd(num_constraints) << (v / num_variables - sp_para_->volfrac) | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#ifndef MECH_ONLY | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                , T(v_dof).array() / sp_para_->T_limit - 1 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#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, 1.0 / sp_para_->T_limit * dT_dx).finished().transpose() * | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                << 1.0 / num_variables * dv | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#ifndef MECH_ONLY | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                , 1.0 / sp_para_->T_limit * dT_dx | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					#endif | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                               ).finished().transpose() * | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					                               SENSITIVITY_SCALE_COEF; | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					
 | 
				
			
			
		
	
		
			
				
					 | 
					 | 
				
				 | 
				
					
 | 
				
			
			
		
	
	
		
			
				
					| 
						
							
								
							
						
						
						
					 | 
				
				 | 
				
					
  |