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