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