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 : // C++
20 : #include <cstring> // for "Jacobian" exception test
21 :
22 : /**
23 : * Base class for assembly-like calculations.
24 : */
25 : template <typename RangeType>
26 : class ThreadedElementLoopBase
27 : {
28 : public:
29 : ThreadedElementLoopBase(MooseMesh & mesh);
30 :
31 : ThreadedElementLoopBase(ThreadedElementLoopBase & x, Threads::split split);
32 :
33 : virtual ~ThreadedElementLoopBase();
34 :
35 : virtual void operator()(const RangeType & range, bool bypass_threading = false);
36 :
37 : /**
38 : * Called before the element range loop
39 : */
40 : virtual void pre();
41 :
42 : /**
43 : * Called after the element range loop
44 : */
45 : virtual void post();
46 :
47 : /**
48 : * Assembly of the element (not including surface assembly)
49 : *
50 : * @param elem - active element
51 : */
52 : virtual void onElement(const Elem * elem);
53 :
54 : /**
55 : * Called before the element assembly
56 : *
57 : * @param elem - active element
58 : */
59 : virtual void preElement(const Elem * elem);
60 :
61 : /**
62 : * Called after the element assembly is done (including surface assembling)
63 : *
64 : * @param elem - active element
65 : */
66 : virtual void postElement(const Elem * elem);
67 :
68 : /**
69 : * Called before the boundary assembly
70 : *
71 : * @param elem - The element we are checking is on the boundary.
72 : * @param side - The side of the element in question.
73 : * @param bnd_id - ID of the boundary we are at
74 : * @param lower_d_elem - Lower dimensional element (e.g. Mortar)
75 : */
76 : virtual void preBoundary(const Elem * elem,
77 : unsigned int side,
78 : BoundaryID bnd_id,
79 : const Elem * lower_d_elem = nullptr);
80 :
81 : /**
82 : * Called when doing boundary assembling
83 : *
84 : * @param elem - The element we are checking is on the boundary.
85 : * @param side - The side of the element in question.
86 : * @param bnd_id - ID of the boundary we are at
87 : * @param lower_d_elem - Lower dimensional element (e.g. Mortar)
88 : */
89 : virtual void onBoundary(const Elem * elem,
90 : unsigned int side,
91 : BoundaryID bnd_id,
92 : const Elem * lower_d_elem = nullptr);
93 :
94 : /**
95 : * Called before evaluations on an element internal side
96 : *
97 : * @param elem - Element we are on
98 : * @param side - local side number of the element 'elem'
99 : */
100 : virtual void preInternalSide(const Elem * elem, unsigned int side);
101 :
102 : /**
103 : * Called after evaluations on an element internal side
104 : *
105 : * @param elem - Element we are on
106 : * @param side - local side number of the element 'elem'
107 : */
108 : virtual void postInternalSide(const Elem * elem, unsigned int side);
109 :
110 : /**
111 : * Called when doing internal edge assembling
112 : *
113 : * @param elem - Element we are on
114 : * @param side - local side number of the element 'elem'
115 : */
116 : virtual void onInternalSide(const Elem * elem, unsigned int side);
117 :
118 : /**
119 : * Called when iterating over external sides (no side neighbor)
120 : *
121 : * @param elem - Element we are on
122 : * @param side - local side number of the element 'elem'
123 : */
124 : virtual void onExternalSide(const Elem * elem, unsigned int side);
125 :
126 : /**
127 : * Called when doing interface assembling
128 : *
129 : * @param elem - Element we are on
130 : * @param side - local side number of the element 'elem'
131 : * @param bnd_id - ID of the interface we are at
132 : */
133 : virtual void onInterface(const Elem * elem, unsigned int side, BoundaryID bnd_id);
134 :
135 : /**
136 : * Called every time the current subdomain changes (i.e. the subdomain of _this_ element
137 : * is not the same as the subdomain of the last element). Beware of over-using this!
138 : * You might think that you can do some expensive stuff in here and get away with it...
139 : * but there are applications that have TONS of subdomains....
140 : */
141 : virtual void subdomainChanged();
142 :
143 : /**
144 : * Called every time the neighbor subdomain changes (i.e. the subdomain of _this_ neighbor
145 : * is not the same as the subdomain of the last neighbor). Beware of over-using this!
146 : * You might think that you can do some expensive stuff in here and get away with it...
147 : * but there are applications that have TONS of subdomains....
148 : */
149 : virtual void neighborSubdomainChanged();
150 :
151 : /**
152 : * Called if a MooseException is caught anywhere during the computation.
153 : * The single input parameter taken is a MooseException object.
154 : */
155 0 : virtual void caughtMooseException(MooseException &) {};
156 :
157 : /**
158 : * Whether or not the loop should continue.
159 : *
160 : * @return true to keep going, false to stop.
161 : */
162 120424 : virtual bool keepGoing() { return true; }
163 :
164 : protected:
165 : MooseMesh & _mesh;
166 : THREAD_ID _tid;
167 :
168 : /// The subdomain for the current element
169 : SubdomainID _subdomain;
170 :
171 : /// The subdomain for the last element
172 : SubdomainID _old_subdomain;
173 :
174 : /// The subdomain for the current neighbor
175 : SubdomainID _neighbor_subdomain;
176 :
177 : /// The subdomain for the last neighbor
178 : SubdomainID _old_neighbor_subdomain;
179 :
180 : /// Print information about the loop ordering
181 88577 : virtual void printGeneralExecutionInformation() const {}
182 :
183 : /// Print information about the particular ordering of objects on each block
184 386858 : virtual void printBlockExecutionInformation() const {}
185 :
186 : /// Print information about the particular ordering of objects on each boundary
187 12938294 : virtual void printBoundaryExecutionInformation(const unsigned int /*bid*/) const {}
188 :
189 : /// Keep track of which blocks were visited
190 : mutable std::set<SubdomainID> _blocks_exec_printed;
191 :
192 : /// Keep track of which boundaries were visited
193 : mutable std::set<BoundaryID> _boundaries_exec_printed;
194 :
195 : /// Resets the set of blocks and boundaries visited
196 : void resetExecPrintedSets() const;
197 :
198 : /**
199 : * Whether to compute the internal side for the provided element-neighbor pair. Typically this
200 : * will return true if the element id is less than the neighbor id when the elements are equal
201 : * level, or when the element is more refined than the neighbor, and then false otherwise. One
202 : * type of loop where the logic will be different is when projecting stateful material properties
203 : */
204 : virtual bool shouldComputeInternalSide(const Elem & elem, const Elem & neighbor) const;
205 : };
206 :
207 : template <typename RangeType>
208 3949533 : ThreadedElementLoopBase<RangeType>::ThreadedElementLoopBase(MooseMesh & mesh) : _mesh(mesh)
209 : {
210 3949533 : }
211 :
212 : template <typename RangeType>
213 61 : ThreadedElementLoopBase<RangeType>::ThreadedElementLoopBase(ThreadedElementLoopBase & x,
214 : Threads::split /*split*/)
215 61 : : _mesh(x._mesh)
216 : {
217 61 : }
218 :
219 : template <typename RangeType>
220 4307849 : ThreadedElementLoopBase<RangeType>::~ThreadedElementLoopBase()
221 : {
222 4307849 : }
223 :
224 : template <typename RangeType>
225 : void
226 4173104 : ThreadedElementLoopBase<RangeType>::operator()(const RangeType & range, bool bypass_threading)
227 : {
228 : try
229 : {
230 : try
231 : {
232 4173104 : ParallelUniqueId puid;
233 4173104 : _tid = bypass_threading ? 0 : puid.id;
234 :
235 4173104 : pre();
236 4173104 : printGeneralExecutionInformation();
237 :
238 4173104 : _subdomain = Moose::INVALID_BLOCK_ID;
239 4173104 : _neighbor_subdomain = Moose::INVALID_BLOCK_ID;
240 4173104 : typename RangeType::const_iterator el = range.begin();
241 761884948 : for (el = range.begin(); el != range.end(); ++el)
242 : {
243 378872985 : if (!keepGoing())
244 10 : break;
245 :
246 378872975 : const Elem * elem = *el;
247 :
248 378872975 : preElement(elem);
249 :
250 378872975 : _old_subdomain = _subdomain;
251 378872975 : _subdomain = elem->subdomain_id();
252 378872975 : if (_subdomain != _old_subdomain)
253 : {
254 5483929 : subdomainChanged();
255 5483929 : printBlockExecutionInformation();
256 : }
257 :
258 378872975 : onElement(elem);
259 :
260 757729402 : if (_mesh.interiorLowerDBlocks().count(elem->subdomain_id()) > 0 ||
261 757729402 : _mesh.boundaryLowerDBlocks().count(elem->subdomain_id()) > 0)
262 : {
263 33664 : postElement(elem);
264 33664 : continue;
265 : }
266 :
267 378839096 : const auto elem_boundary_ids = _mesh.getBoundaryIDs(elem);
268 1949002430 : for (unsigned int side = 0; side < elem->n_sides(); side++)
269 : {
270 1570163340 : const auto & boundary_ids = elem_boundary_ids[side];
271 1570163340 : const Elem * lower_d_elem = _mesh.getLowerDElem(elem, side);
272 :
273 1687782689 : for (const auto bnd_id : boundary_ids)
274 : {
275 117619355 : preBoundary(elem, side, bnd_id, lower_d_elem);
276 117619355 : printBoundaryExecutionInformation(bnd_id);
277 117619355 : onBoundary(elem, side, bnd_id, lower_d_elem);
278 : }
279 :
280 1570163334 : const Elem * neighbor = elem->neighbor_ptr(side);
281 1570163334 : if (neighbor)
282 : {
283 1450991432 : preInternalSide(elem, side);
284 :
285 1450991432 : _old_neighbor_subdomain = _neighbor_subdomain;
286 1450991432 : _neighbor_subdomain = neighbor->subdomain_id();
287 1450991432 : if (_neighbor_subdomain != _old_neighbor_subdomain)
288 14000149 : neighborSubdomainChanged();
289 :
290 1450991432 : if (shouldComputeInternalSide(*elem, *neighbor))
291 90795762 : onInternalSide(elem, side);
292 :
293 1451797826 : for (const auto bnd_id : boundary_ids)
294 806394 : onInterface(elem, side, bnd_id);
295 :
296 1450991432 : postInternalSide(elem, side);
297 : }
298 : else
299 119171902 : onExternalSide(elem, side);
300 : } // sides
301 :
302 378839090 : postElement(elem);
303 : } // range
304 :
305 4172883 : post();
306 4172883 : resetExecPrintedSets();
307 4173040 : }
308 157 : catch (MetaPhysicL::LogicError & e)
309 : {
310 0 : moose::translateMetaPhysicLError(e);
311 : }
312 314 : catch (std::exception & e)
313 : {
314 : // Continue if we find a libMesh degenerate map exception, but
315 : // just re-throw for any real error
316 271 : if (!strstr(e.what(), "Jacobian") && !strstr(e.what(), "singular") &&
317 114 : !strstr(e.what(), "det != 0"))
318 114 : throw; // not "throw e;" - that destroys type info!
319 :
320 43 : mooseException("We caught a libMesh degeneracy exception in ThreadedElementLoopBase:\n",
321 : e.what());
322 : }
323 : }
324 314 : catch (MooseException & e)
325 : {
326 157 : caughtMooseException(e);
327 : }
328 4173040 : }
329 :
330 : template <typename RangeType>
331 : void
332 4137806 : ThreadedElementLoopBase<RangeType>::pre()
333 : {
334 4137806 : }
335 :
336 : template <typename RangeType>
337 : void
338 78751 : ThreadedElementLoopBase<RangeType>::post()
339 : {
340 78751 : }
341 :
342 : template <typename RangeType>
343 : void
344 0 : ThreadedElementLoopBase<RangeType>::onElement(const Elem * /*elem*/)
345 : {
346 0 : }
347 :
348 : template <typename RangeType>
349 : void
350 120424 : ThreadedElementLoopBase<RangeType>::preElement(const Elem * /*elem*/)
351 : {
352 120424 : }
353 :
354 : template <typename RangeType>
355 : void
356 41100447 : ThreadedElementLoopBase<RangeType>::postElement(const Elem * /*elem*/)
357 : {
358 41100447 : }
359 :
360 : template <typename RangeType>
361 : void
362 38480 : ThreadedElementLoopBase<RangeType>::preBoundary(const Elem * /*elem*/,
363 : unsigned int /*side*/,
364 : BoundaryID /*bnd_id*/,
365 : const Elem * /*lower_d_elem = nullptr*/)
366 : {
367 38480 : }
368 :
369 : template <typename RangeType>
370 : void
371 5956030 : ThreadedElementLoopBase<RangeType>::onBoundary(const Elem * /*elem*/,
372 : unsigned int /*side*/,
373 : BoundaryID /*bnd_id*/,
374 : const Elem * /*lower_d_elem = nullptr*/)
375 : {
376 5956030 : }
377 :
378 : template <typename RangeType>
379 : void
380 551844 : ThreadedElementLoopBase<RangeType>::preInternalSide(const Elem * /*elem*/, unsigned int /*side*/)
381 : {
382 551844 : }
383 :
384 : template <typename RangeType>
385 : void
386 1450893384 : ThreadedElementLoopBase<RangeType>::postInternalSide(const Elem * /*elem*/, unsigned int /*side*/)
387 : {
388 1450893384 : }
389 :
390 : template <typename RangeType>
391 : void
392 45450500 : ThreadedElementLoopBase<RangeType>::onInternalSide(const Elem * /*elem*/, unsigned int /*side*/)
393 : {
394 45450500 : }
395 :
396 : template <typename RangeType>
397 : void
398 7052096 : ThreadedElementLoopBase<RangeType>::onExternalSide(const Elem * /*elem*/, unsigned int /*side*/)
399 : {
400 7052096 : }
401 :
402 : template <typename RangeType>
403 : void
404 250943 : ThreadedElementLoopBase<RangeType>::onInterface(const Elem * /*elem*/,
405 : unsigned int /*side*/,
406 : BoundaryID /*bnd_id*/)
407 : {
408 250943 : }
409 :
410 : template <typename RangeType>
411 : void
412 298822 : ThreadedElementLoopBase<RangeType>::subdomainChanged()
413 : {
414 298822 : }
415 :
416 : template <typename RangeType>
417 : void
418 3935 : ThreadedElementLoopBase<RangeType>::neighborSubdomainChanged()
419 : {
420 3935 : }
421 :
422 : template <typename RangeType>
423 : bool
424 177820683 : ThreadedElementLoopBase<RangeType>::shouldComputeInternalSide(const Elem & elem,
425 : const Elem & neighbor) const
426 : {
427 : // If we're going to compute the internal side with this elem-neighbor pair, then they must both
428 : // be active. Note that if elem is an active coarse element at an interface with finer elements,
429 : // then its neighbor will be an equal-level inactive element, hence the following
430 : // neighbor.active() check. In that case we'll catch current 'elem' when we come around and
431 : // examine one of the finer elements as 'elem'
432 : mooseAssert(elem.active(), "This method should never be called with an inactive element");
433 177820683 : if (!neighbor.active())
434 486373 : return false;
435 :
436 : //
437 : // Define an ordering: first prefer finer by h (higher level()), then prefer finer by p (higher
438 : // p_level()), then prefer smaller id()"
439 : //
440 :
441 177334310 : const auto elem_level = elem.level(), neighbor_level = neighbor.level();
442 177334310 : if (elem_level != neighbor_level)
443 1135176 : return elem_level > neighbor_level;
444 :
445 176199134 : const auto elem_p_level = elem.p_level(), neighbor_p_level = neighbor.p_level();
446 176199134 : if (elem_p_level != neighbor_p_level)
447 85418 : return elem_p_level > neighbor_p_level;
448 :
449 176113716 : return elem.id() < neighbor.id();
450 : }
451 :
452 : template <typename RangeType>
453 : void
454 4172883 : ThreadedElementLoopBase<RangeType>::resetExecPrintedSets() const
455 : {
456 4172883 : _blocks_exec_printed.clear();
457 4172883 : _boundaries_exec_printed.clear();
458 4172883 : }
|
