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 "GeneralizedPlaneStrainOffDiagOSPD.h"
11 : #include "MooseVariableScalar.h"
12 : #include "PeridynamicsMesh.h"
13 :
14 : registerMooseObject("PeridynamicsApp", GeneralizedPlaneStrainOffDiagOSPD);
15 :
16 : InputParameters
17 134 : GeneralizedPlaneStrainOffDiagOSPD::validParams()
18 : {
19 134 : InputParameters params = MechanicsBasePD::validParams();
20 134 : params.addClassDescription(
21 : "Class for calculating the off-diagonal Jacobian corresponding to "
22 : "coupling between displacements (or temperature) and the scalar out-of-plane strain for "
23 : "the generalized plane strain using the OSPD formulation");
24 :
25 268 : params.addCoupledVar("scalar_out_of_plane_strain",
26 : "Scalar variable for strain in the out-of-plane direction");
27 :
28 134 : return params;
29 0 : }
30 :
31 98 : GeneralizedPlaneStrainOffDiagOSPD::GeneralizedPlaneStrainOffDiagOSPD(
32 98 : const InputParameters & parameters)
33 : : MechanicsBasePD(parameters),
34 98 : _bond_local_dfdE(getMaterialProperty<Real>("bond_local_dfdE")),
35 196 : _bond_nonlocal_dfdE(getMaterialProperty<Real>("bond_nonlocal_dfdE")),
36 196 : _alpha(getMaterialProperty<Real>("thermal_expansion_coeff")),
37 196 : _Cijkl(getMaterialProperty<RankFourTensor>("elasticity_tensor")),
38 196 : _scalar_out_of_plane_strain_var_num(coupledScalar("scalar_out_of_plane_strain"))
39 : {
40 : // Consistency check
41 98 : if (_disp_var.size() != 2)
42 0 : mooseError("GeneralizedPlaneStrain only works for two dimensional case!");
43 98 : }
44 :
45 : void
46 25320 : GeneralizedPlaneStrainOffDiagOSPD::computeOffDiagJacobianScalar(unsigned int jvar_num)
47 : {
48 25320 : if (jvar_num == _scalar_out_of_plane_strain_var_num)
49 : {
50 25320 : prepare();
51 :
52 25320 : if (_var.number() == _disp_var[0]->number())
53 12366 : if (_use_full_jacobian)
54 1176 : computeDispFullOffDiagJacobianScalar(0, jvar_num);
55 : else
56 11190 : computeDispPartialOffDiagJacobianScalar(0, jvar_num);
57 12954 : else if (_var.number() == _disp_var[1]->number())
58 12366 : if (_use_full_jacobian)
59 1176 : computeDispFullOffDiagJacobianScalar(1, jvar_num);
60 : else
61 11190 : computeDispPartialOffDiagJacobianScalar(1, jvar_num);
62 588 : else if (_temp_coupled ? _var.number() == _temp_var->number() : 0)
63 588 : computeTempOffDiagJacobianScalar(jvar_num);
64 : }
65 25320 : }
66 :
67 : void
68 2352 : GeneralizedPlaneStrainOffDiagOSPD::computeDispFullOffDiagJacobianScalar(unsigned int component,
69 : unsigned int jvar_num)
70 : {
71 2352 : prepareMatrixTag(_assembly, _var.number(), jvar_num, _ken);
72 2352 : prepareMatrixTag(_assembly, jvar_num, _var.number(), _kne);
73 2352 : MooseVariableScalar & jvar = _sys.getScalarVariable(_tid, jvar_num);
74 :
75 : // LOCAL jacobian contribution
76 : // fill in the column corresponding to the scalar variable from bond ij
77 7056 : for (unsigned int i = 0; i < _nnodes; ++i)
78 9408 : for (unsigned int j = 0; j < jvar.order(); ++j)
79 4704 : _ken(i, j) +=
80 7056 : (i == j ? -1 : 1) * _current_unit_vec(component) * _bond_local_dfdE[0] * _bond_status;
81 :
82 : // NONLOCAL jacobian contribution
83 2352 : std::vector<RankTwoTensor> shape(_nnodes), dgrad(_nnodes);
84 :
85 : // for off-diagonal low components
86 2352 : if (_bond_status > 0.5)
87 : {
88 6768 : for (unsigned int nd = 0; nd < _nnodes; nd++)
89 : {
90 4512 : if (_dim == 2)
91 4512 : shape[nd](2, 2) = dgrad[nd](2, 2) = 1.0;
92 : // calculation of jacobian contribution to current node's neighbors
93 4512 : std::vector<dof_id_type> ivardofs(_nnodes);
94 4512 : ivardofs[nd] = _ivardofs[nd];
95 4512 : std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_elem->node_id(nd));
96 4512 : std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
97 :
98 : Real vol_nb, weight_nb;
99 : RealGradient origin_vec_nb, current_vec_nb;
100 :
101 61168 : for (unsigned int nb = 0; nb < neighbors.size(); nb++)
102 56656 : if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[nb])) > 0.5)
103 : {
104 55424 : ivardofs[1 - nd] =
105 55424 : _pdmesh.nodePtr(neighbors[nb])->dof_number(_sys.number(), _var.number(), 0);
106 55424 : vol_nb = _pdmesh.getNodeVolume(neighbors[nb]);
107 :
108 : // obtain bond nb's origin length and current orientation
109 55424 : origin_vec_nb = _pdmesh.getNodeCoord(neighbors[nb]) -
110 55424 : _pdmesh.getNodeCoord(_current_elem->node_id(nd));
111 :
112 166272 : for (unsigned int k = 0; k < _dim; k++)
113 110848 : current_vec_nb(k) = origin_vec_nb(k) +
114 221696 : _disp_var[k]->getNodalValue(*_pdmesh.nodePtr(neighbors[nb])) -
115 110848 : _disp_var[k]->getNodalValue(*_current_elem->node_ptr(nd));
116 :
117 55424 : weight_nb = _horizon_radius[nd] / origin_vec_nb.norm();
118 : // prepare shape tensor and deformation gradient tensor for current node
119 166272 : for (unsigned int k = 0; k < _dim; k++)
120 332544 : for (unsigned int l = 0; l < _dim; l++)
121 : {
122 221696 : shape[nd](k, l) += weight_nb * origin_vec_nb(k) * origin_vec_nb(l) * vol_nb;
123 221696 : dgrad[nd](k, l) += weight_nb * current_vec_nb(k) * origin_vec_nb(l) * vol_nb;
124 : }
125 :
126 : // cache the nonlocal jacobian contribution
127 55424 : _local_ke.resize(_ken.m(), _ken.n());
128 : _local_ke.zero();
129 83168 : _local_ke(0, 0) = (nd == 0 ? -1 : 1) * current_vec_nb(component) / current_vec_nb.norm() *
130 55424 : _bond_nonlocal_dfdE[nd] / origin_vec_nb.norm() * vol_nb * _bond_status;
131 83168 : _local_ke(1, 0) = (nd == 0 ? 1 : -1) * current_vec_nb(component) / current_vec_nb.norm() *
132 55424 : _bond_nonlocal_dfdE[nd] / origin_vec_nb.norm() * vol_nb * _bond_status;
133 :
134 55424 : addJacobian(_assembly, _local_ke, ivardofs, jvar.dofIndices(), _var.scalingFactor());
135 : }
136 :
137 : // finalize the shape tensor and deformation gradient tensor for node_i
138 4512 : if (shape[nd].det() == 0.)
139 0 : mooseError("Singular shape tensor is detected in GeneralizedPlaneStrainOffDiagOSPD! Use "
140 : "SingularShapeTensorEliminatorUserObjectPD to avoid singular shape tensor");
141 :
142 4512 : shape[nd] = shape[nd].inverse();
143 4512 : dgrad[nd] = dgrad[nd] * shape[nd];
144 4512 : }
145 : }
146 :
147 : // off-diagonal jacobian entries on the row
148 : Real dEidUi =
149 2352 : -_node_vol[1] * _horizon_radius[0] / _origin_vec.norm() *
150 2352 : (_Cijkl[0](2, 2, 0, 0) * (_origin_vec(0) * shape[0](0, 0) + _origin_vec(1) * shape[0](1, 0)) *
151 2352 : dgrad[0](component, 0) +
152 2352 : _Cijkl[0](2, 2, 1, 1) * (_origin_vec(0) * shape[0](0, 1) + _origin_vec(1) * shape[0](1, 1)) *
153 2352 : dgrad[0](component, 1));
154 : Real dEjdUj =
155 2352 : _node_vol[0] * _horizon_radius[1] / _origin_vec.norm() *
156 2352 : (_Cijkl[0](2, 2, 0, 0) * (_origin_vec(0) * shape[1](0, 0) + _origin_vec(1) * shape[1](1, 0)) *
157 2352 : dgrad[1](component, 0) +
158 2352 : _Cijkl[0](2, 2, 1, 1) * (_origin_vec(0) * shape[1](0, 1) + _origin_vec(1) * shape[1](1, 1)) *
159 2352 : dgrad[1](component, 1));
160 :
161 : // fill in the row corresponding to the scalar variable
162 2352 : _kne(0, 0) += (dEidUi * _node_vol[0] - dEjdUj * _node_vol[1]) * _bond_status; // node i
163 2352 : _kne(0, 1) += (dEjdUj * _node_vol[1] - dEidUi * _node_vol[0]) * _bond_status; // node j
164 :
165 2352 : accumulateTaggedLocalMatrix(_assembly, _var.number(), jvar_num, _ken);
166 2352 : accumulateTaggedLocalMatrix(_assembly, jvar_num, _var.number(), _kne);
167 2352 : }
168 :
169 : void
170 22380 : GeneralizedPlaneStrainOffDiagOSPD::computeDispPartialOffDiagJacobianScalar(unsigned int component,
171 : unsigned int jvar_num)
172 : {
173 22380 : prepareMatrixTag(_assembly, _var.number(), jvar_num, _ken);
174 22380 : prepareMatrixTag(_assembly, jvar_num, _var.number(), _kne);
175 22380 : MooseVariableScalar & jvar = _sys.getScalarVariable(_tid, jvar_num);
176 :
177 : // fill in the column corresponding to the scalar variable from bond ij
178 67140 : for (unsigned int i = 0; i < _nnodes; ++i)
179 89520 : for (unsigned int j = 0; j < jvar.order(); ++j)
180 : {
181 44760 : _ken(i, j) +=
182 67140 : (i == j ? -1 : 1) * _current_unit_vec(component) * _bond_local_dfdE[0] * _bond_status;
183 44760 : _kne(j, i) +=
184 44760 : (i == j ? -1 : 1) * _current_unit_vec(component) * _bond_local_dfdE[0] * _bond_status;
185 : }
186 :
187 22380 : accumulateTaggedLocalMatrix(_assembly, _var.number(), jvar_num, _ken);
188 22380 : accumulateTaggedLocalMatrix(_assembly, jvar_num, _var.number(), _kne);
189 22380 : }
190 :
191 : void
192 588 : GeneralizedPlaneStrainOffDiagOSPD::computeTempOffDiagJacobianScalar(unsigned int jvar_num)
193 : {
194 : // off-diagonal jacobian entries on the row
195 588 : prepareMatrixTag(_assembly, jvar_num, _var.number(), _kne);
196 :
197 : // calculate number of active neighbors for node i and j
198 588 : std::vector<unsigned int> active_neighbors(_nnodes, 0);
199 1764 : for (unsigned int nd = 0; nd < _nnodes; nd++)
200 : {
201 1176 : std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_elem->node_id(nd));
202 15888 : for (unsigned int nb = 0; nb < bonds.size(); ++nb)
203 14712 : if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[nb])) > 0.5)
204 14164 : active_neighbors[nd]++;
205 :
206 1176 : if (active_neighbors[nd] == 0) // avoid dividing by zero
207 20 : active_neighbors[nd] = 1;
208 1176 : }
209 :
210 : // one-way coupling between the strain_zz and temperature. fill in the row corresponding to the
211 : // scalar variable strain_zz
212 588 : _kne(0, 0) += -_alpha[0] *
213 588 : (_Cijkl[0](2, 2, 0, 0) + _Cijkl[0](2, 2, 1, 1) + _Cijkl[0](2, 2, 2, 2)) *
214 588 : _node_vol[0] / active_neighbors[0] * _bond_status; // node i
215 588 : _kne(0, 1) += -_alpha[0] *
216 588 : (_Cijkl[0](2, 2, 0, 0) + _Cijkl[0](2, 2, 1, 1) + _Cijkl[0](2, 2, 2, 2)) *
217 588 : _node_vol[1] / active_neighbors[1] * _bond_status; // node j
218 :
219 588 : accumulateTaggedLocalMatrix(_assembly, jvar_num, _var.number(), _kne);
220 588 : }
|
