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 "HorizonStabilizedFormIISmallStrainMechanicsNOSPD.h"
11 : #include "PeridynamicsMesh.h"
12 :
13 : registerMooseObject("PeridynamicsApp", HorizonStabilizedFormIISmallStrainMechanicsNOSPD);
14 :
15 : InputParameters
16 47 : HorizonStabilizedFormIISmallStrainMechanicsNOSPD::validParams()
17 : {
18 47 : InputParameters params = MechanicsBaseNOSPD::validParams();
19 47 : params.addClassDescription(
20 : "Class for calculating the residual and the Jacobian for Form II "
21 : "of the horizon-stabilized peridynamic correspondence model under small strain assumptions");
22 :
23 94 : 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 47 : return params;
28 0 : }
29 :
30 35 : HorizonStabilizedFormIISmallStrainMechanicsNOSPD::HorizonStabilizedFormIISmallStrainMechanicsNOSPD(
31 35 : const InputParameters & parameters)
32 70 : : MechanicsBaseNOSPD(parameters), _component(getParam<unsigned int>("component"))
33 : {
34 35 : }
35 :
36 : void
37 193269 : HorizonStabilizedFormIISmallStrainMechanicsNOSPD::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 579807 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
46 1159614 : for (unsigned int i = 0; i < _nnodes; ++i)
47 1159614 : _local_re(i) += (i == 0 ? -1 : 1) * _multi[nd] * _horizon_radius[nd] / _origin_vec.norm() *
48 1546152 : (_stress[nd] * _shape2[nd].inverse()).row(_component) * _origin_vec *
49 773076 : _node_vol[1 - nd] * _bond_status;
50 193269 : }
51 :
52 : void
53 193269 : HorizonStabilizedFormIISmallStrainMechanicsNOSPD::computeNonlocalResidual()
54 : {
55 579807 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
56 : {
57 : // calculation of residual contribution to current_node's neighbors
58 386538 : std::vector<dof_id_type> ivardofs(_nnodes);
59 386538 : ivardofs[0] = _ivardofs[nd];
60 386538 : std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
61 386538 : std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
62 :
63 : dof_id_type nb_index =
64 386538 : std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - nd)) -
65 386538 : neighbors.begin();
66 : std::vector<dof_id_type> dg_neighbors =
67 386538 : _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(nd), nb_index);
68 :
69 : RealGradient origin_vec_nb;
70 : Real node_vol_nb;
71 :
72 5836674 : for (unsigned int nb = 0; nb < dg_neighbors.size(); ++nb)
73 10513734 : if (neighbors[dg_neighbors[nb]] != _current_elem->node_id(1 - nd) &&
74 5063598 : _bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb]])) > 0.5)
75 : {
76 5063598 : ivardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])
77 5063598 : ->dof_number(_sys.number(), _var.number(), 0);
78 5063598 : origin_vec_nb = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb]]) -
79 5063598 : _pdmesh.getNodeCoord(_current_elem->node_id(nd));
80 5063598 : node_vol_nb = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb]]);
81 :
82 15190794 : for (unsigned int i = 0; i < _nnodes; ++i)
83 15190794 : _local_re(i) = (i == 0 ? -1 : 1) * _multi[nd] * _horizon_radius[nd] /
84 10127196 : origin_vec_nb.norm() *
85 10127196 : (_stress[nd] * _shape2[nd].inverse()).row(_component) * origin_vec_nb *
86 10127196 : node_vol_nb * _bond_status;
87 :
88 : // cache the residual contribution
89 5063598 : addResiduals(_assembly, _local_re, ivardofs, _var.scalingFactor());
90 :
91 5063598 : if (_has_save_in)
92 : {
93 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
94 0 : for (unsigned int i = 0; i < _save_in.size(); ++i)
95 : {
96 0 : std::vector<dof_id_type> save_in_dofs(_nnodes);
97 0 : save_in_dofs[0] = _current_elem->node_ptr(nd)->dof_number(
98 : _save_in[i]->sys().number(), _save_in[i]->number(), 0);
99 0 : save_in_dofs[1] =
100 0 : _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])
101 0 : ->dof_number(_save_in[i]->sys().number(), _save_in[i]->number(), 0);
102 :
103 0 : _save_in[i]->sys().solution().add_vector(_local_re, save_in_dofs);
104 0 : }
105 : }
106 : }
107 386538 : }
108 193269 : }
109 :
110 : void
111 2841 : HorizonStabilizedFormIISmallStrainMechanicsNOSPD::computeLocalJacobian()
112 : {
113 8523 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
114 : {
115 5682 : std::vector<dof_id_type> ivardofs(_nnodes);
116 5682 : ivardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), _var.number(), 0);
117 5682 : std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
118 5682 : std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
119 :
120 : dof_id_type nb_index =
121 5682 : std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - nd)) -
122 5682 : neighbors.begin();
123 : std::vector<dof_id_type> dg_neighbors =
124 5682 : _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(nd), nb_index);
125 :
126 : RealGradient origin_vec_nb;
127 : Real node_vol_nb;
128 :
129 83994 : for (unsigned int nb = 0; nb < dg_neighbors.size(); ++nb)
130 78312 : if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb]])) > 0.5)
131 : {
132 78312 : ivardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])
133 78312 : ->dof_number(_sys.number(), _var.number(), 0);
134 78312 : origin_vec_nb = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb]]) -
135 78312 : _pdmesh.getNodeCoord(_current_elem->node_id(nd));
136 78312 : node_vol_nb = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb]]);
137 :
138 234936 : for (unsigned int i = 0; i < _nnodes; ++i)
139 469872 : for (unsigned int j = 0; j < _nnodes; ++j)
140 313248 : _local_ke(i, j) =
141 626496 : (i == 0 ? -1 : 1) * (j == 0 ? 1 : 0) * _multi[nd] * _horizon_radius[nd] /
142 313248 : origin_vec_nb.norm() *
143 626496 : (computeDSDU(_component, nd) * _shape2[nd].inverse()).row(_component) *
144 313248 : origin_vec_nb * node_vol_nb * _bond_status;
145 :
146 78312 : addJacobian(_assembly, _local_ke, ivardofs, ivardofs, _var.scalingFactor());
147 : }
148 5682 : }
149 : _local_ke.zero();
150 2841 : }
151 :
152 : void
153 2841 : HorizonStabilizedFormIISmallStrainMechanicsNOSPD::computeNonlocalJacobian()
154 : {
155 : // includes dTi/dUj and dTj/dUi contributions
156 8523 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
157 : {
158 : // calculation of jacobian contribution to current_node's neighbors
159 5682 : std::vector<dof_id_type> ivardofs(_nnodes), jvardofs(_nnodes);
160 5682 : ivardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), _var.number(), 0);
161 5682 : jvardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), _var.number(), 0);
162 5682 : std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
163 5682 : std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
164 :
165 : dof_id_type nb_index =
166 5682 : std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - nd)) -
167 5682 : neighbors.begin();
168 : std::vector<dof_id_type> dg_neighbors =
169 5682 : _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(nd), nb_index);
170 :
171 : RealGradient origin_vec_nb1;
172 : Real node_vol_nb1;
173 :
174 83994 : for (unsigned int nb1 = 0; nb1 < dg_neighbors.size(); ++nb1)
175 78312 : if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb1]])) > 0.5)
176 : {
177 78312 : ivardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb1]])
178 78312 : ->dof_number(_sys.number(), _var.number(), 0);
179 78312 : origin_vec_nb1 = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb1]]) -
180 78312 : _pdmesh.getNodeCoord(_current_elem->node_id(nd));
181 78312 : node_vol_nb1 = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb1]]);
182 :
183 : Real vol_nb2;
184 : RealGradient origin_vec_nb2;
185 78312 : RankTwoTensor dFdUk, dPxdUkx;
186 :
187 1230396 : for (unsigned int nb2 = 0; nb2 < dg_neighbors.size(); ++nb2)
188 1152084 : if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb2]])) > 0.5)
189 : {
190 1152084 : jvardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb2]])
191 1152084 : ->dof_number(_sys.number(), _var.number(), 0);
192 1152084 : vol_nb2 = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb2]]);
193 :
194 1152084 : origin_vec_nb2 = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb2]]) -
195 2304168 : _pdmesh.getNodeCoord(_current_elem->node_id(nd));
196 :
197 : dFdUk.zero();
198 4232448 : for (unsigned int i = 0; i < _dim; ++i)
199 3080364 : dFdUk(_component, i) =
200 3080364 : _horizon_radius[nd] / origin_vec_nb2.norm() * origin_vec_nb2(i) * vol_nb2;
201 :
202 1152084 : dFdUk *= _shape2[nd].inverse();
203 2304168 : dPxdUkx = _Jacobian_mult[nd] * 0.5 * (dFdUk.transpose() + dFdUk);
204 :
205 3456252 : for (unsigned int i = 0; i < _nnodes; ++i)
206 6912504 : for (unsigned int j = 0; j < _nnodes; ++j)
207 6912504 : _local_ke(i, j) = (i == 0 ? -1 : 1) * (j == 0 ? 0 : 1) * _multi[nd] *
208 4608336 : _horizon_radius[nd] / origin_vec_nb1.norm() *
209 9216672 : (dPxdUkx * _shape2[nd].inverse()).row(_component) *
210 4608336 : origin_vec_nb1 * node_vol_nb1 * _bond_status;
211 :
212 1152084 : addJacobian(_assembly, _local_ke, ivardofs, jvardofs, _var.scalingFactor());
213 :
214 1152084 : if (_has_diag_save_in)
215 : {
216 0 : unsigned int rows = _nnodes;
217 0 : DenseVector<Real> diag(rows);
218 0 : for (unsigned int i = 0; i < rows; ++i)
219 0 : diag(i) = _local_ke(i, i);
220 :
221 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
222 0 : for (unsigned int i = 0; i < _diag_save_in.size(); ++i)
223 : {
224 0 : std::vector<dof_id_type> diag_save_in_dofs(2);
225 0 : diag_save_in_dofs[0] = _current_elem->node_ptr(nd)->dof_number(
226 : _diag_save_in[i]->sys().number(), _diag_save_in[i]->number(), 0);
227 0 : diag_save_in_dofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb2]])
228 0 : ->dof_number(_diag_save_in[i]->sys().number(),
229 : _diag_save_in[i]->number(),
230 : 0);
231 :
232 0 : _diag_save_in[i]->sys().solution().add_vector(diag, diag_save_in_dofs);
233 0 : }
234 : }
235 : }
236 : }
237 5682 : }
238 2841 : }
239 :
240 : void
241 5094 : HorizonStabilizedFormIISmallStrainMechanicsNOSPD::computeLocalOffDiagJacobian(
242 : unsigned int jvar_num, unsigned int coupled_component)
243 : {
244 5094 : if (_temp_coupled && jvar_num == _temp_var->number()) // temperature is coupled
245 : {
246 588 : std::vector<RankTwoTensor> dSdT(_nnodes);
247 1764 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
248 2352 : for (unsigned int es = 0; es < _deigenstrain_dT.size(); ++es)
249 1176 : dSdT[nd] = -_Jacobian_mult[nd] * (*_deigenstrain_dT[es])[nd];
250 :
251 1764 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
252 : {
253 1176 : std::vector<dof_id_type> ivardofs(_nnodes), jvardofs(_nnodes);
254 1176 : ivardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), _var.number(), 0);
255 1176 : jvardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), jvar_num, 0);
256 1176 : std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
257 1176 : std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
258 :
259 : dof_id_type nb_index =
260 1176 : std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - nd)) -
261 1176 : neighbors.begin();
262 : std::vector<dof_id_type> dg_neighbors =
263 1176 : _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(nd), nb_index);
264 :
265 : RealGradient origin_vec_nb;
266 : Real node_vol_nb;
267 :
268 15888 : for (unsigned int nb = 0; nb < dg_neighbors.size(); ++nb)
269 14712 : if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb]])) > 0.5)
270 : {
271 14712 : ivardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])
272 14712 : ->dof_number(_sys.number(), _var.number(), 0);
273 14712 : jvardofs[1] =
274 14712 : _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])->dof_number(_sys.number(), jvar_num, 0);
275 14712 : origin_vec_nb = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb]]) -
276 14712 : _pdmesh.getNodeCoord(_current_elem->node_id(nd));
277 14712 : node_vol_nb = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb]]);
278 :
279 44136 : for (unsigned int i = 0; i < _nnodes; ++i)
280 88272 : for (unsigned int j = 0; j < _nnodes; ++j)
281 147120 : _local_ke(i, j) = (i == 0 ? -1 : 1) * (j == 0 ? 1 : 0) * _multi[nd] *
282 58848 : _horizon_radius[nd] / origin_vec_nb.norm() *
283 58848 : (dSdT[nd] * _shape2[nd].inverse()).row(_component) * origin_vec_nb *
284 58848 : node_vol_nb * _bond_status;
285 :
286 14712 : addJacobian(_assembly, _local_ke, ivardofs, jvardofs, _var.scalingFactor());
287 : }
288 1176 : }
289 : _local_ke.zero();
290 588 : }
291 4506 : else if (_out_of_plane_strain_coupled &&
292 0 : jvar_num == _out_of_plane_strain_var
293 : ->number()) // weak plane stress case, out_of_plane_strain is coupled
294 : {
295 0 : std::vector<RankTwoTensor> dSdE33(_nnodes);
296 0 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
297 0 : for (unsigned int i = 0; i < 3; ++i)
298 0 : for (unsigned int j = 0; j < 3; ++j)
299 0 : dSdE33[nd](i, j) = _Jacobian_mult[nd](i, j, 2, 2);
300 :
301 0 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
302 : {
303 0 : std::vector<dof_id_type> ivardofs(_nnodes), jvardofs(_nnodes);
304 0 : ivardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), _var.number(), 0);
305 0 : jvardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), jvar_num, 0);
306 0 : std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
307 0 : std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
308 :
309 : dof_id_type nb_index =
310 0 : std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - nd)) -
311 0 : neighbors.begin();
312 : std::vector<dof_id_type> dg_neighbors =
313 0 : _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(nd), nb_index);
314 :
315 : RealGradient origin_vec_nb;
316 : Real node_vol_nb;
317 :
318 0 : for (unsigned int nb = 0; nb < dg_neighbors.size(); ++nb)
319 0 : if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb]])) > 0.5)
320 : {
321 0 : ivardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])
322 0 : ->dof_number(_sys.number(), _var.number(), 0);
323 0 : jvardofs[1] =
324 0 : _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])->dof_number(_sys.number(), jvar_num, 0);
325 0 : origin_vec_nb = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb]]) -
326 0 : _pdmesh.getNodeCoord(_current_elem->node_id(nd));
327 0 : node_vol_nb = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb]]);
328 :
329 0 : for (unsigned int i = 0; i < _nnodes; ++i)
330 0 : for (unsigned int j = 0; j < _nnodes; ++j)
331 0 : _local_ke(i, j) = (i == 0 ? -1 : 1) * (j == 0 ? 1 : 0) * _multi[nd] *
332 0 : _horizon_radius[nd] / origin_vec_nb.norm() *
333 0 : (dSdE33[nd] * _shape2[nd].inverse()).row(_component) *
334 0 : origin_vec_nb * node_vol_nb * _bond_status;
335 :
336 0 : addJacobian(_assembly, _local_ke, ivardofs, jvardofs, _var.scalingFactor());
337 : }
338 0 : }
339 : _local_ke.zero();
340 0 : }
341 : else // off-diagonal Jacobian with respect to other displacement variables
342 : {
343 13518 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
344 : {
345 9012 : std::vector<dof_id_type> ivardofs(_nnodes), jvardofs(_nnodes);
346 9012 : ivardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), _var.number(), 0);
347 9012 : jvardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), jvar_num, 0);
348 9012 : std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
349 9012 : std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
350 :
351 : dof_id_type nb_index =
352 9012 : std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - nd)) -
353 9012 : neighbors.begin();
354 : std::vector<dof_id_type> dg_neighbors =
355 9012 : _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(nd), nb_index);
356 :
357 : RealGradient origin_vec_nb;
358 : Real node_vol_nb;
359 :
360 136212 : for (unsigned int nb = 0; nb < dg_neighbors.size(); ++nb)
361 127200 : if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb]])) > 0.5)
362 : {
363 127200 : ivardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])
364 127200 : ->dof_number(_sys.number(), _var.number(), 0);
365 127200 : jvardofs[1] =
366 127200 : _pdmesh.nodePtr(neighbors[dg_neighbors[nb]])->dof_number(_sys.number(), jvar_num, 0);
367 127200 : origin_vec_nb = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb]]) -
368 127200 : _pdmesh.getNodeCoord(_current_elem->node_id(nd));
369 127200 : node_vol_nb = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb]]);
370 :
371 381600 : for (unsigned int i = 0; i < _nnodes; ++i)
372 763200 : for (unsigned int j = 0; j < _nnodes; ++j)
373 508800 : _local_ke(i, j) =
374 1017600 : (i == 0 ? -1 : 1) * (j == 0 ? 1 : 0) * _multi[nd] * _horizon_radius[nd] /
375 508800 : origin_vec_nb.norm() *
376 1017600 : (computeDSDU(coupled_component, nd) * _shape2[nd].inverse()).row(_component) *
377 508800 : origin_vec_nb * node_vol_nb * _bond_status;
378 :
379 127200 : addJacobian(_assembly, _local_ke, ivardofs, jvardofs, _var.scalingFactor());
380 : }
381 9012 : }
382 : _local_ke.zero();
383 : }
384 5094 : }
385 :
386 : void
387 5094 : HorizonStabilizedFormIISmallStrainMechanicsNOSPD::computePDNonlocalOffDiagJacobian(
388 : unsigned int jvar_num, unsigned int coupled_component)
389 : {
390 5094 : if (_temp_coupled && jvar_num == _temp_var->number())
391 : {
392 : // no nonlocal contribution from temperature
393 : }
394 4506 : else if (_out_of_plane_strain_coupled && jvar_num == _out_of_plane_strain_var->number())
395 : {
396 : // no nonlocal contribution from out of plane strain
397 : }
398 : else
399 : {
400 13518 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
401 : {
402 : // calculation of jacobian contribution to current_node's neighbors
403 : // NOT including the contribution to nodes i and j, which is considered as local off-diagonal
404 9012 : std::vector<dof_id_type> ivardofs(_nnodes), jvardofs(_nnodes);
405 9012 : ivardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), _var.number(), 0);
406 9012 : jvardofs[0] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), jvar_num, 0);
407 9012 : std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
408 9012 : std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
409 :
410 : dof_id_type nb_index =
411 9012 : std::find(neighbors.begin(), neighbors.end(), _current_elem->node_id(1 - nd)) -
412 9012 : neighbors.begin();
413 : std::vector<dof_id_type> dg_neighbors =
414 9012 : _pdmesh.getBondDeformationGradientNeighbors(_current_elem->node_id(nd), nb_index);
415 :
416 : RealGradient origin_vec_nb1;
417 : Real node_vol_nb1;
418 :
419 136212 : for (unsigned int nb1 = 0; nb1 < dg_neighbors.size(); ++nb1)
420 127200 : if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb1]])) > 0.5)
421 : {
422 127200 : ivardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb1]])
423 127200 : ->dof_number(_sys.number(), _var.number(), 0);
424 127200 : origin_vec_nb1 = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb1]]) -
425 127200 : _pdmesh.getNodeCoord(_current_elem->node_id(nd));
426 127200 : node_vol_nb1 = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb1]]);
427 :
428 : Real vol_nb2;
429 : RealGradient origin_vec_nb2;
430 127200 : RankTwoTensor dFdUk, dPxdUky;
431 :
432 2055480 : for (unsigned int nb2 = 0; nb2 < dg_neighbors.size(); ++nb2)
433 1928280 : if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[dg_neighbors[nb2]])) >
434 : 0.5)
435 : {
436 1928280 : jvardofs[1] = _pdmesh.nodePtr(neighbors[dg_neighbors[nb2]])
437 1928280 : ->dof_number(_sys.number(), jvar_num, 0);
438 1928280 : vol_nb2 = _pdmesh.getNodeVolume(neighbors[dg_neighbors[nb2]]);
439 :
440 1928280 : origin_vec_nb2 = _pdmesh.getNodeCoord(neighbors[dg_neighbors[nb2]]) -
441 3856560 : _pdmesh.getNodeCoord(_current_elem->node_id(nd));
442 :
443 : dFdUk.zero();
444 7337232 : for (unsigned int i = 0; i < _dim; ++i)
445 5408952 : dFdUk(coupled_component, i) =
446 5408952 : _horizon_radius[nd] / origin_vec_nb2.norm() * origin_vec_nb2(i) * vol_nb2;
447 :
448 1928280 : dFdUk *= _shape2[nd].inverse();
449 3856560 : dPxdUky = _Jacobian_mult[nd] * 0.5 * (dFdUk.transpose() + dFdUk);
450 :
451 5784840 : for (unsigned int i = 0; i < _nnodes; ++i)
452 11569680 : for (unsigned int j = 0; j < _nnodes; ++j)
453 11569680 : _local_ke(i, j) = (i == 0 ? -1 : 1) * (j == 0 ? 0 : 1) * _multi[nd] *
454 7713120 : _horizon_radius[nd] / origin_vec_nb1.norm() *
455 15426240 : (dPxdUky * _shape2[nd].inverse()).row(_component) *
456 7713120 : origin_vec_nb1 * node_vol_nb1 * _bond_status;
457 :
458 1928280 : addJacobian(_assembly, _local_ke, ivardofs, jvardofs, _var.scalingFactor());
459 : }
460 : }
461 9012 : }
462 : }
463 5094 : }
|
