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 : #pragma once
11 :
12 : #include "ParallelUniqueId.h"
13 : #include "MooseMesh.h"
14 : #include "MooseTypes.h"
15 : #include "MooseException.h"
16 : #include "libmesh/libmesh_exceptions.h"
17 : #include "libmesh/elem.h"
18 :
19 : /**
20 : * Base class for assembly-like calculations.
21 : */
22 : template <typename RangeType>
23 : class ThreadedElementLoopBase
24 : {
25 : public:
26 : ThreadedElementLoopBase(MooseMesh & mesh);
27 :
28 : ThreadedElementLoopBase(ThreadedElementLoopBase & x, Threads::split split);
29 :
30 : virtual ~ThreadedElementLoopBase();
31 :
32 : virtual void operator()(const RangeType & range, bool bypass_threading = false);
33 :
34 : /**
35 : * Called before the element range loop
36 : */
37 : virtual void pre();
38 :
39 : /**
40 : * Called after the element range loop
41 : */
42 : virtual void post();
43 :
44 : /**
45 : * Assembly of the element (not including surface assembly)
46 : *
47 : * @param elem - active element
48 : */
49 : virtual void onElement(const Elem * elem);
50 :
51 : /**
52 : * Called before the element assembly
53 : *
54 : * @param elem - active element
55 : */
56 : virtual void preElement(const Elem * elem);
57 :
58 : /**
59 : * Called after the element assembly is done (including surface assembling)
60 : *
61 : * @param elem - active element
62 : */
63 : virtual void postElement(const Elem * elem);
64 :
65 : /**
66 : * Called before the boundary assembly
67 : *
68 : * @param elem - The element we are checking is on the boundary.
69 : * @param side - The side of the element in question.
70 : * @param bnd_id - ID of the boundary we are at
71 : * @param lower_d_elem - Lower dimensional element (e.g. Mortar)
72 : */
73 : virtual void preBoundary(const Elem * elem,
74 : unsigned int side,
75 : BoundaryID bnd_id,
76 : const Elem * lower_d_elem = nullptr);
77 :
78 : /**
79 : * Called when doing boundary assembling
80 : *
81 : * @param elem - The element we are checking is on the boundary.
82 : * @param side - The side of the element in question.
83 : * @param bnd_id - ID of the boundary we are at
84 : * @param lower_d_elem - Lower dimensional element (e.g. Mortar)
85 : */
86 : virtual void onBoundary(const Elem * elem,
87 : unsigned int side,
88 : BoundaryID bnd_id,
89 : const Elem * lower_d_elem = nullptr);
90 :
91 : /**
92 : * Called before evaluations on an element internal side
93 : *
94 : * @param elem - Element we are on
95 : * @param side - local side number of the element 'elem'
96 : */
97 : virtual void preInternalSide(const Elem * elem, unsigned int side);
98 :
99 : /**
100 : * Called after evaluations on an element internal side
101 : *
102 : * @param elem - Element we are on
103 : * @param side - local side number of the element 'elem'
104 : */
105 : virtual void postInternalSide(const Elem * elem, unsigned int side);
106 :
107 : /**
108 : * Called when doing internal edge assembling
109 : *
110 : * @param elem - Element we are on
111 : * @param side - local side number of the element 'elem'
112 : */
113 : virtual void onInternalSide(const Elem * elem, unsigned int side);
114 :
115 : /**
116 : * Called when iterating over external sides (no side neighbor)
117 : *
118 : * @param elem - Element we are on
119 : * @param side - local side number of the element 'elem'
120 : */
121 : virtual void onExternalSide(const Elem * elem, unsigned int side);
122 :
123 : /**
124 : * Called when doing interface assembling
125 : *
126 : * @param elem - Element we are on
127 : * @param side - local side number of the element 'elem'
128 : * @param bnd_id - ID of the interface we are at
129 : */
130 : virtual void onInterface(const Elem * elem, unsigned int side, BoundaryID bnd_id);
131 :
132 : /**
133 : * Called every time the current subdomain changes (i.e. the subdomain of _this_ element
134 : * is not the same as the subdomain of the last element). Beware of over-using this!
135 : * You might think that you can do some expensive stuff in here and get away with it...
136 : * but there are applications that have TONS of subdomains....
137 : */
138 : virtual void subdomainChanged();
139 :
140 : /**
141 : * Called every time the neighbor subdomain changes (i.e. the subdomain of _this_ neighbor
142 : * is not the same as the subdomain of the last neighbor). Beware of over-using this!
143 : * You might think that you can do some expensive stuff in here and get away with it...
144 : * but there are applications that have TONS of subdomains....
145 : */
146 : virtual void neighborSubdomainChanged();
147 :
148 : /**
149 : * Called if a MooseException is caught anywhere during the computation.
150 : * The single input parameter taken is a MooseException object.
151 : */
152 0 : virtual void caughtMooseException(MooseException &){};
153 :
154 : /**
155 : * Whether or not the loop should continue.
156 : *
157 : * @return true to keep going, false to stop.
158 : */
159 130106 : virtual bool keepGoing() { return true; }
160 :
161 : protected:
162 : MooseMesh & _mesh;
163 : THREAD_ID _tid;
164 :
165 : /// The subdomain for the current element
166 : SubdomainID _subdomain;
167 :
168 : /// The subdomain for the last element
169 : SubdomainID _old_subdomain;
170 :
171 : /// The subdomain for the current neighbor
172 : SubdomainID _neighbor_subdomain;
173 :
174 : /// The subdomain for the last neighbor
175 : SubdomainID _old_neighbor_subdomain;
176 :
177 : /// Print information about the loop ordering
178 82655 : virtual void printGeneralExecutionInformation() const {}
179 :
180 : /// Print information about the particular ordering of objects on each block
181 281093 : virtual void printBlockExecutionInformation() const {}
182 :
183 : /// Print information about the particular ordering of objects on each boundary
184 11949927 : virtual void printBoundaryExecutionInformation(const unsigned int /*bid*/) const {}
185 :
186 : /// Keep track of which blocks were visited
187 : mutable std::set<SubdomainID> _blocks_exec_printed;
188 :
189 : /// Keep track of which boundaries were visited
190 : mutable std::set<BoundaryID> _boundaries_exec_printed;
191 :
192 : /// Resets the set of blocks and boundaries visited
193 : void resetExecPrintedSets() const;
194 :
195 : /**
196 : * Whether to compute the internal side for the provided element-neighbor pair. Typically this
197 : * will return true if the element id is less than the neighbor id when the elements are equal
198 : * level, or when the element is more refined than the neighbor, and then false otherwise. One
199 : * type of loop where the logic will be different is when projecting stateful material properties
200 : */
201 : virtual bool shouldComputeInternalSide(const Elem & elem, const Elem & neighbor) const;
202 : };
203 :
204 : template <typename RangeType>
205 3923358 : ThreadedElementLoopBase<RangeType>::ThreadedElementLoopBase(MooseMesh & mesh) : _mesh(mesh)
206 : {
207 3923358 : }
208 :
209 : template <typename RangeType>
210 61 : ThreadedElementLoopBase<RangeType>::ThreadedElementLoopBase(ThreadedElementLoopBase & x,
211 : Threads::split /*split*/)
212 61 : : _mesh(x._mesh)
213 : {
214 61 : }
215 :
216 : template <typename RangeType>
217 4265769 : ThreadedElementLoopBase<RangeType>::~ThreadedElementLoopBase()
218 : {
219 4265769 : }
220 :
221 : template <typename RangeType>
222 : void
223 4206318 : ThreadedElementLoopBase<RangeType>::operator()(const RangeType & range, bool bypass_threading)
224 : {
225 : try
226 : {
227 : try
228 : {
229 4206318 : ParallelUniqueId puid;
230 4206318 : _tid = bypass_threading ? 0 : puid.id;
231 :
232 4206318 : pre();
233 4206318 : printGeneralExecutionInformation();
234 :
235 4206318 : _subdomain = Moose::INVALID_BLOCK_ID;
236 4206318 : _neighbor_subdomain = Moose::INVALID_BLOCK_ID;
237 4206318 : typename RangeType::const_iterator el = range.begin();
238 392823626 : for (el = range.begin(); el != range.end(); ++el)
239 : {
240 388617543 : if (!keepGoing())
241 5 : break;
242 :
243 388617538 : const Elem * elem = *el;
244 :
245 388617538 : preElement(elem);
246 :
247 388617538 : _old_subdomain = _subdomain;
248 388617538 : _subdomain = elem->subdomain_id();
249 388617538 : if (_subdomain != _old_subdomain)
250 : {
251 5296735 : subdomainChanged();
252 5296735 : printBlockExecutionInformation();
253 : }
254 :
255 388617538 : onElement(elem);
256 :
257 777218520 : if (_mesh.interiorLowerDBlocks().count(elem->subdomain_id()) > 0 ||
258 777218520 : _mesh.boundaryLowerDBlocks().count(elem->subdomain_id()) > 0)
259 : {
260 33640 : postElement(elem);
261 33640 : continue;
262 : }
263 :
264 1984866154 : for (unsigned int side = 0; side < elem->n_sides(); side++)
265 : {
266 1596282486 : std::vector<BoundaryID> boundary_ids = _mesh.getBoundaryIDs(elem, side);
267 1596282486 : const Elem * lower_d_elem = _mesh.getLowerDElem(elem, side);
268 :
269 1596282486 : if (boundary_ids.size() > 0)
270 116013272 : for (std::vector<BoundaryID>::iterator it = boundary_ids.begin();
271 232365470 : it != boundary_ids.end();
272 116352198 : ++it)
273 : {
274 116352205 : preBoundary(elem, side, *it, lower_d_elem);
275 116352205 : printBoundaryExecutionInformation(*it);
276 116352205 : onBoundary(elem, side, *it, lower_d_elem);
277 : }
278 :
279 1596282479 : const Elem * neighbor = elem->neighbor_ptr(side);
280 1596282479 : if (neighbor)
281 : {
282 1478472052 : preInternalSide(elem, side);
283 :
284 1478472052 : _old_neighbor_subdomain = _neighbor_subdomain;
285 1478472052 : _neighbor_subdomain = neighbor->subdomain_id();
286 1478472052 : if (_neighbor_subdomain != _old_neighbor_subdomain)
287 12818419 : neighborSubdomainChanged();
288 :
289 1478472052 : if (shouldComputeInternalSide(*elem, *neighbor))
290 739362288 : onInternalSide(elem, side);
291 :
292 1478472052 : if (boundary_ids.size() > 0)
293 769094 : for (std::vector<BoundaryID>::iterator it = boundary_ids.begin();
294 1543444 : it != boundary_ids.end();
295 774350 : ++it)
296 774350 : onInterface(elem, side, *it);
297 :
298 1478472052 : postInternalSide(elem, side);
299 : }
300 : else
301 117810427 : onExternalSide(elem, side);
302 : } // sides
303 :
304 388583668 : postElement(elem);
305 : } // range
306 :
307 4206088 : post();
308 4206088 : resetExecPrintedSets();
309 4206254 : }
310 194 : catch (libMesh::LogicError & e)
311 : {
312 28 : mooseException("We caught a libMesh error in ThreadedElementLoopBase:", e.what());
313 : }
314 0 : catch (MetaPhysicL::LogicError & e)
315 : {
316 0 : moose::translateMetaPhysicLError(e);
317 : }
318 : }
319 332 : catch (MooseException & e)
320 : {
321 166 : caughtMooseException(e);
322 : }
323 4206254 : }
324 :
325 : template <typename RangeType>
326 : void
327 4170768 : ThreadedElementLoopBase<RangeType>::pre()
328 : {
329 4170768 : }
330 :
331 : template <typename RangeType>
332 : void
333 74497 : ThreadedElementLoopBase<RangeType>::post()
334 : {
335 74497 : }
336 :
337 : template <typename RangeType>
338 : void
339 0 : ThreadedElementLoopBase<RangeType>::onElement(const Elem * /*elem*/)
340 : {
341 0 : }
342 :
343 : template <typename RangeType>
344 : void
345 130106 : ThreadedElementLoopBase<RangeType>::preElement(const Elem * /*elem*/)
346 : {
347 130106 : }
348 :
349 : template <typename RangeType>
350 : void
351 40317165 : ThreadedElementLoopBase<RangeType>::postElement(const Elem * /*elem*/)
352 : {
353 40317165 : }
354 :
355 : template <typename RangeType>
356 : void
357 40420 : ThreadedElementLoopBase<RangeType>::preBoundary(const Elem * /*elem*/,
358 : unsigned int /*side*/,
359 : BoundaryID /*bnd_id*/,
360 : const Elem * /*lower_d_elem = nullptr*/)
361 : {
362 40420 : }
363 :
364 : template <typename RangeType>
365 : void
366 5844747 : ThreadedElementLoopBase<RangeType>::onBoundary(const Elem * /*elem*/,
367 : unsigned int /*side*/,
368 : BoundaryID /*bnd_id*/,
369 : const Elem * /*lower_d_elem = nullptr*/)
370 : {
371 5844747 : }
372 :
373 : template <typename RangeType>
374 : void
375 596774 : ThreadedElementLoopBase<RangeType>::preInternalSide(const Elem * /*elem*/, unsigned int /*side*/)
376 : {
377 596774 : }
378 :
379 : template <typename RangeType>
380 : void
381 1478384292 : ThreadedElementLoopBase<RangeType>::postInternalSide(const Elem * /*elem*/, unsigned int /*side*/)
382 : {
383 1478384292 : }
384 :
385 : template <typename RangeType>
386 : void
387 45362165 : ThreadedElementLoopBase<RangeType>::onInternalSide(const Elem * /*elem*/, unsigned int /*side*/)
388 : {
389 45362165 : }
390 :
391 : template <typename RangeType>
392 : void
393 112629251 : ThreadedElementLoopBase<RangeType>::onExternalSide(const Elem * /*elem*/, unsigned int /*side*/)
394 : {
395 112629251 : }
396 :
397 : template <typename RangeType>
398 : void
399 252232 : ThreadedElementLoopBase<RangeType>::onInterface(const Elem * /*elem*/,
400 : unsigned int /*side*/,
401 : BoundaryID /*bnd_id*/)
402 : {
403 252232 : }
404 :
405 : template <typename RangeType>
406 : void
407 240959 : ThreadedElementLoopBase<RangeType>::subdomainChanged()
408 : {
409 240959 : }
410 :
411 : template <typename RangeType>
412 : void
413 4015 : ThreadedElementLoopBase<RangeType>::neighborSubdomainChanged()
414 : {
415 4015 : }
416 :
417 : template <typename RangeType>
418 : bool
419 1471845104 : ThreadedElementLoopBase<RangeType>::shouldComputeInternalSide(const Elem & elem,
420 : const Elem & neighbor) const
421 : {
422 4415535312 : auto level = [this](const auto & elem_arg)
423 : {
424 2943690208 : if (_mesh.doingPRefinement())
425 6865316 : return elem_arg.p_level();
426 : else
427 2936824892 : return elem_arg.level();
428 : };
429 1471845104 : const auto elem_id = elem.id(), neighbor_id = neighbor.id();
430 1471845104 : const auto elem_level = level(elem), neighbor_level = level(neighbor);
431 :
432 : // When looping over elements and then sides, we need to make sure that we do not duplicate
433 : // effort, e.g. if a face is shared by element 1 and element 2, then we do not want to do compute
434 : // work both when we are visiting element 1 *and* then later when visiting element 2. Our rule is
435 : // to only compute when we are visiting the element that has the lower element id when element and
436 : // neighbor are of the same adaptivity level, and then if they are not of the same level, then
437 : // we only compute when we are visiting the finer element
438 1471845104 : return (neighbor.active() && (neighbor_level == elem_level) && (elem_id < neighbor_id)) ||
439 1471845104 : (neighbor_level < elem_level);
440 : }
441 :
442 : template <typename RangeType>
443 : void
444 4206088 : ThreadedElementLoopBase<RangeType>::resetExecPrintedSets() const
445 : {
446 4206088 : _blocks_exec_printed.clear();
447 4206088 : _boundaries_exec_printed.clear();
448 4206088 : }
|
