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 : // App includes
11 : #include "AugmentSparsityOnInterface.h"
12 : #include "Executioner.h"
13 : #include "FEProblemBase.h"
14 : #include "MooseApp.h"
15 :
16 : // libMesh includes
17 : #include "libmesh/elem.h"
18 : #include "libmesh/mesh_base.h"
19 : #include "libmesh/boundary_info.h"
20 :
21 : registerMooseObject("MooseApp", AugmentSparsityOnInterface);
22 :
23 : using namespace libMesh;
24 :
25 : InputParameters
26 23537 : AugmentSparsityOnInterface::validParams()
27 : {
28 23537 : InputParameters params = RelationshipManager::validParams();
29 94148 : params.addRequiredParam<BoundaryName>("primary_boundary",
30 : "The name of the primary boundary sideset.");
31 94148 : params.addRequiredParam<BoundaryName>("secondary_boundary",
32 : "The name of the secondary boundary sideset.");
33 94148 : params.addRequiredParam<SubdomainName>("primary_subdomain",
34 : "The name of the primary lower dimensional subdomain.");
35 94148 : params.addRequiredParam<SubdomainName>("secondary_subdomain",
36 : "The name of the secondary lower dimensional subdomain.");
37 70611 : params.addParam<bool>(
38 : "ghost_point_neighbors",
39 47074 : false,
40 : "Whether we should ghost point neighbors of secondary lower-dimensional elements and "
41 : "also their interior parents and mortar interface couples for applications such as mortar "
42 : "nodal auxiliary kernels and nodal-normal contact Jacobians.");
43 70611 : params.addParam<bool>(
44 : "ghost_higher_d_neighbors",
45 47074 : false,
46 : "Whether we should ghost higher-dimensional neighbors. This is necessary when we are doing "
47 : "second order mortar with finite volume primal variables, because in order for the method to "
48 : "be second order we must use cell gradients, which couples in the neighbor cells.");
49 :
50 : // We want to wait until our mortar mesh has been built before trying to delete remote elements.
51 : // And our mortar mesh cannot be built until the entire mesh has been generated. By setting this
52 : // parameter to false we will make sure that any prepare_for_use calls during the mesh generation
53 : // phase will not delete remote elements *and* we will set a flag on the moose mesh saying that we
54 : // need to delete remote elements after the addition of late geometric ghosting functors
55 : // (including this ghosting functor)
56 23537 : params.set<bool>("attach_geometric_early") = false;
57 23537 : return params;
58 0 : }
59 :
60 8580 : AugmentSparsityOnInterface::AugmentSparsityOnInterface(const InputParameters & params)
61 : : RelationshipManager(params),
62 8580 : _primary_boundary_name(getParam<BoundaryName>("primary_boundary")),
63 17160 : _secondary_boundary_name(getParam<BoundaryName>("secondary_boundary")),
64 17160 : _primary_subdomain_name(getParam<SubdomainName>("primary_subdomain")),
65 17160 : _secondary_subdomain_name(getParam<SubdomainName>("secondary_subdomain")),
66 8580 : _is_coupling_functor(isType(Moose::RelationshipManagerType::COUPLING)),
67 17160 : _ghost_point_neighbors(getParam<bool>("ghost_point_neighbors")),
68 25740 : _ghost_higher_d_neighbors(getParam<bool>("ghost_higher_d_neighbors"))
69 : {
70 8580 : }
71 :
72 3288 : AugmentSparsityOnInterface::AugmentSparsityOnInterface(const AugmentSparsityOnInterface & other)
73 : : RelationshipManager(other),
74 3288 : _primary_boundary_name(other._primary_boundary_name),
75 3288 : _secondary_boundary_name(other._secondary_boundary_name),
76 3288 : _primary_subdomain_name(other._primary_subdomain_name),
77 3288 : _secondary_subdomain_name(other._secondary_subdomain_name),
78 3288 : _is_coupling_functor(other._is_coupling_functor),
79 3288 : _ghost_point_neighbors(other._ghost_point_neighbors),
80 3288 : _ghost_higher_d_neighbors(other._ghost_higher_d_neighbors)
81 : {
82 3288 : }
83 :
84 : void
85 3288 : AugmentSparsityOnInterface::internalInitWithMesh(const MeshBase &)
86 : {
87 3288 : }
88 :
89 : std::string
90 0 : AugmentSparsityOnInterface::getInfo() const
91 : {
92 0 : std::ostringstream oss;
93 0 : oss << "AugmentSparsityOnInterface";
94 0 : return oss.str();
95 0 : }
96 :
97 : void
98 490732 : AugmentSparsityOnInterface::ghostMortarInterfaceCouplings(
99 : const processor_id_type p,
100 : const Elem * const elem,
101 : map_type & coupled_elements,
102 : const AutomaticMortarGeneration & amg) const
103 : {
104 : // Look up elem in the mortar_interface_coupling data structure.
105 490732 : const auto & mic = amg.mortarInterfaceCoupling();
106 490732 : auto find_it = mic.find(elem->id());
107 490732 : if (find_it == mic.end())
108 396837 : return;
109 :
110 93895 : const auto & coupled_set = find_it->second;
111 :
112 596279 : for (const auto coupled_elem_id : coupled_set)
113 : {
114 : // Coupling data can precede a mesh rebuild, so IDs absent from the current mesh are ignored.
115 502384 : const Elem * const coupled_elem = _mesh->query_elem_ptr(coupled_elem_id);
116 502384 : if (coupled_elem && coupled_elem->processor_id() != p)
117 402817 : coupled_elements.emplace(coupled_elem, _null_mat);
118 : }
119 : }
120 :
121 : void
122 350 : AugmentSparsityOnInterface::ghostLowerDSecondaryElemPointNeighbors(
123 : const processor_id_type p,
124 : const Elem * const query_elem,
125 : map_type & coupled_elements,
126 : const SubdomainID secondary_subdomain_id,
127 : const AutomaticMortarGeneration & amg) const
128 : {
129 : // Node ownership can follow a higher-dimensional parent, so the query may be a secondary face,
130 : // its interior parent, or a paired primary parent. Start from every associated secondary face so
131 : // the ghosting graph contains the complete remote nodal-normal support in either row direction.
132 350 : std::vector<const Elem *> secondary_lower_elem_candidates;
133 350 : if (query_elem->subdomain_id() == secondary_subdomain_id)
134 105 : secondary_lower_elem_candidates.push_back(query_elem);
135 :
136 : // Coupling functors can run before mortar data is rebuilt after mesh refinement. Use the
137 : // relationship manager's current mesh for face topology instead of AMG-owned element pointers.
138 1400 : for (const auto side : query_elem->side_index_range())
139 1050 : if (const Elem * const lower_elem = _moose_mesh->getLowerDElem(query_elem, side);
140 1050 : lower_elem && lower_elem->subdomain_id() == secondary_subdomain_id)
141 105 : secondary_lower_elem_candidates.push_back(lower_elem);
142 :
143 350 : const auto & mic = amg.mortarInterfaceCoupling();
144 350 : if (const auto find_it = mic.find(query_elem->id()); find_it != mic.end())
145 1190 : for (const auto coupled_elem_id : find_it->second)
146 910 : if (const Elem * const coupled_elem = _mesh->query_elem_ptr(coupled_elem_id);
147 910 : coupled_elem && coupled_elem->subdomain_id() == secondary_subdomain_id)
148 245 : secondary_lower_elem_candidates.push_back(coupled_elem);
149 :
150 350 : std::set<dof_id_type> secondary_lower_elems_handled;
151 805 : for (const Elem * const candidate : secondary_lower_elem_candidates)
152 : {
153 455 : if (candidate->subactive())
154 0 : continue;
155 :
156 : // Resolve inactive ancestors to the active faces used by the current coupling graph.
157 455 : std::vector<const Elem *> active_secondary_lower_elems;
158 455 : candidate->active_family_tree(active_secondary_lower_elems);
159 910 : for (const Elem * const secondary_lower_elem : active_secondary_lower_elems)
160 : {
161 910 : if (secondary_lower_elem->subdomain_id() != secondary_subdomain_id ||
162 910 : !secondary_lower_elems_handled.insert(secondary_lower_elem->id()).second)
163 105 : continue;
164 :
165 350 : std::set<const Elem *> secondary_lower_elem_point_neighbors;
166 350 : secondary_lower_elem->find_point_neighbors(secondary_lower_elem_point_neighbors);
167 350 : secondary_lower_elem_point_neighbors.insert(secondary_lower_elem);
168 :
169 1190 : for (const Elem * const neigh : secondary_lower_elem_point_neighbors)
170 : {
171 840 : if (!neigh->active() || neigh->subdomain_id() != secondary_subdomain_id)
172 0 : continue;
173 :
174 : // Ghosting functors return only elements not owned by the requesting processor.
175 840 : if (neigh->processor_id() != p)
176 350 : coupled_elements.emplace(neigh, _null_mat);
177 :
178 : // Every point-neighbor face contributes to the smoothed nodal normal, even when it has no
179 : // mortar segment. Its interior parent therefore contributes displacement derivatives.
180 840 : const Elem * const interior_parent = neigh->interior_parent();
181 : mooseAssert(interior_parent,
182 : "A lower-dimensional secondary mortar element must have an interior parent.");
183 840 : if (interior_parent->processor_id() != p)
184 350 : coupled_elements.emplace(interior_parent, _null_mat);
185 :
186 840 : ghostMortarInterfaceCouplings(p, neigh, coupled_elements, amg);
187 : }
188 350 : }
189 455 : }
190 350 : }
191 :
192 : void
193 33516 : AugmentSparsityOnInterface::ghostHigherDNeighbors(const processor_id_type p,
194 : const Elem * const query_elem,
195 : map_type & coupled_elements,
196 : const BoundaryID secondary_boundary_id,
197 : const SubdomainID secondary_subdomain_id,
198 : const AutomaticMortarGeneration & amg) const
199 : {
200 : // The coupling is this for second order FV: secondary lower dimensional elem dofs will depend on
201 : // higher-d secondary elem neighbor primal dofs and higher-d primary elem neighbor primal dofs.
202 : // Higher dimensional primal dofs only depend on the secondary lower dimensional LM dofs for the
203 : // current FV gap heat transfer use cases
204 :
205 33516 : if (query_elem->subdomain_id() != secondary_subdomain_id)
206 33138 : return;
207 :
208 378 : const BoundaryInfo & binfo = _mesh->get_boundary_info();
209 378 : const auto which_side = query_elem->interior_parent()->which_side_am_i(query_elem);
210 378 : if (!binfo.has_boundary_id(query_elem->interior_parent(), which_side, secondary_boundary_id))
211 0 : return;
212 :
213 378 : const auto & mic = amg.mortarInterfaceCoupling();
214 378 : auto find_it = mic.find(query_elem->id());
215 378 : if (find_it == mic.end())
216 : // Perhaps no projection onto primary
217 0 : return;
218 :
219 378 : const auto & lower_d_coupled_set = find_it->second;
220 1512 : for (const auto coupled_elem_id : lower_d_coupled_set)
221 : {
222 1134 : auto * const coupled_elem = _mesh->elem_ptr(coupled_elem_id);
223 1134 : if (coupled_elem->dim() == query_elem->dim())
224 : // lower-d-elem
225 378 : continue;
226 :
227 756 : std::vector<const Elem *> active_neighbors;
228 :
229 3780 : for (auto s : coupled_elem->side_index_range())
230 : {
231 3024 : const Elem * const neigh = coupled_elem->neighbor_ptr(s);
232 :
233 3024 : if (!neigh || neigh == remote_elem)
234 828 : continue;
235 :
236 : // With any kind of neighbor, we need to couple to all the
237 : // active descendants on our side.
238 2196 : neigh->active_family_tree_by_neighbor(active_neighbors, coupled_elem);
239 :
240 4392 : for (const auto & neighbor : active_neighbors)
241 2196 : if (neighbor->processor_id() != p)
242 : {
243 : mooseAssert(
244 : neighbor->subdomain_id() != secondary_subdomain_id,
245 : "Ensure that we aren't missing potential erasures from the secondary-to-msms map");
246 2196 : coupled_elements.emplace(neighbor, _null_mat);
247 : }
248 : }
249 756 : }
250 : }
251 :
252 : void
253 255167 : AugmentSparsityOnInterface::operator()(const MeshBase::const_element_iterator & range_begin,
254 : const MeshBase::const_element_iterator & range_end,
255 : const processor_id_type p,
256 : map_type & coupled_elements)
257 : {
258 : // We ask the user to pass boundary names instead of ids to our constraint object. However, We
259 : // are unable to get the boundary ids from boundary names until we've attached the MeshBase object
260 : // to the MooseMesh
261 255167 : const bool generating_mesh = !_moose_mesh->getMeshPtr();
262 : const auto primary_boundary_id = generating_mesh
263 255167 : ? Moose::INVALID_BOUNDARY_ID
264 255167 : : _moose_mesh->getBoundaryID(_primary_boundary_name);
265 : const auto secondary_boundary_id = generating_mesh
266 255167 : ? Moose::INVALID_BOUNDARY_ID
267 255167 : : _moose_mesh->getBoundaryID(_secondary_boundary_name);
268 : const auto primary_subdomain_id = generating_mesh
269 255167 : ? Moose::INVALID_BLOCK_ID
270 255167 : : _moose_mesh->getSubdomainID(_primary_subdomain_name);
271 : const auto secondary_subdomain_id = generating_mesh
272 255167 : ? Moose::INVALID_BLOCK_ID
273 255167 : : _moose_mesh->getSubdomainID(_secondary_subdomain_name);
274 :
275 : const AutomaticMortarGeneration * const amg =
276 510334 : _app.getExecutioner() ? &_app.getExecutioner()->feProblem().getMortarInterface(
277 255167 : std::make_pair(primary_boundary_id, secondary_boundary_id),
278 0 : std::make_pair(primary_subdomain_id, secondary_subdomain_id),
279 255167 : _use_displaced_mesh)
280 255167 : : nullptr;
281 :
282 : // If we're on a dynamic mesh or we have not yet constructed the mortar mesh, we need to ghost the
283 : // entire interface because we don't know a priori what elements will project onto what. We *do
284 : // not* add the whole interface if we are a coupling functor because it is very expensive. This is
285 : // because when building the sparsity pattern, we call through to the ghosting functors with one
286 : // element at a time (and then below we do a loop over all the mesh's active elements). It's
287 : // perhaps faster in this case to deal with mallocs coming out of MatSetValues, especially if the
288 : // mesh displacements are relatively small
289 255167 : if ((!amg || _use_displaced_mesh) && !_is_coupling_functor)
290 : {
291 127048 : for (const Elem * const elem : _mesh->active_element_ptr_range())
292 : {
293 63088 : if (generating_mesh)
294 : {
295 : // We are still generating the mesh, so it's possible we don't even have the right boundary
296 : // ids created yet! So we actually ghost all boundary elements and all lower dimensional
297 : // elements who have parents on a boundary
298 0 : if (elem->on_boundary())
299 0 : coupled_elements.insert(std::make_pair(elem, _null_mat));
300 0 : else if (const Elem * const ip = elem->interior_parent())
301 : {
302 0 : if (ip->on_boundary())
303 0 : coupled_elements.insert(std::make_pair(elem, _null_mat));
304 : }
305 : }
306 : else
307 : {
308 : // We've finished generating our mesh so we can be selective and only ghost elements lying
309 : // in our lower-dimensional subdomains and their interior parents
310 :
311 : mooseAssert(primary_boundary_id != Moose::INVALID_BOUNDARY_ID,
312 : "Primary boundary id should exist by now.");
313 : mooseAssert(secondary_boundary_id != Moose::INVALID_BOUNDARY_ID,
314 : "Secondary boundary id should exist by now.");
315 : mooseAssert(primary_subdomain_id != Moose::INVALID_BLOCK_ID,
316 : "Primary subdomain id should exist by now.");
317 : mooseAssert(secondary_subdomain_id != Moose::INVALID_BLOCK_ID,
318 : "Secondary subdomain id should exist by now.");
319 :
320 : // Higher-dimensional boundary elements
321 63088 : const BoundaryInfo & binfo = _mesh->get_boundary_info();
322 :
323 258564 : for (auto side : elem->side_index_range())
324 292112 : if ((elem->processor_id() != p) &&
325 96636 : (binfo.has_boundary_id(elem, side, primary_boundary_id) ||
326 94656 : binfo.has_boundary_id(elem, side, secondary_boundary_id)))
327 4080 : coupled_elements.insert(std::make_pair(elem, _null_mat));
328 :
329 : // Lower dimensional subdomain elements
330 92560 : if ((elem->processor_id() != p) && (elem->subdomain_id() == primary_subdomain_id ||
331 29472 : elem->subdomain_id() == secondary_subdomain_id))
332 : {
333 3864 : coupled_elements.insert(std::make_pair(elem, _null_mat));
334 :
335 : #ifndef NDEBUG
336 : // let's do some safety checks
337 : const Elem * const ip = elem->interior_parent();
338 : mooseAssert(ip,
339 : "We should have set interior parents for all of our lower-dimensional mortar "
340 : "subdomains");
341 : auto side = ip->which_side_am_i(elem);
342 : auto bnd_id = elem->subdomain_id() == primary_subdomain_id ? primary_boundary_id
343 : : secondary_boundary_id;
344 : mooseAssert(_mesh->get_boundary_info().has_boundary_id(ip, side, bnd_id),
345 : "The interior parent for the lower-dimensional element does not lie on the "
346 : "boundary");
347 : #endif
348 : }
349 : }
350 872 : }
351 872 : }
352 : // For a static mesh (or for determining a sparsity pattern approximation on a displaced mesh) we
353 : // can just ghost the coupled elements determined during mortar mesh generation
354 254295 : else if (amg)
355 : {
356 1234079 : for (const Elem * const elem : as_range(range_begin, range_end))
357 : {
358 489892 : ghostMortarInterfaceCouplings(p, elem, coupled_elements, *amg);
359 :
360 489892 : if (_ghost_point_neighbors)
361 350 : ghostLowerDSecondaryElemPointNeighbors(
362 : p, elem, coupled_elements, secondary_subdomain_id, *amg);
363 489892 : if (_ghost_higher_d_neighbors)
364 33516 : ghostHigherDNeighbors(
365 : p, elem, coupled_elements, secondary_boundary_id, secondary_subdomain_id, *amg);
366 254295 : } // end for loop over input range
367 : } // end if amg
368 255167 : }
369 :
370 : bool
371 43953 : AugmentSparsityOnInterface::operator>=(const RelationshipManager & other) const
372 : {
373 43953 : if (auto asoi = dynamic_cast<const AugmentSparsityOnInterface *>(&other))
374 : {
375 37006 : if (_primary_boundary_name == asoi->_primary_boundary_name &&
376 29380 : _secondary_boundary_name == asoi->_secondary_boundary_name &&
377 29380 : _primary_subdomain_name == asoi->_primary_subdomain_name &&
378 14690 : _secondary_subdomain_name == asoi->_secondary_subdomain_name &&
379 14690 : _ghost_point_neighbors >= asoi->_ghost_point_neighbors &&
380 37006 : _ghost_higher_d_neighbors >= asoi->_ghost_higher_d_neighbors && baseGreaterEqual(*asoi))
381 6446 : return true;
382 : }
383 37507 : return false;
384 : }
385 :
386 : std::unique_ptr<GhostingFunctor>
387 3288 : AugmentSparsityOnInterface::clone() const
388 : {
389 3288 : return _app.getFactory().copyConstruct(*this);
390 : }
|
