Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "HorizonStabilizedFormISmallStrainMechanicsNOSPD.h"
11 : #include "PeridynamicsMesh.h"
12 :
13 : registerMooseObject("PeridynamicsApp", HorizonStabilizedFormISmallStrainMechanicsNOSPD);
14 :
15 : InputParameters
16 306 : HorizonStabilizedFormISmallStrainMechanicsNOSPD::validParams()
17 : {
18 306 : InputParameters params = MechanicsBaseNOSPD::validParams();
19 306 : params.addClassDescription(
20 : "Class for calculating the residual and the Jacobian for Form I of the horizon-stabilized"
21 : "peridynamic correspondence model under small strain assumptions");
22 :
23 612 : params.addRequiredParam<unsigned int>(
24 : "component",
25 : "An integer corresponding to the variable this kernel acts on (0 for x, 1 for y, 2 for z)");
26 :
27 306 : return params;
28 0 : }
29 :
30 220 : HorizonStabilizedFormISmallStrainMechanicsNOSPD::HorizonStabilizedFormISmallStrainMechanicsNOSPD(
31 220 : const InputParameters & parameters)
32 440 : : MechanicsBaseNOSPD(parameters), _component(getParam<unsigned int>("component"))
33 : {
34 220 : }
35 :
36 : void
37 2097227 : HorizonStabilizedFormISmallStrainMechanicsNOSPD::computeLocalResidual()
38 : {
39 : // For small strain assumptions, stress measures, i.e., Cauchy stress (Sigma), the first
40 : // Piola-Kirchhoff stress (P), and the second Piola-Kirchhoff stress (S) are approximately the
41 : // same. Thus, the nodal force state tensors are calculated using the Cauchy stresses,
42 : // i.e., T = Sigma * inv(Shape) * xi * multi.
43 : // Cauchy stress is calculated as Sigma = C * E in the ComputeSmallStrainNOSPD material class.
44 :
45 6291681 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
46 12583362 : for (unsigned int i = 0; i < _nnodes; ++i)
47 12583362 : _local_re(i) += (i == 0 ? -1 : 1) * _multi[nd] *
48 8388908 : (_stress[nd] * _shape2[nd].inverse()).row(_component) * _origin_vec *
49 8388908 : _bond_status;
50 2097227 : }
51 :
52 : void
53 87923 : HorizonStabilizedFormISmallStrainMechanicsNOSPD::computeLocalJacobian()
54 : {
55 : // excludes dTi/dUj and dTj/dUi contributions, which were considered as nonlocal contribution
56 263769 : for (unsigned int i = 0; i < _nnodes; ++i)
57 527538 : for (unsigned int j = 0; j < _nnodes; ++j)
58 527538 : _local_ke(i, j) += (i == 0 ? -1 : 1) * _multi[j] *
59 703384 : (computeDSDU(_component, j) * _shape2[j].inverse()).row(_component) *
60 351692 : _origin_vec * _bond_status;
61 87923 : }
62 :
63 : void
64 8655 : HorizonStabilizedFormISmallStrainMechanicsNOSPD::computeNonlocalJacobian()
65 : {
66 : // includes dTi/dUj and dTj/dUi contributions
67 : // excludes contributions to node i and j, which were considered as local contributions
68 25965 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
69 : {
70 : // calculation of jacobian contribution to current_node's neighbors
71 17310 : std::vector<dof_id_type> ivardofs(_nnodes);
72 17310 : ivardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), _var.number(), 0);
73 17310 : std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
74 17310 : std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
75 :
76 : dof_id_type nb_index =
77 17310 : std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - nd)) -
78 17310 : neighbors.begin();
79 : std::vector<dof_id_type> dg_neighbors =
80 17310 : _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(nd), nb_index);
81 :
82 : Real vol_nb;
83 : RealGradient origin_vec_nb;
84 17310 : RankTwoTensor dFdUk, dPxdUkx;
85 :
86 245910 : for (unsigned int nb = 0; nb < dg_neighbors.size(); ++nb)
87 228600 : if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb]])) > 0.5)
88 : {
89 221638 : ivardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])
90 221638 : ->dof_number(_sys.number(), _var.number(), 0);
91 221638 : vol_nb = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb]]);
92 :
93 221638 : origin_vec_nb = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb]]) -
94 443276 : _pdmesh.getNodeCoord(_current_elem->node_id(nd));
95 :
96 : dFdUk.zero();
97 744912 : for (unsigned int i = 0; i < _dim; ++i)
98 523274 : dFdUk(_component, i) =
99 523274 : _horizon_radius[nd] / origin_vec_nb.norm() * origin_vec_nb(i) * vol_nb;
100 :
101 221638 : dFdUk *= _shape2[nd].inverse();
102 :
103 443276 : dPxdUkx = _Jacobian_mult[nd] * 0.5 * (dFdUk.transpose() + dFdUk);
104 :
105 664914 : for (unsigned int i = 0; i < _nnodes; ++i)
106 1329828 : for (unsigned int j = 0; j < _nnodes; ++j)
107 1329828 : _local_ke(i, j) = (i == 0 ? -1 : 1) * (j == 0 ? 0 : 1) * _multi[nd] *
108 886552 : (dPxdUkx * _shape2[nd].inverse()).row(_component) * _origin_vec *
109 886552 : _bond_status;
110 :
111 221638 : addJacobian(_assembly, _local_ke, _ivardofs, ivardofs, _var.scalingFactor());
112 :
113 221638 : if (_has_diag_save_in)
114 : {
115 0 : unsigned int rows = _nnodes;
116 0 : DenseVector<Real> diag(rows);
117 0 : for (unsigned int i = 0; i < rows; ++i)
118 0 : diag(i) = _local_ke(i, i);
119 :
120 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
121 0 : for (unsigned int i = 0; i < _diag_save_in.size(); ++i)
122 : {
123 0 : std::vector<dof_id_type> diag_save_in_dofs(2);
124 0 : diag_save_in_dofs[0] = _current_elem->node_ptr(nd)->dof_number(
125 : _diag_save_in[i]->sys().number(), _diag_save_in[i]->number(), 0);
126 0 : diag_save_in_dofs[1] =
127 0 : _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])
128 0 : ->dof_number(_diag_save_in[i]->sys().number(), _diag_save_in[i]->number(), 0);
129 :
130 0 : _diag_save_in[i]->sys().solution().add_vector(diag, diag_save_in_dofs);
131 0 : }
132 : }
133 : }
134 17310 : }
135 8655 : }
136 :
137 : void
138 107922 : HorizonStabilizedFormISmallStrainMechanicsNOSPD::computeLocalOffDiagJacobian(
139 : unsigned int jvar_num, unsigned int coupled_component)
140 : {
141 : _local_ke.zero();
142 107922 : if (_temp_coupled && jvar_num == _temp_var->number()) // temperature is coupled
143 : {
144 6588 : std::vector<RankTwoTensor> dSdT(_nnodes);
145 19764 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
146 26352 : for (unsigned int es = 0; es < _deigenstrain_dT.size(); ++es)
147 13176 : dSdT[nd] = -_Jacobian_mult[nd] * (*_deigenstrain_dT[es])[nd];
148 :
149 19764 : for (unsigned int i = 0; i < _nnodes; ++i)
150 39528 : for (unsigned int j = 0; j < _nnodes; ++j)
151 39528 : _local_ke(i, j) += (i == 0 ? -1 : 1) * _multi[j] *
152 26352 : (dSdT[j] * _shape2[j].inverse()).row(_component) * _origin_vec *
153 26352 : _bond_status;
154 6588 : }
155 101334 : else if (_out_of_plane_strain_coupled &&
156 27392 : jvar_num == _out_of_plane_strain_var
157 : ->number()) // weak plane stress case, out_of_plane_strain is coupled
158 : {
159 13696 : std::vector<RankTwoTensor> dSdE33(_nnodes);
160 41088 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
161 109568 : for (unsigned int i = 0; i < 3; ++i)
162 328704 : for (unsigned int j = 0; j < 3; ++j)
163 246528 : dSdE33[nd](i, j) = _Jacobian_mult[nd](i, j, 2, 2);
164 :
165 41088 : for (unsigned int i = 0; i < _nnodes; ++i)
166 82176 : for (unsigned int j = 0; j < _nnodes; ++j)
167 82176 : _local_ke(i, j) += (i == 0 ? -1 : 1) * _multi[j] *
168 54784 : (dSdE33[j] * _shape2[j].inverse()).row(_component) * _origin_vec *
169 54784 : _bond_status;
170 13696 : }
171 : else // off-diagonal Jacobian with respect to other displacement variables
172 : {
173 : // ONLY consider the contributions to node i and j
174 : // contributions to their neighbors are considered as nonlocal off-diagonal
175 262914 : for (unsigned int i = 0; i < _nnodes; ++i)
176 525828 : for (unsigned int j = 0; j < _nnodes; ++j)
177 350552 : _local_ke(i, j) +=
178 525828 : (i == 0 ? -1 : 1) * _multi[j] *
179 701104 : (computeDSDU(coupled_component, j) * _shape2[j].inverse()).row(_component) *
180 350552 : _origin_vec * _bond_status;
181 : }
182 107922 : }
183 :
184 : void
185 16134 : HorizonStabilizedFormISmallStrainMechanicsNOSPD::computePDNonlocalOffDiagJacobian(
186 : unsigned int jvar_num, unsigned int coupled_component)
187 : {
188 16134 : if (_temp_coupled && jvar_num == _temp_var->number())
189 : {
190 : // no nonlocal contribution from temperature
191 : }
192 12606 : else if (_out_of_plane_strain_coupled && jvar_num == _out_of_plane_strain_var->number())
193 : {
194 : // no nonlocal contribution from out of plane strain
195 : }
196 : else
197 : {
198 34290 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
199 : {
200 : // calculation of jacobian contribution to current_node's neighbors
201 : // NOT including the contribution to nodes i and j, which is considered as local off-diagonal
202 22860 : std::vector<dof_id_type> jvardofs(_nnodes);
203 22860 : jvardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), jvar_num, 0);
204 22860 : std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
205 22860 : std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
206 :
207 : dof_id_type nb_index =
208 22860 : std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - nd)) -
209 22860 : neighbors.begin();
210 : std::vector<dof_id_type> dg_neighbors =
211 22860 : _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(nd), nb_index);
212 :
213 : Real vol_nb;
214 : RealGradient origin_vec_nb;
215 22860 : RankTwoTensor dFdUk, dPxdUky;
216 :
217 332940 : for (unsigned int nb = 0; nb < dg_neighbors.size(); ++nb)
218 310080 : if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb]])) > 0.5)
219 : {
220 301636 : jvardofs[1] =
221 301636 : _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])->dof_number(_sys.number(), jvar_num, 0);
222 301636 : vol_nb = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb]]);
223 :
224 301636 : origin_vec_nb = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb]]) -
225 603272 : _pdmesh.getNodeCoord(_current_elem->node_id(nd));
226 :
227 : dFdUk.zero();
228 1064904 : for (unsigned int i = 0; i < _dim; ++i)
229 763268 : dFdUk(coupled_component, i) =
230 763268 : _horizon_radius[nd] / origin_vec_nb.norm() * origin_vec_nb(i) * vol_nb;
231 :
232 301636 : dFdUk *= _shape2[nd].inverse();
233 :
234 603272 : dPxdUky = _Jacobian_mult[nd] * 0.5 * (dFdUk.transpose() + dFdUk);
235 :
236 904908 : for (unsigned int i = 0; i < _nnodes; ++i)
237 1809816 : for (unsigned int j = 0; j < _nnodes; ++j)
238 1809816 : _local_ke(i, j) = (i == 0 ? -1 : 1) * (j == 0 ? 0 : 1) * _multi[nd] *
239 1206544 : (dPxdUky * _shape2[nd].inverse()).row(_component) * _origin_vec *
240 1206544 : _bond_status;
241 :
242 301636 : addJacobian(_assembly, _local_ke, _ivardofs, jvardofs, _var.scalingFactor());
243 : }
244 22860 : }
245 : }
246 16134 : }
|
