From a83bbddf036f88a560bd397d6fe5b60bbb41212d Mon Sep 17 00:00:00 2001 From: Manav Bhatia Date: Wed, 9 Oct 2019 22:02:41 -0500 Subject: [PATCH] bug-fixes for sensitivity analysis of stress evaluation. --- ...ermal_stress_jacobian_scaling_function.cpp | 2 +- src/elasticity/structural_element_1d.cpp | 19 +++++--- src/elasticity/structural_element_2d.cpp | 43 ++++++++++++++----- 3 files changed, 46 insertions(+), 18 deletions(-) diff --git a/examples/structural/base/thermal_stress_jacobian_scaling_function.cpp b/examples/structural/base/thermal_stress_jacobian_scaling_function.cpp index fdefb9af..857643bd 100644 --- a/examples/structural/base/thermal_stress_jacobian_scaling_function.cpp +++ b/examples/structural/base/thermal_stress_jacobian_scaling_function.cpp @@ -89,7 +89,7 @@ MAST::Examples::ThermalJacobianScaling::operator() (Real& val) const { if (max_res > 0.) max_log = std::log10(max_res); - val = std::fmax(1.-std::exp(_accel_factor*(log-max_log)), low); + val = std::fmax(std::pow(1-res/max_res, _accel_factor), low); //libMesh::out << log << " " << max_log << " " << val << std::endl; } } diff --git a/src/elasticity/structural_element_1d.cpp b/src/elasticity/structural_element_1d.cpp index 9e647405..a790b537 100644 --- a/src/elasticity/structural_element_1d.cpp +++ b/src/elasticity/structural_element_1d.cpp @@ -247,7 +247,9 @@ MAST::StructuralElement1D::calculate_stress(bool request_derivative, strain_3D = RealVectorX::Zero(6), stress_3D = RealVectorX::Zero(6), dstrain_dp = RealVectorX::Zero(n1), - dstress_dp = RealVectorX::Zero(n1); + dstress_dp = RealVectorX::Zero(n1), + vec1 = RealVectorX::Zero(n2), + vec2 = RealVectorX::Zero(n2); MAST::FEMOperatorMatrix Bmat_mem, @@ -433,8 +435,12 @@ MAST::StructuralElement1D::calculate_stress(bool request_derivative, dstress_dX = material_mat * dstrain_dX; // copy to the 3D structure - dstress_dX_3D.row(0) = dstress_dX.row(0); - dstrain_dX_3D.row(0) = dstrain_dX.row(0); + vec1 = dstress_dX.row(0); + this->transform_vector_to_global_system(vec1, vec2); + dstress_dX_3D.row(0) = vec2; + vec1 = dstrain_dX.row(0); + this->transform_vector_to_global_system(vec1, vec2); + dstrain_dX_3D.row(0) = vec2; if (request_derivative) data->set_derivatives(dstress_dX_3D, dstrain_dX_3D); @@ -461,7 +467,7 @@ MAST::StructuralElement1D::calculate_stress(bool request_derivative, temp_func->derivative(*p, xyz[qp_loc_index], _time, dtemp); ref_temp_func->derivative(*p, xyz[qp_loc_index], _time, dref_t); alpha_func->derivative(*p, xyz[qp_loc_index], _time, dalpha); - dstrain_dp(0) -= alpha*(dtemp-dref_t) - dalpha*(temp-ref_t); + dstrain_dp(0) -= alpha*(dtemp-dref_t) + dalpha*(temp-ref_t); } // include the dependence of strain on the thickness @@ -2098,7 +2104,7 @@ MAST::StructuralElement1D::thermal_residual_sensitivity (const MAST::FunctionBas &temp_func = bc.get >("temperature"), &ref_temp_func = bc.get >("ref_temperature"); - Real t, t0, t_sens; + Real t, t0, t_sens, t0_sens; for (unsigned int qp=0; qptransform_vector_to_global_system(vec1, vec2); + dstress_dX_3D.row(0) = vec2; + // sigma-yy + vec1 = dstress_dX.row(1); + this->transform_vector_to_global_system(vec1, vec2); + dstress_dX_3D.row(1) = vec2; + // tau-xy + vec1 = dstress_dX.row(2); + this->transform_vector_to_global_system(vec1, vec2); + dstress_dX_3D.row(3) = vec2; + // epsilon-xx + vec1 = dstrain_dX.row(0); + this->transform_vector_to_global_system(vec1, vec2); + dstrain_dX_3D.row(0) = vec2; + // epsilon-yy + vec1 = dstrain_dX.row(1); + this->transform_vector_to_global_system(vec1, vec2); + dstrain_dX_3D.row(1) = vec2; + // gamma-xy + vec1 = dstrain_dX.row(2); + this->transform_vector_to_global_system(vec1, vec2); + dstrain_dX_3D.row(3) = vec2; if (request_derivative) data->set_derivatives(dstress_dX_3D, dstrain_dX_3D); @@ -559,8 +579,8 @@ MAST::StructuralElement2D::calculate_stress(bool request_derivative, temp_func->derivative(*p, xyz[qp_loc_index], _time, dtemp); ref_temp_func->derivative(*p, xyz[qp_loc_index], _time, dref_t); alpha_func->derivative(*p, xyz[qp_loc_index], _time, dalpha); - dstrain_dp(0) -= alpha*(dtemp-dref_t) - dalpha*(temp-ref_t); // epsilon-xx - dstrain_dp(1) -= alpha*(dtemp-dref_t) - dalpha*(temp-ref_t); // epsilon-yy + dstrain_dp(0) -= alpha*(dtemp-dref_t) + dalpha*(temp-ref_t); // epsilon-xx + dstrain_dp(1) -= alpha*(dtemp-dref_t) + dalpha*(temp-ref_t); // epsilon-yy } @@ -2715,7 +2735,7 @@ thermal_residual_sensitivity (const MAST::FunctionBase& p, &temp_func = bc.get >("temperature"), &ref_temp_func = bc.get >("ref_temperature"); - Real t, t0, t_sens; + Real t, t0, t_sens, t0_sens; for (unsigned int qp=0; qp