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 "ElementDeletionGeneratorBase.h"
11 : #include "libmesh/remote_elem.h"
12 :
13 : #include "MooseMeshUtils.h"
14 : #include "CastUniquePointer.h"
15 :
16 : InputParameters
17 12877 : ElementDeletionGeneratorBase::validParams()
18 : {
19 12877 : InputParameters params = MeshGenerator::validParams();
20 :
21 51508 : params.addRequiredParam<MeshGeneratorName>("input", "The mesh we want to modify");
22 38631 : params.addParam<bool>("delete_exteriors",
23 25754 : true,
24 : "Whether to delete lower-d elements whose interior parents are deleted");
25 38631 : params.addParam<BoundaryName>("new_boundary",
26 : "optional boundary name to assign to the cut surface");
27 :
28 12877 : return params;
29 0 : }
30 :
31 315 : ElementDeletionGeneratorBase::ElementDeletionGeneratorBase(const InputParameters & parameters)
32 : : MeshGenerator(parameters),
33 315 : _input(getMesh("input")),
34 630 : _assign_boundary(isParamValid("new_boundary")),
35 630 : _delete_exteriors(getParam<bool>("delete_exteriors")),
36 684 : _boundary_name(_assign_boundary ? getParam<BoundaryName>("new_boundary") : "")
37 : {
38 315 : }
39 :
40 : std::unique_ptr<MeshBase>
41 296 : ElementDeletionGeneratorBase::generate()
42 : {
43 296 : std::unique_ptr<MeshBase> mesh = std::move(_input);
44 :
45 : // Make sure that the mesh is prepared. We'll just do a full
46 : // prepare_for_use() here, since we expect to be able to use
47 : // neighbor pointers, multiple caches, *and* on distributed meshes a
48 : // load-balanced partitioning.
49 296 : if (!mesh->is_prepared())
50 236 : mesh->prepare_for_use();
51 :
52 : // Elements that the deleter will remove
53 296 : std::set<Elem *> deleteable_elems;
54 :
55 : // First let's figure out which elements need to be deleted
56 191851 : for (auto & elem : mesh->element_ptr_range())
57 : {
58 191555 : if (shouldDelete(elem))
59 153823 : deleteable_elems.insert(elem);
60 296 : }
61 :
62 : // check for dangling interior parents
63 191851 : for (auto & elem : mesh->element_ptr_range())
64 191555 : if (elem->dim() < mesh->mesh_dimension() && deleteable_elems.count(elem->interior_parent()) > 0)
65 : {
66 330 : if (_delete_exteriors)
67 38 : deleteable_elems.insert(elem);
68 : else
69 292 : elem->set_interior_parent(nullptr);
70 296 : }
71 :
72 : /**
73 : * If we are in parallel we'd better have a consistent idea of what
74 : * should be deleted. This can't be checked cheaply.
75 : */
76 : #ifdef DEBUG
77 : dof_id_type pmax_elem_id = mesh->max_elem_id();
78 : mesh->comm().max(pmax_elem_id);
79 :
80 : for (dof_id_type i = 0; i != pmax_elem_id; ++i)
81 : {
82 : Elem * elem = mesh->query_elem_ptr(i);
83 : bool is_deleteable = elem && deleteable_elems.count(elem);
84 :
85 : libmesh_assert(mesh->comm().semiverify(elem ? &is_deleteable : libmesh_nullptr));
86 : }
87 : #endif
88 :
89 : // Get the BoundaryID from the mesh
90 296 : boundary_id_type boundary_id = 0;
91 296 : if (_assign_boundary)
92 156 : boundary_id = MooseMeshUtils::getBoundaryIDs(*mesh, {_boundary_name}, true)[0];
93 :
94 : // Get a reference to our BoundaryInfo object for later use
95 296 : BoundaryInfo & boundary_info = mesh->get_boundary_info();
96 :
97 : /**
98 : * Delete all of the elements
99 : *
100 : * TODO: We need to sort these not because they have to be deleted in a certain order in libMesh,
101 : * but because the order of deletion might impact what happens to any existing sidesets or
102 : * nodesets.
103 : */
104 154157 : for (auto & elem : deleteable_elems)
105 : {
106 : // On distributed meshes, we'll need neighbor links to be useable
107 : // shortly, so we can't just leave dangling pointers.
108 : //
109 : // FIXME - this could be made AMR-aware and refactored into
110 : // libMesh - roystgnr
111 153861 : unsigned int n_sides = elem->n_sides();
112 790170 : for (unsigned int n = 0; n != n_sides; ++n)
113 : {
114 636309 : Elem * neighbor = elem->neighbor_ptr(n);
115 636309 : if (!neighbor || neighbor == remote_elem)
116 322468 : continue;
117 :
118 313841 : const unsigned int return_side = neighbor->which_neighbor_am_i(elem);
119 :
120 313841 : if (neighbor->neighbor_ptr(return_side) == elem)
121 : {
122 313841 : neighbor->set_neighbor(return_side, nullptr);
123 :
124 : // assign cut surface boundary
125 313841 : if (_assign_boundary)
126 20996 : boundary_info.add_side(neighbor, return_side, boundary_id);
127 : }
128 : }
129 :
130 153861 : mesh->delete_elem(elem);
131 : }
132 :
133 : /**
134 : * If we are on a distributed mesh, we may have deleted elements
135 : * which are remote_elem links on other processors, and we need
136 : * to make neighbor and interior parent links into NULL pointers
137 : * (i.e. domain boundaries in the former case) instead.
138 : */
139 296 : if (!mesh->is_serial())
140 : {
141 46 : const processor_id_type my_n_proc = mesh->n_processors();
142 46 : const processor_id_type my_proc_id = mesh->processor_id();
143 : typedef std::vector<std::pair<dof_id_type, unsigned int>> vec_type;
144 46 : std::vector<vec_type> queries(my_n_proc);
145 :
146 : // Loop over the elements looking for those with remote neighbors
147 : // or interior parents.
148 : // The ghost_elements iterators in libMesh need to be updated
149 : // before we can use them safely here, so we'll test for
150 : // ghost-vs-local manually.
151 3039 : for (const auto & elem : mesh->element_ptr_range())
152 : {
153 2993 : const processor_id_type pid = elem->processor_id();
154 2993 : if (pid == my_proc_id)
155 2583 : continue;
156 :
157 410 : const unsigned int n_sides = elem->n_sides();
158 2640 : for (unsigned int n = 0; n != n_sides; ++n)
159 2230 : if (elem->neighbor_ptr(n) == remote_elem)
160 405 : queries[pid].push_back(std::make_pair(elem->id(), n));
161 :
162 : // Use an OOB side index to encode "interior_parent". We will use this OOB index later
163 410 : if (elem->interior_parent() == remote_elem)
164 0 : queries[pid].push_back(std::make_pair(elem->id(), n_sides));
165 46 : }
166 :
167 46 : const auto queries_tag = mesh->comm().get_unique_tag(),
168 46 : replies_tag = mesh->comm().get_unique_tag();
169 :
170 138 : std::vector<Parallel::Request> query_requests(my_n_proc - 1), reply_requests(my_n_proc - 1);
171 :
172 : // Make all requests
173 138 : for (processor_id_type p = 0; p != my_n_proc; ++p)
174 : {
175 92 : if (p == my_proc_id)
176 46 : continue;
177 :
178 46 : Parallel::Request & request = query_requests[p - (p > my_proc_id)];
179 :
180 46 : mesh->comm().send(p, queries[p], request, queries_tag);
181 : }
182 :
183 : // Reply to all requests
184 46 : std::vector<vec_type> responses(my_n_proc - 1);
185 :
186 92 : for (processor_id_type p = 1; p != my_n_proc; ++p)
187 : {
188 46 : vec_type query;
189 :
190 46 : Parallel::Status status(mesh->comm().probe(Parallel::any_source, queries_tag));
191 46 : const processor_id_type source_pid = cast_int<processor_id_type>(status.source());
192 :
193 46 : mesh->comm().receive(source_pid, query, queries_tag);
194 :
195 46 : Parallel::Request & request = reply_requests[p - 1];
196 :
197 451 : for (const auto & q : query)
198 : {
199 405 : const Elem * elem = mesh->elem_ptr(q.first);
200 405 : const unsigned int side = q.second;
201 : const Elem * target =
202 405 : (side >= elem->n_sides()) ? elem->interior_parent() : elem->neighbor_ptr(side);
203 :
204 405 : if (target == nullptr) // linked element was deleted!
205 23 : responses[p - 1].push_back(std::make_pair(elem->id(), side));
206 : }
207 :
208 46 : mesh->comm().send(source_pid, responses[p - 1], request, replies_tag);
209 46 : }
210 :
211 : // Process all incoming replies
212 92 : for (processor_id_type p = 1; p != my_n_proc; ++p)
213 : {
214 46 : Parallel::Status status(this->comm().probe(Parallel::any_source, replies_tag));
215 46 : const processor_id_type source_pid = cast_int<processor_id_type>(status.source());
216 :
217 46 : vec_type response;
218 :
219 46 : this->comm().receive(source_pid, response, replies_tag);
220 :
221 69 : for (const auto & r : response)
222 : {
223 23 : Elem * elem = mesh->elem_ptr(r.first);
224 23 : const unsigned int side = r.second;
225 :
226 23 : if (side < elem->n_sides())
227 : {
228 : mooseAssert(elem->neighbor_ptr(side) == remote_elem, "element neighbor != remote_elem");
229 :
230 23 : elem->set_neighbor(side, nullptr);
231 :
232 : // assign cut surface boundary
233 23 : if (_assign_boundary)
234 0 : boundary_info.add_side(elem, side, boundary_id);
235 : }
236 : else
237 : {
238 : mooseAssert(side == elem->n_sides(), "internal communication error");
239 : mooseAssert(elem->interior_parent() == remote_elem, "interior parent != remote_elem");
240 :
241 0 : elem->set_interior_parent(nullptr);
242 : }
243 : }
244 46 : }
245 :
246 46 : Parallel::wait(query_requests);
247 46 : Parallel::wait(reply_requests);
248 46 : }
249 :
250 296 : if (_assign_boundary)
251 : {
252 52 : boundary_info.sideset_name(boundary_id) = _boundary_name;
253 52 : boundary_info.nodeset_name(boundary_id) = _boundary_name;
254 : }
255 :
256 : /**
257 : * If we are on a ReplicatedMesh, deleting nodes and elements leaves
258 : * NULLs in the mesh datastructure. We ought to get rid of those.
259 : * For now, we'll call contract and notify the SetupMeshComplete
260 : * Action that we need to re-prepare the mesh.
261 : */
262 296 : mesh->contract();
263 296 : mesh->prepare_for_use();
264 :
265 592 : return dynamic_pointer_cast<MeshBase>(mesh);
266 348 : }
|
