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 42983 : NodeFaceConstraint::validParams()
24 : {
25 42983 : MooseEnum orders("FIRST SECOND THIRD FOURTH", "FIRST");
26 42983 : InputParameters params = Constraint::validParams();
27 42983 : params.addParam<BoundaryName>("secondary", "The boundary ID associated with the secondary side");
28 42983 : params.addParam<BoundaryName>("primary", "The boundary ID associated with the primary side");
29 42983 : params.addParam<Real>("tangential_tolerance",
30 : "Tangential distance to extend edges of contact surfaces");
31 42983 : params.addParam<Real>(
32 : "normal_smoothing_distance",
33 : "Distance from edge in parametric coordinates over which to smooth contact normal");
34 42983 : params.addParam<std::string>("normal_smoothing_method",
35 : "Method to use to smooth normals (edge_based|nodal_normal_based)");
36 42983 : params.addParam<MooseEnum>("order", orders, "The finite element order used for projections");
37 :
38 42983 : params.addCoupledVar("primary_variable", "The variable on the primary side of the domain");
39 :
40 85966 : return params;
41 42983 : }
42 :
43 93 : 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 93 : _secondary(_mesh.getBoundaryID(getParam<BoundaryName>("secondary"))),
51 93 : _primary(_mesh.getBoundaryID(getParam<BoundaryName>("primary"))),
52 93 : _var(_sys.getFieldVariable<Real>(_tid, parameters.get<NonlinearVariableName>("variable"))),
53 :
54 93 : _primary_q_point(_assembly.qPointsFace()),
55 93 : _primary_qrule(_assembly.qRuleFace()),
56 :
57 93 : _penetration_locator(
58 93 : getPenetrationLocator(getParam<BoundaryName>("primary"),
59 : getParam<BoundaryName>("secondary"),
60 186 : Utility::string_to_enum<Order>(getParam<MooseEnum>("order")))),
61 :
62 93 : _current_node(_var.node()),
63 93 : _current_primary(_var.neighbor()),
64 93 : _u_secondary(_var.dofValues()),
65 93 : _phi_secondary(1), // One entry
66 93 : _test_secondary(1), // One entry
67 :
68 93 : _primary_var(*getVar("primary_variable", 0)),
69 93 : _primary_var_num(_primary_var.number()),
70 :
71 93 : _phi_primary(_assembly.phiFaceNeighbor(_primary_var)),
72 93 : _grad_phi_primary(_assembly.gradPhiFaceNeighbor(_primary_var)),
73 :
74 93 : _test_primary(_var.phiFaceNeighbor()),
75 93 : _grad_test_primary(_var.gradPhiFaceNeighbor()),
76 :
77 93 : _u_primary(_primary_var.slnNeighbor()),
78 93 : _grad_u_primary(_primary_var.gradSlnNeighbor()),
79 :
80 93 : _dof_map(_sys.dofMap()),
81 93 : _node_to_elem_map(_mesh.nodeToElemMap()),
82 :
83 93 : _overwrite_secondary_residual(true),
84 279 : _primary_JxW(_assembly.JxWNeighbor())
85 : {
86 93 : addMooseVariableDependency(&_var);
87 :
88 93 : if (parameters.isParamValid("tangential_tolerance"))
89 : {
90 0 : _penetration_locator.setTangentialTolerance(getParam<Real>("tangential_tolerance"));
91 : }
92 93 : if (parameters.isParamValid("normal_smoothing_distance"))
93 : {
94 0 : _penetration_locator.setNormalSmoothingDistance(getParam<Real>("normal_smoothing_distance"));
95 : }
96 93 : 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 93 : _test_secondary[0].push_back(1);
104 93 : }
105 :
106 186 : NodeFaceConstraint::~NodeFaceConstraint()
107 : {
108 93 : _phi_secondary.release();
109 93 : _test_secondary.release();
110 93 : }
111 :
112 : void
113 1880 : NodeFaceConstraint::computeSecondaryValue(NumericVector<Number> & current_solution)
114 : {
115 1880 : const dof_id_type & dof_idx = _var.nodalDofIndex();
116 1880 : _qp = 0;
117 1880 : current_solution.set(dof_idx, computeQpSecondaryValue());
118 1880 : }
119 :
120 : void
121 12822 : NodeFaceConstraint::residualSetup()
122 : {
123 12822 : _secondary_residual_computed = false;
124 12822 : }
125 :
126 : Real
127 39276 : 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 39276 : return _secondary_residual;
132 : }
133 :
134 : void
135 39276 : NodeFaceConstraint::computeResidual()
136 : {
137 39276 : _qp = 0;
138 :
139 39276 : prepareVectorTagNeighbor(_assembly, _primary_var.number());
140 196380 : for (_i = 0; _i < _test_primary.size(); _i++)
141 157104 : _local_re(_i) += computeQpResidual(Moose::Primary);
142 39276 : accumulateTaggedLocalResidual();
143 :
144 39276 : prepareVectorTag(_assembly, _var.number());
145 39276 : _i = 0;
146 39276 : _secondary_residual = _local_re(0) = computeQpResidual(Moose::Secondary);
147 39276 : assignTaggedLocalResidual();
148 39276 : _secondary_residual_computed = true;
149 39276 : }
150 :
151 : void
152 3696 : NodeFaceConstraint::computeJacobian()
153 : {
154 3696 : getConnectedDofIndices(_var.number());
155 :
156 3696 : _Kee.resize(_test_secondary.size(), _connected_dof_indices.size());
157 3696 : _Kne.resize(_test_primary.size(), _connected_dof_indices.size());
158 :
159 3696 : _phi_secondary.resize(_connected_dof_indices.size());
160 :
161 3696 : _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 23600 : for (unsigned int j = 0; j < _connected_dof_indices.size(); j++)
166 : {
167 19904 : _phi_secondary[j].resize(1);
168 :
169 19904 : if (_connected_dof_indices[j] == _var.nodalDofIndex())
170 3696 : _phi_secondary[j][_qp] = 1.0;
171 : else
172 16208 : _phi_secondary[j][_qp] = 0.0;
173 : }
174 :
175 7392 : for (_i = 0; _i < _test_secondary.size(); _i++)
176 : // Loop over the connected dof indices so we can get all the jacobian contributions
177 23600 : for (_j = 0; _j < _connected_dof_indices.size(); _j++)
178 19904 : _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 3696 : prepareMatrixTagNeighbor(_assembly, _var.number(), _var.number(), Moose::ElementNeighbor, _Ken);
183 3696 : if (_Ken.m() && _Ken.n())
184 7264 : for (_i = 0; _i < _test_secondary.size(); _i++)
185 18160 : for (_j = 0; _j < _phi_primary.size(); _j++)
186 14528 : _Ken(_i, _j) += computeQpJacobian(Moose::SecondaryPrimary);
187 : // Don't accumulate here
188 :
189 18480 : for (_i = 0; _i < _test_primary.size(); _i++)
190 : // Loop over the connected dof indices so we can get all the jacobian contributions
191 94400 : for (_j = 0; _j < _connected_dof_indices.size(); _j++)
192 79616 : _Kne(_i, _j) += computeQpJacobian(Moose::PrimarySecondary);
193 :
194 3696 : prepareMatrixTagNeighbor(
195 3696 : _assembly, _primary_var.number(), _var.number(), Moose::NeighborNeighbor);
196 3696 : if (_local_ke.m() && _local_ke.n())
197 18160 : for (_i = 0; _i < _test_primary.size(); _i++)
198 72640 : for (_j = 0; _j < _phi_primary.size(); _j++)
199 58112 : _local_ke(_i, _j) += computeQpJacobian(Moose::PrimaryPrimary);
200 3696 : accumulateTaggedLocalMatrix();
201 3696 : }
202 :
203 : void
204 64 : NodeFaceConstraint::computeOffDiagJacobian(const unsigned int jvar_num)
205 : {
206 64 : getConnectedDofIndices(jvar_num);
207 :
208 64 : _Kee.resize(_test_secondary.size(), _connected_dof_indices.size());
209 64 : _Kne.resize(_test_primary.size(), _connected_dof_indices.size());
210 :
211 64 : _phi_secondary.resize(_connected_dof_indices.size());
212 :
213 64 : _qp = 0;
214 :
215 64 : 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 64 : 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 128 : for (_i = 0; _i < _test_secondary.size(); _i++)
230 : // Loop over the connected dof indices so we can get all the jacobian contributions
231 64 : 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 64 : prepareMatrixTagNeighbor(_assembly, _var.number(), jvar_num, Moose::ElementNeighbor, _Ken);
237 128 : for (_i = 0; _i < _test_secondary.size(); _i++)
238 320 : for (_j = 0; _j < primary_jsize; _j++)
239 256 : _Ken(_i, _j) += computeQpOffDiagJacobian(Moose::SecondaryPrimary, jvar_num);
240 : // Don't accumulate here
241 :
242 64 : 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 64 : prepareMatrixTagNeighbor(_assembly, _primary_var.number(), jvar_num, Moose::NeighborNeighbor);
249 320 : for (_i = 0; _i < _test_primary.size(); _i++)
250 1280 : for (_j = 0; _j < primary_jsize; _j++)
251 1024 : _local_ke(_i, _j) += computeQpOffDiagJacobian(Moose::PrimaryPrimary, jvar_num);
252 64 : accumulateTaggedLocalMatrix();
253 64 : }
254 :
255 : void
256 3760 : NodeFaceConstraint::getConnectedDofIndices(unsigned int var_num)
257 : {
258 3760 : MooseVariableFEBase & var = _sys.getVariable(0, var_num);
259 :
260 3760 : _connected_dof_indices.clear();
261 3760 : std::set<dof_id_type> unique_dof_indices;
262 :
263 3760 : 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 3760 : 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 10112 : for (const auto & cur_elem : elems)
269 : {
270 6352 : std::vector<dof_id_type> dof_indices;
271 :
272 6352 : var.getDofIndices(_mesh.elemPtr(cur_elem), dof_indices);
273 :
274 31376 : for (const auto & dof : dof_indices)
275 25024 : unique_dof_indices.insert(dof);
276 6352 : }
277 :
278 23664 : for (const auto & dof : unique_dof_indices)
279 19904 : _connected_dof_indices.push_back(dof);
280 3760 : }
281 :
282 : bool
283 46668 : NodeFaceConstraint::overwriteSecondaryResidual()
284 : {
285 46668 : 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 12 : NodeFaceConstraint::getSecondaryConnectedBlocks() const
298 : {
299 12 : return _mesh.getBoundaryConnectedBlocks(_secondary);
300 : }
|
