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 "NodeFaceConstraint.h"
11 :
12 : // MOOSE includes
13 : #include "Assembly.h"
14 : #include "MooseEnum.h"
15 : #include "MooseMesh.h"
16 : #include "MooseVariableFE.h"
17 : #include "PenetrationLocator.h"
18 : #include "SystemBase.h"
19 :
20 : #include "libmesh/string_to_enum.h"
21 :
22 : InputParameters
23 42999 : NodeFaceConstraint::validParams()
24 : {
25 42999 : MooseEnum orders("FIRST SECOND THIRD FOURTH", "FIRST");
26 42999 : InputParameters params = Constraint::validParams();
27 42999 : params.addParam<BoundaryName>("secondary", "The boundary ID associated with the secondary side");
28 42999 : params.addParam<BoundaryName>("primary", "The boundary ID associated with the primary side");
29 42999 : params.addParam<Real>("tangential_tolerance",
30 : "Tangential distance to extend edges of contact surfaces");
31 42999 : params.addParam<Real>(
32 : "normal_smoothing_distance",
33 : "Distance from edge in parametric coordinates over which to smooth contact normal");
34 42999 : params.addParam<std::string>("normal_smoothing_method",
35 : "Method to use to smooth normals (edge_based|nodal_normal_based)");
36 42999 : params.addParam<MooseEnum>("order", orders, "The finite element order used for projections");
37 :
38 42999 : params.addCoupledVar("primary_variable", "The variable on the primary side of the domain");
39 :
40 85998 : return params;
41 42999 : }
42 :
43 101 : NodeFaceConstraint::NodeFaceConstraint(const InputParameters & parameters)
44 : : Constraint(parameters),
45 : // The secondary side is at nodes (hence passing 'true'). The neighbor side is the primary side
46 : // and it is not at nodes (so passing false)
47 : NeighborCoupleableMooseVariableDependencyIntermediateInterface(this, true, false),
48 : NeighborMooseVariableInterface<Real>(
49 : this, true, Moose::VarKindType::VAR_SOLVER, Moose::VarFieldType::VAR_FIELD_STANDARD),
50 101 : _secondary(_mesh.getBoundaryID(getParam<BoundaryName>("secondary"))),
51 101 : _primary(_mesh.getBoundaryID(getParam<BoundaryName>("primary"))),
52 101 : _var(_sys.getFieldVariable<Real>(_tid, parameters.get<NonlinearVariableName>("variable"))),
53 :
54 101 : _primary_q_point(_assembly.qPointsFace()),
55 101 : _primary_qrule(_assembly.qRuleFace()),
56 :
57 101 : _penetration_locator(
58 101 : getPenetrationLocator(getParam<BoundaryName>("primary"),
59 : getParam<BoundaryName>("secondary"),
60 202 : Utility::string_to_enum<Order>(getParam<MooseEnum>("order")))),
61 :
62 101 : _current_node(_var.node()),
63 101 : _current_primary(_var.neighbor()),
64 101 : _u_secondary(_var.dofValues()),
65 101 : _phi_secondary(1), // One entry
66 101 : _test_secondary(1), // One entry
67 :
68 101 : _primary_var(*getVar("primary_variable", 0)),
69 101 : _primary_var_num(_primary_var.number()),
70 :
71 101 : _phi_primary(_assembly.phiFaceNeighbor(_primary_var)),
72 101 : _grad_phi_primary(_assembly.gradPhiFaceNeighbor(_primary_var)),
73 :
74 101 : _test_primary(_var.phiFaceNeighbor()),
75 101 : _grad_test_primary(_var.gradPhiFaceNeighbor()),
76 :
77 101 : _u_primary(_primary_var.slnNeighbor()),
78 101 : _grad_u_primary(_primary_var.gradSlnNeighbor()),
79 :
80 101 : _dof_map(_sys.dofMap()),
81 101 : _node_to_elem_map(_mesh.nodeToElemMap()),
82 :
83 101 : _overwrite_secondary_residual(true),
84 303 : _primary_JxW(_assembly.JxWNeighbor())
85 : {
86 101 : addMooseVariableDependency(&_var);
87 :
88 101 : if (parameters.isParamValid("tangential_tolerance"))
89 : {
90 0 : _penetration_locator.setTangentialTolerance(getParam<Real>("tangential_tolerance"));
91 : }
92 101 : if (parameters.isParamValid("normal_smoothing_distance"))
93 : {
94 0 : _penetration_locator.setNormalSmoothingDistance(getParam<Real>("normal_smoothing_distance"));
95 : }
96 101 : if (parameters.isParamValid("normal_smoothing_method"))
97 : {
98 0 : _penetration_locator.setNormalSmoothingMethod(
99 0 : parameters.get<std::string>("normal_smoothing_method"));
100 : }
101 : // Put a "1" into test_secondary
102 : // will always only have one entry that is 1
103 101 : _test_secondary[0].push_back(1);
104 101 : }
105 :
106 202 : NodeFaceConstraint::~NodeFaceConstraint()
107 : {
108 101 : _phi_secondary.release();
109 101 : _test_secondary.release();
110 101 : }
111 :
112 : void
113 2128 : NodeFaceConstraint::computeSecondaryValue(NumericVector<Number> & current_solution)
114 : {
115 2128 : const dof_id_type & dof_idx = _var.nodalDofIndex();
116 2128 : _qp = 0;
117 2128 : current_solution.set(dof_idx, computeQpSecondaryValue());
118 2128 : }
119 :
120 : void
121 14004 : NodeFaceConstraint::residualSetup()
122 : {
123 14004 : _secondary_residual_computed = false;
124 14004 : }
125 :
126 : Real
127 44234 : NodeFaceConstraint::secondaryResidual() const
128 : {
129 : mooseAssert(_secondary_residual_computed,
130 : "The secondary residual has not yet been computed, so the value will be garbage!");
131 44234 : return _secondary_residual;
132 : }
133 :
134 : void
135 44234 : NodeFaceConstraint::computeResidual()
136 : {
137 44234 : _qp = 0;
138 :
139 44234 : prepareVectorTagNeighbor(_assembly, _primary_var.number());
140 221170 : for (_i = 0; _i < _test_primary.size(); _i++)
141 176936 : _local_re(_i) += computeQpResidual(Moose::Primary);
142 44234 : accumulateTaggedLocalResidual();
143 :
144 44234 : prepareVectorTag(_assembly, _var.number());
145 44234 : _i = 0;
146 44234 : _secondary_residual = _local_re(0) = computeQpResidual(Moose::Secondary);
147 44234 : assignTaggedLocalResidual();
148 44234 : _secondary_residual_computed = true;
149 44234 : }
150 :
151 : void
152 4184 : NodeFaceConstraint::computeJacobian()
153 : {
154 4184 : getConnectedDofIndices(_var.number());
155 :
156 4184 : _Kee.resize(_test_secondary.size(), _connected_dof_indices.size());
157 4184 : _Kne.resize(_test_primary.size(), _connected_dof_indices.size());
158 :
159 4184 : _phi_secondary.resize(_connected_dof_indices.size());
160 :
161 4184 : _qp = 0;
162 :
163 : // Fill up _phi_secondary so that it is 1 when j corresponds to this dof and 0 for every other dof
164 : // This corresponds to evaluating all of the connected shape functions at _this_ node
165 26716 : for (unsigned int j = 0; j < _connected_dof_indices.size(); j++)
166 : {
167 22532 : _phi_secondary[j].resize(1);
168 :
169 22532 : if (_connected_dof_indices[j] == _var.nodalDofIndex())
170 4184 : _phi_secondary[j][_qp] = 1.0;
171 : else
172 18348 : _phi_secondary[j][_qp] = 0.0;
173 : }
174 :
175 8368 : for (_i = 0; _i < _test_secondary.size(); _i++)
176 : // Loop over the connected dof indices so we can get all the jacobian contributions
177 26716 : for (_j = 0; _j < _connected_dof_indices.size(); _j++)
178 22532 : _Kee(_i, _j) += computeQpJacobian(Moose::SecondarySecondary);
179 :
180 : // Just do a direct assignment here because the Jacobian coming from assembly has already been
181 : // properly sized according to the neighbor _var dof indices. It has also been zeroed
182 4184 : prepareMatrixTagNeighbor(_assembly, _var.number(), _var.number(), Moose::ElementNeighbor, _Ken);
183 4184 : if (_Ken.m() && _Ken.n())
184 8224 : for (_i = 0; _i < _test_secondary.size(); _i++)
185 20560 : for (_j = 0; _j < _phi_primary.size(); _j++)
186 16448 : _Ken(_i, _j) += computeQpJacobian(Moose::SecondaryPrimary);
187 : // Don't accumulate here
188 :
189 20920 : for (_i = 0; _i < _test_primary.size(); _i++)
190 : // Loop over the connected dof indices so we can get all the jacobian contributions
191 106864 : for (_j = 0; _j < _connected_dof_indices.size(); _j++)
192 90128 : _Kne(_i, _j) += computeQpJacobian(Moose::PrimarySecondary);
193 :
194 4184 : prepareMatrixTagNeighbor(
195 4184 : _assembly, _primary_var.number(), _var.number(), Moose::NeighborNeighbor);
196 4184 : if (_local_ke.m() && _local_ke.n())
197 20560 : for (_i = 0; _i < _test_primary.size(); _i++)
198 82240 : for (_j = 0; _j < _phi_primary.size(); _j++)
199 65792 : _local_ke(_i, _j) += computeQpJacobian(Moose::PrimaryPrimary);
200 4184 : accumulateTaggedLocalMatrix();
201 4184 : }
202 :
203 : void
204 72 : NodeFaceConstraint::computeOffDiagJacobian(const unsigned int jvar_num)
205 : {
206 72 : getConnectedDofIndices(jvar_num);
207 :
208 72 : _Kee.resize(_test_secondary.size(), _connected_dof_indices.size());
209 72 : _Kne.resize(_test_primary.size(), _connected_dof_indices.size());
210 :
211 72 : _phi_secondary.resize(_connected_dof_indices.size());
212 :
213 72 : _qp = 0;
214 :
215 72 : auto primary_jsize = getVariable(jvar_num).dofIndicesNeighbor().size();
216 :
217 : // Fill up _phi_secondary so that it is 1 when j corresponds to this dof and 0 for every other dof
218 : // This corresponds to evaluating all of the connected shape functions at _this_ node
219 72 : for (unsigned int j = 0; j < _connected_dof_indices.size(); j++)
220 : {
221 0 : _phi_secondary[j].resize(1);
222 :
223 0 : if (_connected_dof_indices[j] == _var.nodalDofIndex())
224 0 : _phi_secondary[j][_qp] = 1.0;
225 : else
226 0 : _phi_secondary[j][_qp] = 0.0;
227 : }
228 :
229 144 : for (_i = 0; _i < _test_secondary.size(); _i++)
230 : // Loop over the connected dof indices so we can get all the jacobian contributions
231 72 : for (_j = 0; _j < _connected_dof_indices.size(); _j++)
232 0 : _Kee(_i, _j) += computeQpOffDiagJacobian(Moose::SecondarySecondary, jvar_num);
233 :
234 : // Just do a direct assignment here because the Jacobian coming from assembly has already been
235 : // properly sized according to the jvar neighbor dof indices. It has also been zeroed
236 72 : prepareMatrixTagNeighbor(_assembly, _var.number(), jvar_num, Moose::ElementNeighbor, _Ken);
237 144 : for (_i = 0; _i < _test_secondary.size(); _i++)
238 360 : for (_j = 0; _j < primary_jsize; _j++)
239 288 : _Ken(_i, _j) += computeQpOffDiagJacobian(Moose::SecondaryPrimary, jvar_num);
240 : // Don't accumulate here
241 :
242 72 : if (_Kne.m() && _Kne.n())
243 0 : for (_i = 0; _i < _test_primary.size(); _i++)
244 : // Loop over the connected dof indices so we can get all the jacobian contributions
245 0 : for (_j = 0; _j < _connected_dof_indices.size(); _j++)
246 0 : _Kne(_i, _j) += computeQpOffDiagJacobian(Moose::PrimarySecondary, jvar_num);
247 :
248 72 : prepareMatrixTagNeighbor(_assembly, _primary_var.number(), jvar_num, Moose::NeighborNeighbor);
249 360 : for (_i = 0; _i < _test_primary.size(); _i++)
250 1440 : for (_j = 0; _j < primary_jsize; _j++)
251 1152 : _local_ke(_i, _j) += computeQpOffDiagJacobian(Moose::PrimaryPrimary, jvar_num);
252 72 : accumulateTaggedLocalMatrix();
253 72 : }
254 :
255 : void
256 4256 : NodeFaceConstraint::getConnectedDofIndices(unsigned int var_num)
257 : {
258 4256 : MooseVariableFEBase & var = _sys.getVariable(0, var_num);
259 :
260 4256 : _connected_dof_indices.clear();
261 4256 : std::set<dof_id_type> unique_dof_indices;
262 :
263 4256 : auto node_to_elem_pair = _node_to_elem_map.find(_current_node->id());
264 : mooseAssert(node_to_elem_pair != _node_to_elem_map.end(), "Missing entry in node to elem map");
265 4256 : const std::vector<dof_id_type> & elems = node_to_elem_pair->second;
266 :
267 : // Get the dof indices from each elem connected to the node
268 11446 : for (const auto & cur_elem : elems)
269 : {
270 7190 : std::vector<dof_id_type> dof_indices;
271 :
272 7190 : var.getDofIndices(_mesh.elemPtr(cur_elem), dof_indices);
273 :
274 35518 : for (const auto & dof : dof_indices)
275 28328 : unique_dof_indices.insert(dof);
276 7190 : }
277 :
278 26788 : for (const auto & dof : unique_dof_indices)
279 22532 : _connected_dof_indices.push_back(dof);
280 4256 : }
281 :
282 : bool
283 52602 : NodeFaceConstraint::overwriteSecondaryResidual()
284 : {
285 52602 : return _overwrite_secondary_residual;
286 : }
287 :
288 : const std::set<BoundaryID> &
289 0 : NodeFaceConstraint::buildBoundaryIDs()
290 : {
291 0 : _boundary_ids.insert(_primary);
292 0 : _boundary_ids.insert(_secondary);
293 0 : return _boundary_ids;
294 : }
295 :
296 : std::set<SubdomainID>
297 13 : NodeFaceConstraint::getSecondaryConnectedBlocks() const
298 : {
299 13 : return _mesh.getBoundaryConnectedBlocks(_secondary);
300 : }
|
