https://mooseframework.inl.gov
Loading...
Searching...
No Matches
HorizonStabilizedFormIIFiniteStrainMechanicsNOSPD.C
Go to the documentation of this file.
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 "HorizonStabilizedFormIIFiniteStrainMechanicsNOSPD.h"
11#include "PeridynamicsMesh.h"
12
13registerMooseObject("PeridynamicsApp", HorizonStabilizedFormIIFiniteStrainMechanicsNOSPD);
14
15InputParameters
16HorizonStabilizedFormIIFiniteStrainMechanicsNOSPD::validParams()
17{
18 InputParameters params = MechanicsFiniteStrainBaseNOSPD::validParams();
19 params.addClassDescription(
20 "Class for calculating the residual and the Jacobian for Form II "
21 "of the horizon-stabilized peridynamic correspondence model under finite strain assumptions");
22
23 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 return params;
28}
29
35
36void
37HorizonStabilizedFormIIFiniteStrainMechanicsNOSPD::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 for (unsigned int nd = 0; nd < _nnodes; ++nd)
47 for (unsigned int i = 0; i < _nnodes; ++i)
48 _local_re(i) += (i == 0 ? -1 : 1) * _multi[nd] * _horizon_radius[nd] / _origin_vec.norm() *
49 ((_dgrad[nd].det() * _stress[nd] * _dgrad[nd].inverse().transpose()) *
50 _shape2[nd].inverse())
51 .row(_component) *
52 _origin_vec * _node_vol[1 - nd] * _bond_status;
53}
54
55void
56HorizonStabilizedFormIIFiniteStrainMechanicsNOSPD::computeNonlocalResidual()
57{
58 for (unsigned int nd = 0; nd < _nnodes; ++nd)
59 {
60 // calculation of residual contribution to current_node's neighbors
61 std::vector<dof_id_type> ivardofs(_nnodes);
62 ivardofs[0] = _ivardofs[nd];
63 std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
64 std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
65
66 dof_id_type nb_index =
67 std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - nd)) -
68 neighbors.begin();
69 std::vector<dof_id_type> dg_neighbors =
70 _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(nd), nb_index);
71
72 RealGradient origin_vec_nb;
73 Real node_vol_nb;
74
75 for (unsigned int nb = 0; nb < dg_neighbors.size(); ++nb)
76 if (neighbors[dg_neighbors[nb]] != _current_elem->node_id(1 - nd) &&
77 _bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb]])) > 0.5)
78 {
79 ivardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])
80 ->dof_number(_sys.number(), _var.number(), 0);
81 origin_vec_nb = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb]]) -
82 _pdmesh.getNodeCoord(_current_elem->node_id(nd));
83 node_vol_nb = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb]]);
84
85 for (unsigned int i = 0; i < _nnodes; ++i)
86 _local_re(i) = (i == 0 ? -1 : 1) * _multi[nd] * _horizon_radius[nd] /
87 origin_vec_nb.norm() *
88 ((_dgrad[nd].det() * _stress[nd] * _dgrad[nd].inverse().transpose()) *
89 _shape2[nd].inverse())
90 .row(_component) *
91 origin_vec_nb * node_vol_nb * _bond_status;
92
93 // cache the residual contribution
94 addResiduals(_assembly, _local_re, ivardofs, _var.scalingFactor());
95
96 if (_has_save_in)
97 {
98 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
99 for (unsigned int i = 0; i < _save_in.size(); ++i)
100 {
101 std::vector<dof_id_type> save_in_dofs(_nnodes);
102 save_in_dofs[0] = _current_elem->node_ptr(nd)->dof_number(
103 _save_in[i]->sys().number(), _save_in[i]->number(), 0);
104 save_in_dofs[1] =
105 _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])
106 ->dof_number(_save_in[i]->sys().number(), _save_in[i]->number(), 0);
107
108 _save_in[i]->sys().solution().add_vector(_local_re, save_in_dofs);
109 }
110 }
111 }
112 }
113}
114
115void
116HorizonStabilizedFormIIFiniteStrainMechanicsNOSPD::computeLocalJacobian()
117{
118 for (unsigned int nd = 0; nd < _nnodes; ++nd)
119 {
120 std::vector<dof_id_type> ivardofs(_nnodes);
121 ivardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), _var.number(), 0);
122 std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
123 std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
124
125 dof_id_type nb_index =
126 std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - nd)) -
127 neighbors.begin();
128 std::vector<dof_id_type> dg_neighbors =
129 _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(nd), nb_index);
130
131 RankTwoTensor dPxdUx =
132 computeDJDU(_component, nd) * _stress[nd] * _dgrad[nd].inverse().transpose() +
133 _dgrad[nd].det() * computeDSDU(_component, nd) * _dgrad[nd].inverse().transpose() +
134 _dgrad[nd].det() * _stress[nd] * computeDinvFTDU(_component, nd);
135
136 RealGradient origin_vec_nb;
137 Real node_vol_nb;
138
139 for (unsigned int nb = 0; nb < dg_neighbors.size(); ++nb)
140 if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb]])) > 0.5)
141 {
142 ivardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])
143 ->dof_number(_sys.number(), _var.number(), 0);
144 origin_vec_nb = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb]]) -
145 _pdmesh.getNodeCoord(_current_elem->node_id(nd));
146 node_vol_nb = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb]]);
147
148 for (unsigned int i = 0; i < _nnodes; ++i)
149 for (unsigned int j = 0; j < _nnodes; ++j)
150 _local_ke(i, j) = (i == 0 ? -1 : 1) * (j == 0 ? 1 : 0) * _multi[nd] *
151 _horizon_radius[nd] / origin_vec_nb.norm() *
152 (dPxdUx * _shape2[nd].inverse()).row(_component) * origin_vec_nb *
153 node_vol_nb * _bond_status;
154
155 addJacobian(_assembly, _local_ke, ivardofs, ivardofs, _var.scalingFactor());
156 }
157 }
158 _local_ke.zero();
159}
160
161void
162HorizonStabilizedFormIIFiniteStrainMechanicsNOSPD::computeNonlocalJacobian()
163{
164 // includes dTi/dUj and dTj/dUi contributions
165 for (unsigned int nd = 0; nd < _nnodes; ++nd)
166 {
167 RankFourTensor dSdFhat = computeDSDFhat(nd);
168 RankTwoTensor invF = _dgrad[nd].inverse();
169 Real detF = _dgrad[nd].det();
170 // calculation of jacobian contribution to current_node's neighbors
171 std::vector<dof_id_type> ivardofs(_nnodes), jvardofs(_nnodes);
172 ivardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), _var.number(), 0);
173 jvardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), _var.number(), 0);
174 std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
175 std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
176
177 dof_id_type nb_index =
178 std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - nd)) -
179 neighbors.begin();
180 std::vector<dof_id_type> dg_neighbors =
181 _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(nd), nb_index);
182
183 RealGradient origin_vec_nb1;
184 Real node_vol_nb1;
185
186 for (unsigned int nb1 = 0; nb1 < dg_neighbors.size(); ++nb1)
187 if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb1]])) > 0.5)
188 {
189 ivardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb1]])
190 ->dof_number(_sys.number(), _var.number(), 0);
191 origin_vec_nb1 = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb1]]) -
192 _pdmesh.getNodeCoord(_current_elem->node_id(nd));
193 node_vol_nb1 = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb1]]);
194
195 Real vol_nb2, dJdU;
196 RealGradient origin_vec_nb2;
197 RankTwoTensor dFdUk, dPxdUkx, dSdU, dinvFTdU;
198
199 for (unsigned int nb2 = 0; nb2 < dg_neighbors.size(); ++nb2)
200 if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb2]])) > 0.5)
201 {
202 ivardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb2]])
203 ->dof_number(_sys.number(), _var.number(), 0);
204 vol_nb2 = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb2]]);
205
206 origin_vec_nb2 = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb2]]) -
207 _pdmesh.getNodeCoord(_current_elem->node_id(nd));
208
209 dFdUk.zero();
210 for (unsigned int i = 0; i < _dim; ++i)
211 dFdUk(_component, i) =
212 _horizon_radius[nd] / origin_vec_nb2.norm() * origin_vec_nb2(i) * vol_nb2;
213
214 dFdUk *= _shape2[nd].inverse();
215
216 // calculate dJ/du
217 dJdU = 0.0;
218 for (unsigned int i = 0; i < 3; ++i)
219 for (unsigned int j = 0; j < 3; ++j)
220 dJdU += detF * invF(j, i) * dFdUk(i, j);
221
222 // calculate dS/du
223 dSdU = dSdFhat * dFdUk * _dgrad_old[nd].inverse();
224
225 // calculate dinv(F)Tdu
226 dinvFTdU.zero();
227 for (unsigned int i = 0; i < 3; ++i)
228 for (unsigned int J = 0; J < 3; ++J)
229 for (unsigned int k = 0; k < 3; ++k)
230 for (unsigned int L = 0; L < 3; ++L)
231 dinvFTdU(i, J) += -invF(J, k) * invF(L, i) * dFdUk(k, L);
232
233 // calculate the derivative of first Piola-Kirchhoff stress w.r.t displacements
234 dPxdUkx = dJdU * _stress[nd] * invF.transpose() + detF * dSdU * invF.transpose() +
235 detF * _stress[nd] * dinvFTdU;
236
237 for (unsigned int i = 0; i < _nnodes; ++i)
238 for (unsigned int j = 0; j < _nnodes; ++j)
239 _local_ke(i, j) = (i == 0 ? -1 : 1) * (j == 0 ? 0 : 1) * _multi[nd] *
240 _horizon_radius[nd] / origin_vec_nb1.norm() *
241 (dPxdUkx * _shape2[nd].inverse()).row(_component) *
242 origin_vec_nb1 * node_vol_nb1 * _bond_status;
243
244 addJacobian(_assembly, _local_ke, ivardofs, jvardofs, _var.scalingFactor());
245
246 if (_has_diag_save_in)
247 {
248 unsigned int rows = _nnodes;
249 DenseVector<Real> diag(rows);
250 for (unsigned int i = 0; i < rows; ++i)
251 diag(i) = _local_ke(i, i);
252
253 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
254 for (unsigned int i = 0; i < _diag_save_in.size(); ++i)
255 {
256 std::vector<dof_id_type> diag_save_in_dofs(2);
257 diag_save_in_dofs[0] = _current_elem->node_ptr(nd)->dof_number(
258 _diag_save_in[i]->sys().number(), _diag_save_in[i]->number(), 0);
259 diag_save_in_dofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb2]])
260 ->dof_number(_diag_save_in[i]->sys().number(),
261 _diag_save_in[i]->number(),
262 0);
263
264 _diag_save_in[i]->sys().solution().add_vector(diag, diag_save_in_dofs);
265 }
266 }
267 }
268 }
269 }
270}
271
272void
273HorizonStabilizedFormIIFiniteStrainMechanicsNOSPD::computeLocalOffDiagJacobian(
274 unsigned int jvar_num, unsigned int coupled_component)
275{
276 if (_temp_coupled && jvar_num == _temp_var->number()) // temperature is coupled
277 {
278 std::vector<RankTwoTensor> dSdT(_nnodes);
279 for (unsigned int nd = 0; nd < _nnodes; ++nd)
280 for (unsigned int es = 0; es < _deigenstrain_dT.size(); ++es)
281 dSdT[nd] = -_dgrad[nd].det() * _Jacobian_mult[nd] * (*_deigenstrain_dT[es])[nd] *
282 _dgrad[nd].inverse().transpose();
283
284 for (unsigned int nd = 0; nd < _nnodes; ++nd)
285 {
286 std::vector<dof_id_type> ivardofs(_nnodes), jvardofs(_nnodes);
287 ivardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), _var.number(), 0);
288 jvardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), jvar_num, 0);
289 std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
290 std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
291
292 dof_id_type nb_index =
293 std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - nd)) -
294 neighbors.begin();
295 std::vector<dof_id_type> dg_neighbors =
296 _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(nd), nb_index);
297
298 RealGradient origin_vec_nb;
299 Real node_vol_nb;
300
301 for (unsigned int nb = 0; nb < dg_neighbors.size(); ++nb)
302 if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb]])) > 0.5)
303 {
304 ivardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])
305 ->dof_number(_sys.number(), _var.number(), 0);
306 jvardofs[1] =
307 _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])->dof_number(_sys.number(), jvar_num, 0);
308 origin_vec_nb = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb]]) -
309 _pdmesh.getNodeCoord(_current_elem->node_id(nd));
310 node_vol_nb = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb]]);
311
312 for (unsigned int i = 0; i < _nnodes; ++i)
313 for (unsigned int j = 0; j < _nnodes; ++j)
314 _local_ke(i, j) = (i == 0 ? -1 : 1) * (j == 0 ? 1 : 0) * _multi[nd] *
315 _horizon_radius[nd] / origin_vec_nb.norm() *
316 (dSdT[nd] * _shape2[nd].inverse()).row(_component) * origin_vec_nb *
317 node_vol_nb * _bond_status;
318
319 addJacobian(_assembly, _local_ke, ivardofs, jvardofs, _var.scalingFactor());
320 }
321 }
322 _local_ke.zero();
323 }
324 else if (_out_of_plane_strain_coupled &&
325 jvar_num == _out_of_plane_strain_var
326 ->number()) // weak plane stress case, out_of_plane_strain is coupled
327 {
328 std::vector<RankTwoTensor> dSdE33(_nnodes);
329 for (unsigned int nd = 0; nd < _nnodes; ++nd)
330 {
331 for (unsigned int i = 0; i < 3; ++i)
332 for (unsigned int j = 0; j < 3; ++j)
333 dSdE33[nd](i, j) = _Jacobian_mult[nd](i, j, 2, 2);
334
335 dSdE33[nd] = _dgrad[nd].det() * dSdE33[nd] * _dgrad[nd].inverse().transpose();
336 }
337
338 for (unsigned int nd = 0; nd < _nnodes; ++nd)
339 {
340 std::vector<dof_id_type> ivardofs(_nnodes), jvardofs(_nnodes);
341 ivardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), _var.number(), 0);
342 jvardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), jvar_num, 0);
343 std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
344 std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
345
346 dof_id_type nb_index =
347 std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - nd)) -
348 neighbors.begin();
349 std::vector<dof_id_type> dg_neighbors =
350 _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(nd), nb_index);
351
352 RealGradient origin_vec_nb;
353 Real node_vol_nb;
354
355 for (unsigned int nb = 0; nb < dg_neighbors.size(); ++nb)
356 if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb]])) > 0.5)
357 {
358 ivardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])
359 ->dof_number(_sys.number(), _var.number(), 0);
360 jvardofs[1] =
361 _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])->dof_number(_sys.number(), jvar_num, 0);
362 origin_vec_nb = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb]]) -
363 _pdmesh.getNodeCoord(_current_elem->node_id(nd));
364 node_vol_nb = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb]]);
365
366 for (unsigned int i = 0; i < _nnodes; ++i)
367 for (unsigned int j = 0; j < _nnodes; ++j)
368 _local_ke(i, j) = (i == 0 ? -1 : 1) * (j == 0 ? 1 : 0) * _multi[nd] *
369 _horizon_radius[nd] / origin_vec_nb.norm() *
370 (dSdE33[nd] * _shape2[nd].inverse()).row(_component) *
371 origin_vec_nb * node_vol_nb * _bond_status;
372
373 addJacobian(_assembly, _local_ke, ivardofs, jvardofs, _var.scalingFactor());
374 }
375 }
376 _local_ke.zero();
377 }
378 else // off-diagonal Jacobian with respect to other displacement variables
379 {
380 for (unsigned int nd = 0; nd < _nnodes; ++nd)
381 {
382 std::vector<dof_id_type> ivardofs(_nnodes), jvardofs(_nnodes);
383 ivardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), _var.number(), 0);
384 jvardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), jvar_num, 0);
385 std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
386 std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
387
388 dof_id_type nb_index =
389 std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - nd)) -
390 neighbors.begin();
391 std::vector<dof_id_type> dg_neighbors =
392 _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(nd), nb_index);
393
394 RankTwoTensor dPxdUy =
395 computeDJDU(coupled_component, nd) * _stress[nd] * _dgrad[nd].inverse().transpose() +
396 _dgrad[nd].det() * computeDSDU(coupled_component, nd) * _dgrad[nd].inverse().transpose() +
397 _dgrad[nd].det() * _stress[nd] * computeDinvFTDU(coupled_component, nd);
398
399 RealGradient origin_vec_nb;
400 Real node_vol_nb;
401
402 for (unsigned int nb = 0; nb < dg_neighbors.size(); ++nb)
403 if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb]])) > 0.5)
404 {
405 ivardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])
406 ->dof_number(_sys.number(), _var.number(), 0);
407 jvardofs[1] =
408 _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])->dof_number(_sys.number(), jvar_num, 0);
409 origin_vec_nb = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb]]) -
410 _pdmesh.getNodeCoord(_current_elem->node_id(nd));
411 node_vol_nb = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb]]);
412
413 for (unsigned int i = 0; i < _nnodes; ++i)
414 for (unsigned int j = 0; j < _nnodes; ++j)
415 _local_ke(i, j) = (i == 0 ? -1 : 1) * (j == 0 ? 1 : 0) * _multi[nd] *
416 _horizon_radius[nd] / origin_vec_nb.norm() *
417 (dPxdUy * _shape2[nd].inverse()).row(_component) * origin_vec_nb *
418 node_vol_nb * _bond_status;
419
420 addJacobian(_assembly, _local_ke, ivardofs, jvardofs, _var.scalingFactor());
421 }
422 }
423 _local_ke.zero();
424 }
425}
426
427void
428HorizonStabilizedFormIIFiniteStrainMechanicsNOSPD::computePDNonlocalOffDiagJacobian(
429 unsigned int jvar_num, unsigned int coupled_component)
430{
431 if (_temp_coupled && jvar_num == _temp_var->number())
432 {
433 // no nonlocal contribution from temperature
434 }
435 else if (_out_of_plane_strain_coupled && jvar_num == _out_of_plane_strain_var->number())
436 {
437 // no nonlocal contribution from out of plane strain
438 }
439 else
440 {
441 for (unsigned int nd = 0; nd < _nnodes; ++nd)
442 {
443 RankFourTensor dSdFhat = computeDSDFhat(nd);
444 RankTwoTensor invF = _dgrad[nd].inverse();
445 Real detF = _dgrad[nd].det();
446 // calculation of jacobian contribution to current_node's neighbors
447 // NOT including the contribution to nodes i and j, which is considered as local off-diagonal
448 std::vector<dof_id_type> ivardofs(_nnodes), jvardofs(_nnodes);
449 ivardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), _var.number(), 0);
450 jvardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), jvar_num, 0);
451 std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
452 std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
453
454 dof_id_type nb_index =
455 std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - nd)) -
456 neighbors.begin();
457 std::vector<dof_id_type> dg_neighbors =
458 _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(nd), nb_index);
459
460 RealGradient origin_vec_nb1;
461 Real node_vol_nb1;
462
463 for (unsigned int nb1 = 0; nb1 < dg_neighbors.size(); ++nb1)
464 if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb1]])) > 0.5)
465 {
466 ivardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb1]])
467 ->dof_number(_sys.number(), _var.number(), 0);
468 origin_vec_nb1 = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb1]]) -
469 _pdmesh.getNodeCoord(_current_elem->node_id(nd));
470 node_vol_nb1 = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb1]]);
471
472 Real vol_nb2, dJdU;
473 RealGradient origin_vec_nb2;
474 RankTwoTensor dFdUk, dPxdUky, dSdU, dinvFTdU;
475
476 for (unsigned int nb2 = 0; nb2 < dg_neighbors.size(); ++nb2)
477 if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb2]])) >
478 0.5)
479 {
480 jvardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb2]])
481 ->dof_number(_sys.number(), jvar_num, 0);
482 vol_nb2 = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb2]]);
483
484 origin_vec_nb2 = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb2]]) -
485 _pdmesh.getNodeCoord(_current_elem->node_id(nd));
486
487 dFdUk.zero();
488 for (unsigned int i = 0; i < _dim; ++i)
489 dFdUk(coupled_component, i) =
490 _horizon_radius[nd] / origin_vec_nb2.norm() * origin_vec_nb2(i) * vol_nb2;
491
492 dFdUk *= _shape2[nd].inverse();
493
494 // calculate dJ/du
495 dJdU = 0.0;
496 for (unsigned int i = 0; i < 3; ++i)
497 for (unsigned int j = 0; j < 3; ++j)
498 dJdU += detF * invF(j, i) * dFdUk(i, j);
499
500 // calculate dS/du
501 dSdU = dSdFhat * dFdUk * _dgrad_old[nd].inverse();
502
503 // calculate dinv(F)Tdu
504 dinvFTdU.zero();
505 for (unsigned int i = 0; i < 3; ++i)
506 for (unsigned int J = 0; J < 3; ++J)
507 for (unsigned int k = 0; k < 3; ++k)
508 for (unsigned int L = 0; L < 3; ++L)
509 dinvFTdU(i, J) += -invF(J, k) * invF(L, i) * dFdUk(k, L);
510
511 // calculate the derivative of first Piola-Kirchhoff stress w.r.t displacements
512 dPxdUky = dJdU * _stress[nd] * invF.transpose() + detF * dSdU * invF.transpose() +
513 detF * _stress[nd] * dinvFTdU;
514
515 for (unsigned int i = 0; i < _nnodes; ++i)
516 for (unsigned int j = 0; j < _nnodes; ++j)
517 _local_ke(i, j) = (i == 0 ? -1 : 1) * (j == 0 ? 0 : 1) * _multi[nd] *
518 _horizon_radius[nd] / origin_vec_nb1.norm() *
519 (dPxdUky * _shape2[nd].inverse()).row(_component) *
520 origin_vec_nb1 * node_vol_nb1 * _bond_status;
521
522 addJacobian(_assembly, _local_ke, ivardofs, jvardofs, _var.scalingFactor());
523 }
524 }
525 }
526 }
527}
registerMooseObject
registerMooseObject("PeridynamicsApp", HorizonStabilizedFormIIFiniteStrainMechanicsNOSPD)
int
void ErrorVector unsigned int
HorizonStabilizedFormIIFiniteStrainMechanicsNOSPD
Kernel class for Form II of the horizon-stabilized peridynamic correspondence model for finite strain...
Definition HorizonStabilizedFormIIFiniteStrainMechanicsNOSPD.h:21
HorizonStabilizedFormIIFiniteStrainMechanicsNOSPD::computeLocalOffDiagJacobian
virtual void computeLocalOffDiagJacobian(unsigned int jvar_num, unsigned int coupled_component) override
Function to compute local contribution to the off-diagonal Jacobian at the current nodes.
Definition HorizonStabilizedFormIIFiniteStrainMechanicsNOSPD.C:273
HorizonStabilizedFormIIFiniteStrainMechanicsNOSPD::computePDNonlocalOffDiagJacobian
virtual void computePDNonlocalOffDiagJacobian(unsigned int jvar_num, unsigned int coupled_component) override
Function to compute nonlocal contribution to the off-diagonal Jacobian at the current nodes.
Definition HorizonStabilizedFormIIFiniteStrainMechanicsNOSPD.C:428
InputParameters::addRequiredParam
void addRequiredParam(const std::string &name, const std::string &doc_string)
InputParameters::addClassDescription
void addClassDescription(const std::string &doc_string)
MechanicsBaseNOSPD::_shape2
const MaterialProperty< RankTwoTensor > & _shape2
Definition MechanicsBaseNOSPD.h:37
MechanicsBaseNOSPD::_stress
const MaterialProperty< RankTwoTensor > & _stress
Definition MechanicsBaseNOSPD.h:36
MechanicsBaseNOSPD::_deigenstrain_dT
std::vector< const MaterialProperty< RankTwoTensor > * > _deigenstrain_dT
Definition MechanicsBaseNOSPD.h:44
MechanicsBaseNOSPD::_multi
const MaterialProperty< Real > & _multi
Material point based material properties.
Definition MechanicsBaseNOSPD.h:35
MechanicsBaseNOSPD::_Jacobian_mult
const MaterialProperty< RankFourTensor > & _Jacobian_mult
Definition MechanicsBaseNOSPD.h:42
MechanicsBaseNOSPD::_dgrad
const MaterialProperty< RankTwoTensor > & _dgrad
Definition MechanicsBaseNOSPD.h:38
MechanicsBasePD::_out_of_plane_strain_coupled
const bool _out_of_plane_strain_coupled
Parameters for out-of-plane strain in weak plane stress formulation.
Definition MechanicsBasePD.h:62
MechanicsBasePD::_temp_var
MooseVariable * _temp_var
Definition MechanicsBasePD.h:55
MechanicsBasePD::_out_of_plane_strain_var
MooseVariable * _out_of_plane_strain_var
Definition MechanicsBasePD.h:63
MechanicsBasePD::_ivardofs
std::vector< dof_id_type > _ivardofs
Current variable dof numbers for nodes i and j.
Definition MechanicsBasePD.h:70
MechanicsBasePD::_temp_coupled
const bool _temp_coupled
Temperature variable.
Definition MechanicsBasePD.h:54
MechanicsFiniteStrainBaseNOSPD
Base kernel class for finite strain correspondence models.
Definition MechanicsFiniteStrainBaseNOSPD.h:18
MechanicsFiniteStrainBaseNOSPD::computeDJDU
Real computeDJDU(unsigned int component, unsigned int nd)
Function to compute derivative of determinant of deformation gradient with respect to displacements.
Definition MechanicsFiniteStrainBaseNOSPD.C:110
MechanicsFiniteStrainBaseNOSPD::computeDinvFTDU
RankTwoTensor computeDinvFTDU(unsigned int component, unsigned int nd)
Function to compute derivative of deformation gradient inverse with respect to displacements.
Definition MechanicsFiniteStrainBaseNOSPD.C:134
MechanicsFiniteStrainBaseNOSPD::computeDSDFhat
RankFourTensor computeDSDFhat(unsigned int nd)
Function to compute derivative of stress with respect to derived deformation gradient.
Definition MechanicsFiniteStrainBaseNOSPD.C:61
MechanicsFiniteStrainBaseNOSPD::computeDSDU
virtual RankTwoTensor computeDSDU(unsigned int component, unsigned int nd) override
Function to compute derivative of stress with respect to displacements for small strain problems.
Definition MechanicsFiniteStrainBaseNOSPD.C:36
MechanicsFiniteStrainBaseNOSPD::_dgrad_old
const MaterialProperty< RankTwoTensor > & _dgrad_old
Material point based material property.
Definition MechanicsFiniteStrainBaseNOSPD.h:52
MooseVariableFE::number
unsigned int number() const
RankTwoTensorTempl::inverse
RankTwoTensorTempl< T > inverse() const
RankTwoTensorTempl::transpose
RankTwoTensorTempl< T > transpose() const
Generated on Fri Sep 11 2026 20:36:24 for https://mooseframework.inl.gov by doxygen 1.9.8