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 "HorizonStabilizedFormIFiniteStrainMechanicsNOSPD.h"
11 : #include "PeridynamicsMesh.h"
12 :
13 : registerMooseObject("PeridynamicsApp", HorizonStabilizedFormIFiniteStrainMechanicsNOSPD);
14 :
15 : InputParameters
16 17 : HorizonStabilizedFormIFiniteStrainMechanicsNOSPD::validParams()
17 : {
18 17 : InputParameters params = MechanicsFiniteStrainBaseNOSPD::validParams();
19 17 : params.addClassDescription(
20 : "Class for calculating the residual and the Jacobian for Form I "
21 : "of the horizon-stabilized peridynamic correspondence model under finite strain assumptions");
22 :
23 34 : 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 17 : return params;
28 0 : }
29 :
30 12 : HorizonStabilizedFormIFiniteStrainMechanicsNOSPD::HorizonStabilizedFormIFiniteStrainMechanicsNOSPD(
31 12 : const InputParameters & parameters)
32 24 : : MechanicsFiniteStrainBaseNOSPD(parameters), _component(getParam<unsigned int>("component"))
33 : {
34 12 : }
35 :
36 : void
37 146268 : HorizonStabilizedFormIFiniteStrainMechanicsNOSPD::computeLocalResidual()
38 : {
39 : // For finite strain formulation, the _stress tensor gotten from material class is the
40 : // Cauchy stress (Sigma). the first Piola-Kirchhoff stress (P) is then obtained as
41 : // P = J * Sigma * inv(F)^T.
42 : // Nodal force states are based on the first Piola-Kirchhoff stress tensors (P).
43 : // i.e., T = (J * Sigma * inv(F)^T) * inv(Shape) * xi * multi.
44 : // Cauchy stress is calculated as Sigma_n+1 = Sigma_n + R * (C * dt * D) * R^T
45 :
46 438804 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
47 877608 : for (unsigned int i = 0; i < _nnodes; ++i)
48 877608 : _local_re(i) += (i == 0 ? -1 : 1) * _multi[nd] *
49 1170144 : ((_dgrad[nd].det() * _stress[nd] * _dgrad[nd].inverse().transpose()) *
50 585072 : _shape2[nd].inverse())
51 1170144 : .row(_component) *
52 585072 : _origin_vec * _bond_status;
53 146268 : }
54 :
55 : void
56 9180 : HorizonStabilizedFormIFiniteStrainMechanicsNOSPD::computeLocalJacobian()
57 : {
58 : // excludes dTi/dUj and dTj/dUi contribution which was considered as nonlocal contribution
59 9180 : std::vector<RankTwoTensor> dPxdUx(_nnodes);
60 27540 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
61 18360 : dPxdUx[nd] = computeDJDU(_component, nd) * _stress[nd] * _dgrad[nd].inverse().transpose() +
62 36720 : _dgrad[nd].det() * computeDSDU(_component, nd) * _dgrad[nd].inverse().transpose() +
63 36720 : _dgrad[nd].det() * _stress[nd] * computeDinvFTDU(_component, nd);
64 :
65 27540 : for (unsigned int i = 0; i < _nnodes; ++i)
66 55080 : for (unsigned int j = 0; j < _nnodes; ++j)
67 55080 : _local_ke(i, j) += (i == 0 ? -1 : 1) * _multi[j] *
68 36720 : (dPxdUx[j] * _shape2[j].inverse()).row(_component) * _origin_vec *
69 36720 : _bond_status;
70 9180 : }
71 :
72 : void
73 0 : HorizonStabilizedFormIFiniteStrainMechanicsNOSPD::computeNonlocalJacobian()
74 : {
75 : // includes dTi/dUj and dTj/dUi contributions
76 0 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
77 : {
78 0 : RankFourTensor dSdFhat = computeDSDFhat(nd);
79 0 : RankTwoTensor invF = _dgrad[nd].inverse();
80 0 : Real detF = _dgrad[nd].det();
81 : // calculation of jacobian contribution to current_node's neighbors
82 0 : std::vector<dof_id_type> ivardofs(_nnodes);
83 0 : ivardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), _var.number(), 0);
84 0 : std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
85 0 : std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
86 :
87 : dof_id_type nb_index =
88 0 : std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - nd)) -
89 0 : neighbors.begin();
90 : std::vector<dof_id_type> dg_neighbors =
91 0 : _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(nd), nb_index);
92 :
93 : Real vol_nb, dJdU;
94 : RealGradient origin_vec_nb;
95 0 : RankTwoTensor dFdUk, dPxdUkx, dSdU, dinvFTdU;
96 :
97 0 : for (unsigned int nb = 0; nb < dg_neighbors.size(); ++nb)
98 0 : if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb]])) > 0.5)
99 : {
100 0 : ivardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])
101 0 : ->dof_number(_sys.number(), _var.number(), 0);
102 0 : vol_nb = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb]]);
103 :
104 0 : origin_vec_nb = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb]]) -
105 0 : _pdmesh.getNodeCoord(_current_elem->node_id(nd));
106 :
107 : dFdUk.zero();
108 0 : for (unsigned int i = 0; i < _dim; ++i)
109 0 : dFdUk(_component, i) =
110 0 : _horizon_radius[nd] / origin_vec_nb.norm() * origin_vec_nb(i) * vol_nb;
111 :
112 0 : dFdUk *= _shape2[nd].inverse();
113 :
114 : // calculate dJ/du
115 0 : dJdU = 0.0;
116 0 : for (unsigned int i = 0; i < 3; ++i)
117 0 : for (unsigned int j = 0; j < 3; ++j)
118 0 : dJdU += detF * invF(j, i) * dFdUk(i, j);
119 :
120 : // calculate dS/du
121 0 : dSdU = dSdFhat * dFdUk * _dgrad_old[nd].inverse();
122 :
123 : // calculate dinv(F)Tdu
124 : dinvFTdU.zero();
125 0 : for (unsigned int i = 0; i < 3; ++i)
126 0 : for (unsigned int J = 0; J < 3; ++J)
127 0 : for (unsigned int k = 0; k < 3; ++k)
128 0 : for (unsigned int L = 0; L < 3; ++L)
129 0 : dinvFTdU(i, J) += -invF(J, k) * invF(L, i) * dFdUk(k, L);
130 :
131 : // calculate the derivative of first Piola-Kirchhoff stress w.r.t displacements
132 0 : dPxdUkx = dJdU * _stress[nd] * invF.transpose() + detF * dSdU * invF.transpose() +
133 0 : detF * _stress[nd] * dinvFTdU;
134 :
135 0 : for (unsigned int i = 0; i < _nnodes; ++i)
136 0 : for (unsigned int j = 0; j < _nnodes; ++j)
137 0 : _local_ke(i, j) = (i == 0 ? -1 : 1) * (j == 0 ? 0 : 1) * _multi[nd] *
138 0 : (dPxdUkx * _shape2[nd].inverse()).row(_component) * _origin_vec *
139 0 : _bond_status;
140 :
141 0 : addJacobian(_assembly, _local_ke, _ivardofs, ivardofs, _var.scalingFactor());
142 :
143 0 : if (_has_diag_save_in)
144 : {
145 0 : unsigned int rows = _nnodes;
146 0 : DenseVector<Real> diag(rows);
147 0 : for (unsigned int i = 0; i < rows; ++i)
148 0 : diag(i) = _local_ke(i, i);
149 :
150 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
151 0 : for (unsigned int i = 0; i < _diag_save_in.size(); ++i)
152 : {
153 0 : std::vector<dof_id_type> diag_save_in_dofs(2);
154 0 : diag_save_in_dofs[0] = _current_elem->node_ptr(nd)->dof_number(
155 : _diag_save_in[i]->sys().number(), _diag_save_in[i]->number(), 0);
156 0 : diag_save_in_dofs[1] =
157 0 : _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])
158 0 : ->dof_number(_diag_save_in[i]->sys().number(), _diag_save_in[i]->number(), 0);
159 :
160 0 : _diag_save_in[i]->sys().solution().add_vector(diag, diag_save_in_dofs);
161 0 : }
162 : }
163 : }
164 0 : }
165 0 : }
166 :
167 : void
168 9180 : HorizonStabilizedFormIFiniteStrainMechanicsNOSPD::computeLocalOffDiagJacobian(
169 : unsigned int jvar_num, unsigned int coupled_component)
170 : {
171 : _local_ke.zero();
172 9180 : if (_temp_coupled && jvar_num == _temp_var->number())
173 : {
174 0 : std::vector<RankTwoTensor> dSdT(_nnodes);
175 0 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
176 0 : for (unsigned int es = 0; es < _deigenstrain_dT.size(); ++es)
177 0 : dSdT[nd] = -_dgrad[nd].det() * _Jacobian_mult[nd] * (*_deigenstrain_dT[es])[nd] *
178 0 : _dgrad[nd].inverse().transpose();
179 :
180 0 : for (unsigned int i = 0; i < _nnodes; ++i)
181 0 : for (unsigned int j = 0; j < _nnodes; ++j)
182 0 : _local_ke(i, j) += (i == 0 ? -1 : 1) * _multi[j] *
183 0 : (dSdT[j] * _shape2[j].inverse()).row(_component) * _origin_vec *
184 0 : _bond_status;
185 0 : }
186 9180 : else if (_out_of_plane_strain_coupled &&
187 0 : jvar_num == _out_of_plane_strain_var
188 : ->number()) // weak plane stress case, out_of_plane_strain is coupled
189 : {
190 0 : std::vector<RankTwoTensor> dSdE33(_nnodes);
191 0 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
192 : {
193 0 : for (unsigned int i = 0; i < 3; ++i)
194 0 : for (unsigned int j = 0; j < 3; ++j)
195 0 : dSdE33[nd](i, j) = _Jacobian_mult[nd](i, j, 2, 2);
196 :
197 0 : dSdE33[nd] = _dgrad[nd].det() * dSdE33[nd] * _dgrad[nd].inverse().transpose();
198 : }
199 :
200 0 : for (unsigned int i = 0; i < _nnodes; ++i)
201 0 : for (unsigned int j = 0; j < _nnodes; ++j)
202 0 : _local_ke(i, j) += (i == 0 ? -1 : 1) * _multi[j] *
203 0 : (dSdE33[j] * _shape2[j].inverse()).row(_component) * _origin_vec *
204 0 : _bond_status;
205 0 : }
206 : else
207 : {
208 9180 : std::vector<RankTwoTensor> dPxdUy(_nnodes);
209 27540 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
210 18360 : dPxdUy[nd] =
211 18360 : computeDJDU(coupled_component, nd) * _stress[nd] * _dgrad[nd].inverse().transpose() +
212 36720 : _dgrad[nd].det() * computeDSDU(coupled_component, nd) * _dgrad[nd].inverse().transpose() +
213 36720 : _dgrad[nd].det() * _stress[nd] * computeDinvFTDU(coupled_component, nd);
214 :
215 27540 : for (unsigned int i = 0; i < _nnodes; ++i)
216 55080 : for (unsigned int j = 0; j < _nnodes; ++j)
217 55080 : _local_ke(i, j) += (i == 0 ? -1 : 1) * _multi[j] *
218 36720 : (dPxdUy[j] * _shape2[j].inverse()).row(_component) * _origin_vec *
219 36720 : _bond_status;
220 9180 : }
221 9180 : }
222 :
223 : void
224 0 : HorizonStabilizedFormIFiniteStrainMechanicsNOSPD::computePDNonlocalOffDiagJacobian(
225 : unsigned int jvar_num, unsigned int coupled_component)
226 : {
227 0 : if (_temp_coupled && jvar_num == _temp_var->number())
228 : {
229 : // no nonlocal contribution from temperature
230 : }
231 0 : else if (_out_of_plane_strain_coupled && jvar_num == _out_of_plane_strain_var->number())
232 : {
233 : // no nonlocal contribution from out of plane strain
234 : }
235 : else
236 : {
237 0 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
238 : {
239 0 : RankFourTensor dSdFhat = computeDSDFhat(nd);
240 0 : RankTwoTensor invF = _dgrad[nd].inverse();
241 0 : Real detF = _dgrad[nd].det();
242 : // calculation of jacobian contribution to current_node's neighbors
243 0 : std::vector<dof_id_type> jvardofs(_nnodes);
244 0 : jvardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), jvar_num, 0);
245 0 : std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
246 0 : std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
247 :
248 : dof_id_type nb_index =
249 0 : std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - nd)) -
250 0 : neighbors.begin();
251 : std::vector<dof_id_type> dg_neighbors =
252 0 : _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(nd), nb_index);
253 :
254 : Real vol_nb, dJdU;
255 : RealGradient origin_vec_nb;
256 0 : RankTwoTensor dFdUk, dPxdUky, dSdU, dinvFTdU;
257 :
258 0 : for (unsigned int nb = 0; nb < dg_neighbors.size(); ++nb)
259 0 : if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb]])) > 0.5)
260 : {
261 0 : jvardofs[1] =
262 0 : _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])->dof_number(_sys.number(), jvar_num, 0);
263 0 : vol_nb = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb]]);
264 :
265 0 : origin_vec_nb = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb]]) -
266 0 : _pdmesh.getNodeCoord(_current_elem->node_id(nd));
267 :
268 : dFdUk.zero();
269 0 : for (unsigned int i = 0; i < _dim; ++i)
270 0 : dFdUk(coupled_component, i) =
271 0 : _horizon_radius[nd] / origin_vec_nb.norm() * origin_vec_nb(i) * vol_nb;
272 :
273 0 : dFdUk *= _shape2[nd].inverse();
274 :
275 : // calculate dJ/du
276 0 : dJdU = 0.0;
277 0 : for (unsigned int i = 0; i < 3; ++i)
278 0 : for (unsigned int j = 0; j < 3; ++j)
279 0 : dJdU += detF * invF(j, i) * dFdUk(i, j);
280 :
281 : // calculate dS/du
282 0 : dSdU = dSdFhat * dFdUk * _dgrad_old[nd].inverse();
283 :
284 : // calculate dinv(F)Tdu
285 : dinvFTdU.zero();
286 0 : for (unsigned int i = 0; i < 3; ++i)
287 0 : for (unsigned int J = 0; J < 3; ++J)
288 0 : for (unsigned int k = 0; k < 3; ++k)
289 0 : for (unsigned int L = 0; L < 3; ++L)
290 0 : dinvFTdU(i, J) += -invF(J, k) * invF(L, i) * dFdUk(k, L);
291 :
292 : // calculate the derivative of first Piola-Kirchhoff stress w.r.t displacements
293 0 : dPxdUky = dJdU * _stress[nd] * invF.transpose() + detF * dSdU * invF.transpose() +
294 0 : detF * _stress[nd] * dinvFTdU;
295 :
296 0 : for (unsigned int i = 0; i < _nnodes; ++i)
297 0 : for (unsigned int j = 0; j < _nnodes; ++j)
298 0 : _local_ke(i, j) = (i == 0 ? -1 : 1) * (j == 0 ? 0 : 1) * _multi[nd] *
299 0 : (dPxdUky * _shape2[nd].inverse()).row(_component) * _origin_vec *
300 0 : _bond_status;
301 :
302 0 : addJacobian(_assembly, _local_ke, _ivardofs, jvardofs, _var.scalingFactor());
303 : }
304 0 : }
305 : }
306 0 : }
|
