From 1a324631836534a3559af8c03fb43ace33ace791 Mon Sep 17 00:00:00 2001 From: Sudipta Biswas Date: Tue, 31 Mar 2020 16:18:17 -0600 Subject: [PATCH] Update bondslip constraint jacobian #99 --- include/constraints/RebarBondSlipConstraint.h | 13 +- src/constraints/RebarBondSlipConstraint.C | 244 ++++++++++++------ test/tests/rebar_bondslip/RCBeam_constraint.i | 15 +- 3 files changed, 183 insertions(+), 89 deletions(-) diff --git a/include/constraints/RebarBondSlipConstraint.h b/include/constraints/RebarBondSlipConstraint.h index 6bdab89b..3ac60e8c 100644 --- a/include/constraints/RebarBondSlipConstraint.h +++ b/include/constraints/RebarBondSlipConstraint.h @@ -30,6 +30,15 @@ class RebarBondSlipConstraint : public EqualValueEmbeddedConstraint bool shouldApply() override; void reinitConstraint(); + /** + * Determine whether the coupled variable is one of the displacement variables, + * and find its component + * @param var_num The number of the variable to be checked + * @param component The component index computed in this routine + * @return bool indicating whether the coupled variable is one of the displacement variables + */ + bool getCoupledVarComponent(unsigned int var_num, unsigned int & component); + protected: virtual void computeTangent(); virtual Real computeQpResidual(Moose::ConstraintType type) override; @@ -81,10 +90,12 @@ class RebarBondSlipConstraint : public EqualValueEmbeddedConstraint std::vector _transitional_slip; /// constraint force needed to enforce the constraint RealVectorValue _constraint_residual; + RealVectorValue _constraint_jacobian_axial; /// penalty force for the current constraint RealVectorValue _pen_force; - RealVectorValue _slave_tangent; + RealVectorValue _secondary_tangent; Real _current_elem_volume; bool _bond; Real _bond_stress; + Real _bond_stress_deriv; }; diff --git a/src/constraints/RebarBondSlipConstraint.C b/src/constraints/RebarBondSlipConstraint.C index 7d94e2d2..0180ba26 100644 --- a/src/constraints/RebarBondSlipConstraint.C +++ b/src/constraints/RebarBondSlipConstraint.C @@ -2,7 +2,7 @@ //* https://www.mooseframework.org //* //* All rights reserved, see COPYRIGHT for full restrictions -//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* https://github.com/idaholab/moose/blob/primary/COPYRIGHT //* //* Licensed under LGPL 2.1, please see LICENSE for details //* https://www.gnu.org/licenses/lgpl-2.1.html @@ -73,7 +73,7 @@ RebarBondSlipConstraint::RebarBondSlipConstraint(const InputParameters & paramet void RebarBondSlipConstraint::initialSetup() { - for (auto it = _slave_to_master_map.begin(); it != _slave_to_master_map.end(); ++it) + for (auto it = _secondary_to_primary_map.begin(); it != _secondary_to_primary_map.end(); ++it) if (_bondslip.find(it->first) == _bondslip.end()) _bondslip.insert(std::pair(it->first, bondSlipData())); // initialize @@ -82,7 +82,8 @@ RebarBondSlipConstraint::initialSetup() void RebarBondSlipConstraint::timestepSetup() { - for (auto iter = _slave_to_master_map.begin(); iter != _slave_to_master_map.end(); ++iter) + for (auto iter = _secondary_to_primary_map.begin(); iter != _secondary_to_primary_map.end(); + ++iter) { dof_id_type node_id = iter->first; auto it = _bondslip.find(node_id); @@ -99,22 +100,22 @@ bool RebarBondSlipConstraint::shouldApply() { if (_debug) - // if (_current_node->id() == 144) - { - std::cout << "===========================================\n"; - std::cout << "node id: " << _current_node->id() << std::endl; - std::cout << "at coord: " << (Point)*_current_node << std::endl; - } - auto it = _slave_to_master_map.find(_current_node->id()); + if (_current_node->id() == 144) + { + std::cout << "===========================================\n"; + std::cout << "node id: " << _current_node->id() << std::endl; + std::cout << "at coord: " << (Point)*_current_node << std::endl; + } + auto it = _secondary_to_primary_map.find(_current_node->id()); - if (it != _slave_to_master_map.end()) + if (it != _secondary_to_primary_map.end()) { - const Elem * master_elem = _mesh.elemPtr(it->second); + const Elem * primary_elem = _mesh.elemPtr(it->second); std::vector points = {*_current_node}; - // reinit variables on the master element at the slave point - _fe_problem.setNeighborSubdomainID(master_elem, 0); - _fe_problem.reinitNeighborPhys(master_elem, points, 0); + // reinit variables on the primary element at the secondary point + _fe_problem.setNeighborSubdomainID(primary_elem, 0); + _fe_problem.reinitNeighborPhys(primary_elem, points, 0); reinitConstraint(); @@ -126,7 +127,7 @@ RebarBondSlipConstraint::shouldApply() void RebarBondSlipConstraint::computeTangent() { - _slave_tangent *= 0.0; + _secondary_tangent = 0.0; // get normals // get connected elements of the current node @@ -138,7 +139,7 @@ RebarBondSlipConstraint::computeTangent() for (auto & elem : elems) { Elem * elem_ptr = _mesh.elemPtr(elem); - // _assembly.reinit(elem_ptr, 0); + _assembly.reinit(elem_ptr, 0); _current_elem_volume += _assembly.elemVolume(); // calculate phi and dphi for this element FEType fe_type(Utility::string_to_enum("first"), @@ -149,14 +150,15 @@ RebarBondSlipConstraint::computeTangent() unsigned side = 0; fe->reinit(elem_ptr, side); for (unsigned i = 0; i < tangents->size(); i++) - _slave_tangent += (*tangents)[i]; + _secondary_tangent += (*tangents)[i]; } - _slave_tangent /= _slave_tangent.norm(); + _secondary_tangent /= _secondary_tangent.norm(); _current_elem_volume /= elems.size(); if (_debug) - std::cout << "tangent: " << _slave_tangent << std::endl; + if (_current_node->id() == 144) + std::cout << "tangent: " << _secondary_tangent << std::endl; } void @@ -170,14 +172,14 @@ RebarBondSlipConstraint::reinitConstraint() for (unsigned int i = 0; i < _mesh_dimension; ++i) relative_disp(i) = ((_vars[i]->dofValues())[0] - (_vars[i]->slnNeighbor())[0]); - Real slip = relative_disp * _slave_tangent; - RealVectorValue slip_axial = slip * _slave_tangent; + Real slip = relative_disp * _secondary_tangent; + RealVectorValue slip_axial = slip * _secondary_tangent; RealVectorValue slip_normal = relative_disp - slip_axial; Real slip_ratio = std::abs(slip) / _transitional_slip[0]; // Real bond_stress; if (_debug) - // if (_current_node->id() == 144) - std::cout << "Slip = " << slip << ".\n"; + if (_current_node->id() == 144) + std::cout << "Slip = " << slip << ".\n"; const Node * node = _current_node; auto it = _bondslip.find(node->id()); @@ -188,100 +190,115 @@ RebarBondSlipConstraint::reinitConstraint() bond_slip.slip_max = std::max(bond_slip.slip_max_old, slip); if (_debug) - // if (_current_node->id() == 144) - { - std::cout << "Slip_min = " << bond_slip.slip_min << ".\n"; - std::cout << "Slip_min_old = " << bond_slip.slip_min_old << ".\n"; - std::cout << "Slip_max = " << bond_slip.slip_max << ".\n"; - std::cout << "Slip_max_old = " << bond_slip.slip_max_old << ".\n"; - std::cout << "Bondstress_min = " << bond_slip.bondstress_min << ".\n"; - std::cout << "Bondstress_min_old = " << bond_slip.bondstress_min_old << ".\n"; - std::cout << "Bondstress_max = " << bond_slip.bondstress_max << ".\n"; - std::cout << "Bondstress_max_old = " << bond_slip.bondstress_max_old << ".\n"; - } + if (_current_node->id() == 144) + { + std::cout << "Slip_min = " << bond_slip.slip_min << ".\n"; + std::cout << "Slip_min_old = " << bond_slip.slip_min_old << ".\n"; + std::cout << "Slip_max = " << bond_slip.slip_max << ".\n"; + std::cout << "Slip_max_old = " << bond_slip.slip_max_old << ".\n"; + std::cout << "Bondstress_min = " << bond_slip.bondstress_min << ".\n"; + std::cout << "Bondstress_min_old = " << bond_slip.bondstress_min_old << ".\n"; + std::cout << "Bondstress_max = " << bond_slip.bondstress_max << ".\n"; + std::cout << "Bondstress_max_old = " << bond_slip.bondstress_max_old << ".\n"; + } Real slope = 5.0 * _max_bondstress / _transitional_slip[0]; Real plastic_slip_max = bond_slip.slip_max - bond_slip.bondstress_max / slope; Real plastic_slip_min = bond_slip.slip_min - bond_slip.bondstress_min / slope; + _bond_stress_deriv = 0.0; + if (slip >= bond_slip.slip_max || slip <= bond_slip.slip_min) { if (std::abs(slip) < _transitional_slip[0]) { if (_debug) - // if (_current_node->id() == 144) - std::cout << "Calculating bondstress for Case Ia" - << ".\n"; + if (_current_node->id() == 144) + std::cout << "Calculating bondstress for Case Ia" + << ".\n"; _bond_stress = _max_bondstress * MathUtils::sign(slip) * (5.0 * slip_ratio - 4.5 * slip_ratio * slip_ratio + 1.4 * slip_ratio * slip_ratio * slip_ratio); + _bond_stress_deriv = + _max_bondstress * MathUtils::sign(slip) * + (5.0 / _transitional_slip[0] - 4.5 * 2.0 * slip_ratio / _transitional_slip[0] + + 1.4 * 3.0 * slip_ratio * slip_ratio / _transitional_slip[0]); } else if (slip >= _transitional_slip[0] && slip < _ultimate_slip) { if (_debug) - // if (_current_node->id() == 144) - std::cout << "Calculating bondstress for Case Ib" - << ".\n"; + if (_current_node->id() == 144) + std::cout << "Calculating bondstress for Case Ib" + << ".\n"; _bond_stress = 1.9 * _max_bondstress; } else if (slip <= -_transitional_slip[0] && slip > -_ultimate_slip) { if (_debug) - // if (_current_node->id() == 144) - std::cout << "Calculating bondstress for Case Ic" - << ".\n"; + if (_current_node->id() == 144) + std::cout << "Calculating bondstress for Case Ic" + << ".\n"; _bond_stress = -1.9 * _max_bondstress; } else { if (_debug) - // if (_current_node->id() == 144) - std::cout << "Calculating bondstress for Case Id" - << ".\n"; + if (_current_node->id() == 144) + std::cout << "Calculating bondstress for Case Id" + << ".\n"; _bond_stress = _frictional_bondstress * MathUtils::sign(slip); } } else if (slip > plastic_slip_max && slip < bond_slip.slip_max) { if (_debug) - // if (_current_node->id() == 144) - std::cout << "Calculating bondstress for Case II" - << ".\n"; + if (_current_node->id() == 144) + std::cout << "Calculating bondstress for Case II" + << ".\n"; _bond_stress = (slip - plastic_slip_max) * bond_slip.bondstress_max / (bond_slip.slip_max - plastic_slip_max); + + _bond_stress_deriv = bond_slip.bondstress_max / (bond_slip.slip_max - plastic_slip_max); } else if (slip < plastic_slip_min && slip > bond_slip.slip_min) { if (_debug) - // if (_current_node->id() == 144) - std::cout << "Calculating bondstress for Case III" - << ".\n"; + if (_current_node->id() == 144) + std::cout << "Calculating bondstress for Case III" + << ".\n"; _bond_stress = (slip - plastic_slip_min) * bond_slip.bondstress_min / (bond_slip.slip_min - plastic_slip_min); + _bond_stress_deriv = bond_slip.bondstress_min / (bond_slip.slip_min - plastic_slip_min); } else _bond_stress = _frictional_bondstress; if (_debug) - // if (_current_node->id() == 144) - std::cout << "Bondstress = " << _bond_stress << "\n"; + if (_current_node->id() == 144) + { + std::cout << "Bondstress = " << _bond_stress << "\n"; + std::cout << "Bondstress Derivative = " << _bond_stress_deriv << "\n"; + } Real bond_force = 2.0 * libMesh::pi * _bar_radius * _current_elem_volume * _bond_stress; + Real bond_force_deriv = + 2.0 * libMesh::pi * _bar_radius * _current_elem_volume * _bond_stress_deriv; - RealVectorValue constraint_force_axial = bond_force * _slave_tangent; + RealVectorValue constraint_force_axial = bond_force * _secondary_tangent; RealVectorValue constraint_force_normal = _penalty * slip_normal; _constraint_residual = constraint_force_axial + constraint_force_normal; + _constraint_jacobian_axial = bond_force_deriv * _secondary_tangent; if (_debug) - // if (_current_node->id() == 144) - { - std::cout << "Constraint Residual Axial = " << constraint_force_axial << "\n"; - std::cout << "Constraint Residual Normal = " << constraint_force_normal << "\n"; - std::cout << "Constraint Residual = " << _constraint_residual << "\n"; - } + if (_current_node->id() == 144) + { + std::cout << "Constraint Residual Axial = " << constraint_force_axial << "\n"; + std::cout << "Constraint Residual Normal = " << constraint_force_normal << "\n"; + std::cout << "Constraint Residual = " << _constraint_residual << "\n"; + } bond_slip.bondstress_min = std::min(bond_slip.bondstress_min_old, _bond_stress); bond_slip.bondstress_max = std::max(bond_slip.bondstress_max_old, _bond_stress); @@ -302,11 +319,11 @@ RebarBondSlipConstraint::computeQpResidual(Moose::ConstraintType type) switch (type) { - case Moose::Slave: - return _test_slave[_i][_qp] * resid; + case Moose::Secondary: + return resid * _test_secondary[_i][_qp]; - case Moose::Master: - return _test_master[_i][_qp] * -resid; + case Moose::Primary: + return -resid * _test_primary[_i][_qp]; } return 0.0; @@ -315,23 +332,77 @@ RebarBondSlipConstraint::computeQpResidual(Moose::ConstraintType type) Real RebarBondSlipConstraint::computeQpJacobian(Moose::ConstraintJacobianType type) { + Real jac_axial = _constraint_jacobian_axial(_component); + switch (type) { - case Moose::SlaveSlave: - return _phi_slave[_j][_qp] * _penalty * _test_slave[_i][_qp] * - (1.0 - _slave_tangent(_component) * _slave_tangent(_component)); + case Moose::SecondarySecondary: + return _phi_secondary[_j][_qp] * jac_axial * _secondary_tangent(_component) * + _test_secondary[_i][_qp] + + _phi_secondary[_j][_qp] * _penalty * _test_secondary[_i][_qp] * + (1.0 - _secondary_tangent(_component) * _secondary_tangent(_component)); + + case Moose::SecondaryPrimary: + return -_phi_primary[_j][_qp] * jac_axial * _secondary_tangent(_component) * + _test_secondary[_i][_qp] - + _phi_primary[_j][_qp] * _penalty * _test_secondary[_i][_qp] * + (1.0 - _secondary_tangent(_component) * _secondary_tangent(_component)); + + case Moose::PrimarySecondary: + return -_test_primary[_i][_qp] * jac_axial * _secondary_tangent(_component) * + _phi_secondary[_j][_qp] - + _test_primary[_i][_qp] * _penalty * _phi_secondary[_j][_qp] * + (1.0 - _secondary_tangent(_component) * _secondary_tangent(_component)); + + case Moose::PrimaryPrimary: + return _test_primary[_i][_qp] * jac_axial * _secondary_tangent(_component) * + _phi_primary[_j][_qp] + + _test_primary[_i][_qp] * _penalty * _phi_primary[_j][_qp] * + (1.0 - _secondary_tangent(_component) * _secondary_tangent(_component)); + + default: + mooseError("Unsupported type"); + break; + } + return 0.0; +} - case Moose::SlaveMaster: - return -_phi_master[_j][_qp] * _penalty * _test_slave[_i][_qp] * - (1.0 - _slave_tangent(_component) * _slave_tangent(_component)); +Real +RebarBondSlipConstraint::computeQpOffDiagJacobian(Moose::ConstraintJacobianType type, + unsigned int jvar) +{ + Real jac_axial = _constraint_jacobian_axial(_component); - case Moose::MasterSlave: - return -_test_master[_i][_qp] * _penalty * _phi_slave[_j][_qp] * - (1.0 - _slave_tangent(_component) * _slave_tangent(_component)); + unsigned int coupled_component; + Real tangent_component_in_coupled_var_dir = 1.0; + if (getCoupledVarComponent(jvar, coupled_component)) + tangent_component_in_coupled_var_dir = _secondary_tangent(coupled_component); - case Moose::MasterMaster: - return _test_master[_i][_qp] * _penalty * _phi_master[_j][_qp] * - (1.0 - _slave_tangent(_component) * _slave_tangent(_component)); + switch (type) + { + case Moose::SecondarySecondary: + return _phi_secondary[_j][_qp] * jac_axial * tangent_component_in_coupled_var_dir * + _test_secondary[_i][_qp] - + _phi_secondary[_j][_qp] * _penalty * _test_secondary[_i][_qp] * + _secondary_tangent(_component) * tangent_component_in_coupled_var_dir; + + case Moose::SecondaryPrimary: + return -_phi_primary[_j][_qp] * jac_axial * tangent_component_in_coupled_var_dir * + _test_secondary[_i][_qp] + + _phi_primary[_j][_qp] * _penalty * _test_secondary[_i][_qp] * + _secondary_tangent(_component) * tangent_component_in_coupled_var_dir; + + case Moose::PrimarySecondary: + return -_test_primary[_i][_qp] * jac_axial * tangent_component_in_coupled_var_dir * + _phi_secondary[_j][_qp] + + _test_primary[_i][_qp] * _penalty * _phi_secondary[_j][_qp] * + _secondary_tangent(_component) * tangent_component_in_coupled_var_dir; + + case Moose::PrimaryPrimary: + return _test_primary[_i][_qp] * jac_axial * tangent_component_in_coupled_var_dir * + _phi_primary[_j][_qp] - + _test_primary[_i][_qp] * _penalty * _phi_primary[_j][_qp] * + _secondary_tangent(_component) * tangent_component_in_coupled_var_dir; default: mooseError("Unsupported type"); @@ -340,9 +411,18 @@ RebarBondSlipConstraint::computeQpJacobian(Moose::ConstraintJacobianType type) return 0.0; } -Real -RebarBondSlipConstraint::computeQpOffDiagJacobian(Moose::ConstraintJacobianType /*type*/, - unsigned int /*jvar*/) +bool +RebarBondSlipConstraint::getCoupledVarComponent(unsigned int var_num, unsigned int & component) { - return 0.0; + component = std::numeric_limits::max(); + bool coupled_var_is_disp_var = false; + for (unsigned int i = 0; i < LIBMESH_DIM; ++i) + if (var_num == _var_nums[i]) + { + coupled_var_is_disp_var = true; + component = i; + break; + } + + return coupled_var_is_disp_var; } diff --git a/test/tests/rebar_bondslip/RCBeam_constraint.i b/test/tests/rebar_bondslip/RCBeam_constraint.i index c968f9b2..4796f04b 100644 --- a/test/tests/rebar_bondslip/RCBeam_constraint.i +++ b/test/tests/rebar_bondslip/RCBeam_constraint.i @@ -74,29 +74,31 @@ [Constraints] [rebar_x] type = RebarBondSlipConstraint - slave = 2 - master = 1 + secondary = 2 + primary = 1 penalty = 1e6 variable = 'disp_x' - master_variable = 'disp_x' + primary_variable = 'disp_x' component = 0 max_bondstress = 100 transitional_slip_values = 0.0005 ultimate_slip = 0.1 rebar_radius = 7.98e-3 + debug = true [] [rebar_y] type = RebarBondSlipConstraint - slave = 2 - master = 1 + secondary = 2 + primary = 1 penalty = 1e6 variable = 'disp_y' - master_variable = 'disp_y' + primary_variable = 'disp_y' component = 1 max_bondstress = 100 transitional_slip_values = 0.0005 ultimate_slip = 0.1 rebar_radius = 7.98e-3 + debug = true [] [] @@ -255,6 +257,7 @@ end_time = 30 dtmin = 0.00001 + # num_steps = 5 dt = 0.1 []