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 "MechanicsOSPD.h"
11 : #include "PeridynamicsMesh.h"
12 :
13 : registerMooseObject("PeridynamicsApp", MechanicsOSPD);
14 :
15 : InputParameters
16 282 : MechanicsOSPD::validParams()
17 : {
18 282 : InputParameters params = MechanicsBasePD::validParams();
19 282 : params.addClassDescription(
20 : "Class for calculating the residual and Jacobian for the ordinary state-based "
21 : "peridynamic mechanics formulation");
22 :
23 564 : 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 282 : return params;
28 0 : }
29 :
30 200 : MechanicsOSPD::MechanicsOSPD(const InputParameters & parameters)
31 : : MechanicsBasePD(parameters),
32 200 : _bond_local_force(getMaterialProperty<Real>("bond_local_force")),
33 400 : _bond_nonlocal_force(getMaterialProperty<Real>("bond_nonlocal_force")),
34 400 : _bond_local_dfdU(getMaterialProperty<Real>("bond_dfdU")),
35 400 : _bond_nonlocal_dfdU(getMaterialProperty<Real>("bond_nonlocal_dfdU")),
36 400 : _bond_local_dfdT(getMaterialProperty<Real>("bond_dfdT")),
37 400 : _bond_nonlocal_dfdT(getMaterialProperty<Real>("bond_nonlocal_dfdT")),
38 600 : _component(getParam<unsigned int>("component"))
39 : {
40 200 : }
41 :
42 : void
43 883428 : MechanicsOSPD::computeLocalResidual()
44 : {
45 : // H term
46 883428 : _local_re(0) = -_bond_local_force[0] * _current_unit_vec(_component) * _bond_status;
47 883428 : _local_re(1) = -_local_re(0);
48 883428 : }
49 :
50 : void
51 883428 : MechanicsOSPD::computeNonlocalResidual()
52 : {
53 : // P and Q terms
54 2650284 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
55 : {
56 : // calculation of residual contribution to current_node's neighbors
57 1766856 : std::vector<dof_id_type> ivardofs(_nnodes);
58 1766856 : ivardofs[nd] = _ivardofs[nd];
59 1766856 : std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
60 1766856 : std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
61 :
62 : Real vol_nb;
63 : RealGradient origin_vec_nb, current_vec_nb;
64 :
65 28195456 : for (unsigned int nb = 0; nb < neighbors.size(); ++nb)
66 26428600 : if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[nb])) > 0.5)
67 : {
68 26233512 : ivardofs[1 - nd] =
69 26233512 : _pdmesh.nodePtr(neighbors[nb])->dof_number(_sys.number(), _var.number(), 0);
70 26233512 : vol_nb = _pdmesh.getNodeVolume(neighbors[nb]);
71 :
72 26233512 : origin_vec_nb =
73 26233512 : _pdmesh.getNodeCoord(neighbors[nb]) - _pdmesh.getNodeCoord(_current_elem->node_id(nd));
74 :
75 78700536 : for (unsigned int i = 0; i < _dim; ++i)
76 52467024 : current_vec_nb(i) = origin_vec_nb(i) +
77 104934048 : _disp_var[i]->getNodalValue(*_pdmesh.nodePtr(neighbors[nb])) -
78 52467024 : _disp_var[i]->getNodalValue(*_current_elem->node_ptr(nd));
79 :
80 26233512 : current_vec_nb /= current_vec_nb.norm();
81 :
82 39357794 : _local_re(0) = (nd == 0 ? -1 : 1) * _bond_nonlocal_force[nd] * vol_nb /
83 26233512 : origin_vec_nb.norm() * current_vec_nb(_component) * _bond_status;
84 26233512 : _local_re(1) = -_local_re(0);
85 :
86 : // cache the residual contribution to node_i and its neighbor nb using their global dof
87 : // indices
88 26233512 : addResiduals(_assembly, _local_re, ivardofs, _var.scalingFactor());
89 :
90 : // save in the displacement residuals
91 26233512 : 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(2);
97 0 : save_in_dofs[nd] = _current_elem->node_ptr(nd)->dof_number(
98 : _save_in[i]->sys().number(), _save_in[i]->number(), 0);
99 0 : save_in_dofs[1 - nd] =
100 0 : _pdmesh.nodePtr(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 1766856 : }
108 883428 : }
109 :
110 : void
111 69976 : MechanicsOSPD::computeLocalJacobian()
112 : {
113 : const Real val =
114 69976 : _current_unit_vec(_component) * _current_unit_vec(_component) * _bond_local_dfdU[0] +
115 139952 : _bond_local_force[0] * (1.0 - _current_unit_vec(_component) * _current_unit_vec(_component)) /
116 69976 : _current_vec.norm();
117 :
118 209928 : for (unsigned int i = 0; i < _nnodes; ++i)
119 419856 : for (unsigned int j = 0; j < _nnodes; ++j)
120 419856 : _local_ke(i, j) += (i == j ? 1 : -1) * val * _bond_status;
121 69976 : }
122 :
123 : void
124 72328 : MechanicsOSPD::computeLocalOffDiagJacobian(unsigned int jvar_num, unsigned int coupled_component)
125 : {
126 72328 : if (_temp_coupled && jvar_num == _temp_var->number())
127 : {
128 25836 : for (unsigned int i = 0; i < _nnodes; ++i)
129 51672 : for (unsigned int j = 0; j < _nnodes; ++j)
130 34448 : _local_ke(i, j) +=
131 51672 : (i == 1 ? 1 : -1) * _current_unit_vec(_component) * _bond_local_dfdT[0] * _bond_status;
132 : }
133 : else
134 : {
135 : const Real val =
136 63716 : _current_unit_vec(_component) * _current_unit_vec(coupled_component) * _bond_local_dfdU[0] -
137 63716 : _bond_local_force[0] * _current_unit_vec(_component) *
138 63716 : _current_unit_vec(coupled_component) / _current_vec.norm();
139 :
140 191148 : for (unsigned int i = 0; i < _nnodes; ++i)
141 382296 : for (unsigned int j = 0; j < _nnodes; ++j)
142 382296 : _local_ke(i, j) += (i == j ? 1 : -1) * val * _bond_status;
143 : }
144 72328 : }
145 :
146 : void
147 4704 : MechanicsOSPD::computeNonlocalJacobian()
148 : {
149 14112 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
150 : {
151 : // calculation of jacobian contribution to current_node's neighbors
152 9408 : std::vector<dof_id_type> ivardofs(_nnodes);
153 9408 : ivardofs[nd] = _ivardofs[nd];
154 9408 : std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
155 9408 : std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
156 :
157 : Real vol_nb;
158 : RealGradient origin_vec_nb, current_vec_nb;
159 :
160 127104 : for (unsigned int nb = 0; nb < neighbors.size(); ++nb)
161 117696 : if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[nb])) > 0.5)
162 : {
163 113312 : ivardofs[1 - nd] =
164 113312 : _pdmesh.nodePtr(neighbors[nb])->dof_number(_sys.number(), _var.number(), 0);
165 113312 : vol_nb = _pdmesh.getNodeVolume(neighbors[nb]);
166 :
167 113312 : origin_vec_nb =
168 113312 : _pdmesh.getNodeCoord(neighbors[nb]) - _pdmesh.getNodeCoord(_current_elem->node_id(nd));
169 :
170 339936 : for (unsigned int k = 0; k < _dim; ++k)
171 226624 : current_vec_nb(k) = origin_vec_nb(k) +
172 453248 : _disp_var[k]->getNodalValue(*_pdmesh.nodePtr(neighbors[nb])) -
173 226624 : _disp_var[k]->getNodalValue(*_current_elem->node_ptr(nd));
174 :
175 113312 : current_vec_nb /= current_vec_nb.norm();
176 :
177 113312 : const Real val = (nd == 0 ? 1 : -1) * current_vec_nb(_component) *
178 113312 : _current_unit_vec(_component) * _bond_nonlocal_dfdU[nd];
179 :
180 : _local_ke.zero();
181 339936 : for (unsigned int i = 0; i < _nnodes; ++i)
182 679872 : for (unsigned int j = 0; j < _nnodes; ++j)
183 453248 : _local_ke(i, j) +=
184 679872 : (i == j ? 1 : -1) * val / origin_vec_nb.norm() * vol_nb * _bond_status;
185 :
186 113312 : addJacobian(_assembly, _local_ke, ivardofs, _ivardofs, _var.scalingFactor());
187 :
188 113312 : if (_has_diag_save_in)
189 : {
190 0 : unsigned int rows = _nnodes;
191 0 : DenseVector<Real> diag(rows);
192 0 : for (unsigned int i = 0; i < rows; ++i)
193 0 : diag(i) = _local_ke(i, i);
194 :
195 0 : diag(1 - nd) = 0;
196 :
197 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
198 0 : for (unsigned int i = 0; i < _diag_save_in.size(); ++i)
199 : {
200 0 : std::vector<dof_id_type> diag_save_in_dofs(2);
201 0 : diag_save_in_dofs[nd] = _current_elem->node_ptr(nd)->dof_number(
202 : _diag_save_in[i]->sys().number(), _diag_save_in[i]->number(), 0);
203 0 : diag_save_in_dofs[1 - nd] =
204 0 : _pdmesh.nodePtr(neighbors[nb])
205 0 : ->dof_number(_diag_save_in[i]->sys().number(), _diag_save_in[i]->number(), 0);
206 :
207 0 : _diag_save_in[i]->sys().solution().add_vector(diag, diag_save_in_dofs);
208 0 : }
209 : }
210 :
211 113312 : const Real val2 = _bond_nonlocal_force[nd] *
212 113312 : (1.0 - current_vec_nb(_component) * current_vec_nb(_component)) /
213 113312 : current_vec_nb.norm();
214 :
215 : _local_ke.zero();
216 339936 : for (unsigned int i = 0; i < _nnodes; ++i)
217 679872 : for (unsigned int j = 0; j < _nnodes; ++j)
218 453248 : _local_ke(i, j) +=
219 679872 : (i == j ? 1 : -1) * val2 / origin_vec_nb.norm() * vol_nb * _bond_status;
220 :
221 113312 : addJacobian(_assembly, _local_ke, ivardofs, ivardofs, _var.scalingFactor());
222 :
223 113312 : if (_has_diag_save_in)
224 : {
225 0 : unsigned int rows = _nnodes;
226 0 : DenseVector<Real> diag(rows);
227 0 : for (unsigned int i = 0; i < rows; ++i)
228 0 : diag(i) = _local_ke(i, i);
229 :
230 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
231 0 : for (unsigned int i = 0; i < _diag_save_in.size(); ++i)
232 : {
233 0 : std::vector<dof_id_type> diag_save_in_dofs(2);
234 0 : diag_save_in_dofs[nd] = _current_elem->node_ptr(nd)->dof_number(
235 : _diag_save_in[i]->sys().number(), _diag_save_in[i]->number(), 0);
236 0 : diag_save_in_dofs[1 - nd] =
237 0 : _pdmesh.nodePtr(neighbors[nb])
238 0 : ->dof_number(_diag_save_in[i]->sys().number(), _diag_save_in[i]->number(), 0);
239 :
240 0 : _diag_save_in[i]->sys().solution().add_vector(diag, diag_save_in_dofs);
241 0 : }
242 : }
243 : }
244 9408 : }
245 4704 : }
246 :
247 : void
248 7056 : MechanicsOSPD::computePDNonlocalOffDiagJacobian(unsigned int jvar_num,
249 : unsigned int coupled_component)
250 : {
251 7056 : std::vector<dof_id_type> jvardofs_ij(_nnodes);
252 21168 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
253 14112 : jvardofs_ij[nd] = _current_elem->node_ptr(nd)->dof_number(_sys.number(), jvar_num, 0);
254 :
255 21168 : for (unsigned int nd = 0; nd < _nnodes; ++nd)
256 : {
257 : // calculation of jacobian contribution to current_node's neighbors
258 14112 : std::vector<dof_id_type> ivardofs(_nnodes), jvardofs(_nnodes);
259 14112 : ivardofs[nd] = _ivardofs[nd];
260 14112 : jvardofs[nd] = jvardofs_ij[nd];
261 14112 : std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
262 14112 : std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
263 :
264 : Real vol_nb;
265 : RealGradient origin_vec_nb, current_vec_nb;
266 :
267 190656 : for (unsigned int nb = 0; nb < neighbors.size(); ++nb)
268 176544 : if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[nb])) > 0.5)
269 : {
270 169968 : ivardofs[1 - nd] =
271 169968 : _pdmesh.nodePtr(neighbors[nb])->dof_number(_sys.number(), _var.number(), 0);
272 169968 : jvardofs[1 - nd] = _pdmesh.nodePtr(neighbors[nb])->dof_number(_sys.number(), jvar_num, 0);
273 169968 : vol_nb = _pdmesh.getNodeVolume(neighbors[nb]);
274 :
275 169968 : origin_vec_nb =
276 169968 : _pdmesh.getNodeCoord(neighbors[nb]) - _pdmesh.getNodeCoord(_current_elem->node_id(nd));
277 :
278 509904 : for (unsigned int i = 0; i < _dim; ++i)
279 339936 : current_vec_nb(i) = origin_vec_nb(i) +
280 679872 : _disp_var[i]->getNodalValue(*_pdmesh.nodePtr(neighbors[nb])) -
281 339936 : _disp_var[i]->getNodalValue(*_current_elem->node_ptr(nd));
282 :
283 169968 : current_vec_nb /= current_vec_nb.norm();
284 :
285 : _local_ke.zero();
286 169968 : if (_temp_coupled && jvar_num == _temp_var->number())
287 : {
288 : const Real val =
289 56656 : current_vec_nb(_component) * _bond_nonlocal_dfdT[nd] / origin_vec_nb.norm() * vol_nb;
290 :
291 85536 : _local_ke(0, nd) += (nd == 0 ? -1 : 1) * val * _bond_status;
292 85536 : _local_ke(1, nd) += (nd == 0 ? 1 : -1) * val * _bond_status;
293 : }
294 : else
295 : {
296 113312 : const Real val = (nd == 0 ? 1 : -1) * current_vec_nb(_component) *
297 113312 : _current_unit_vec(coupled_component) * _bond_nonlocal_dfdU[nd] /
298 113312 : origin_vec_nb.norm() * vol_nb;
299 :
300 339936 : for (unsigned int i = 0; i < _nnodes; ++i)
301 679872 : for (unsigned int j = 0; j < _nnodes; ++j)
302 679872 : _local_ke(i, j) += (i == j ? 1 : -1) * val * _bond_status;
303 : }
304 :
305 169968 : addJacobian(_assembly, _local_ke, ivardofs, jvardofs_ij, _var.scalingFactor());
306 : }
307 14112 : }
308 7056 : }
|
