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