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 : // Moose includes
11 : #include "DiracKernelBase.h"
12 : #include "Assembly.h"
13 : #include "SystemBase.h"
14 : #include "Problem.h"
15 : #include "MooseMesh.h"
16 :
17 : #include "libmesh/quadrature.h"
18 :
19 : InputParameters
20 181774 : DiracKernelBase::validParams()
21 : {
22 181774 : InputParameters params = ResidualObject::validParams();
23 181774 : params += MaterialPropertyInterface::validParams();
24 181774 : params += BlockRestrictable::validParams();
25 181774 : params += GeometricSearchInterface::validParams();
26 :
27 545322 : params.addParam<bool>("use_displaced_mesh",
28 363548 : false,
29 : "Whether or not this object should use the displaced mesh for computation. "
30 : "Note that in the case this is true but no displacements are provided in "
31 : "the Mesh block the undisplaced mesh will still be used.");
32 :
33 545322 : params.addParam<bool>(
34 : "drop_duplicate_points",
35 363548 : true,
36 : "By default points added to a DiracKernel are dropped if a point at the same location"
37 : "has been added before. If this option is set to false duplicate points are retained"
38 : "and contribute to residual and Jacobian.");
39 :
40 727096 : MooseEnum point_not_found_behavior("ERROR WARNING IGNORE", "IGNORE");
41 727096 : params.addParam<MooseEnum>(
42 : "point_not_found_behavior",
43 : point_not_found_behavior,
44 : "By default (IGNORE), it is ignored if an added point cannot be located in the "
45 : "specified subdomains. If this option is set to ERROR, this situation will result in an "
46 : "error. If this option is set to WARNING, then a warning will be issued.");
47 :
48 545322 : params.addParam<bool>("allow_moving_sources",
49 363548 : false,
50 : "If true, allow Dirac sources to move, even if the mesh does not move, "
51 : "during the simulation.");
52 :
53 545322 : params.addParamNamesToGroup("use_displaced_mesh drop_duplicate_points", "Advanced");
54 545322 : params.declareControllable("enable");
55 181774 : params.registerSystemAttributeName("DiracKernel");
56 :
57 363548 : return params;
58 181774 : }
59 :
60 1024 : DiracKernelBase::DiracKernelBase(const InputParameters & parameters)
61 : : ResidualObject(parameters),
62 : CoupleableMooseVariableDependencyIntermediateInterface(this, false),
63 : MaterialPropertyInterface(this, Moose::EMPTY_BLOCK_IDS, Moose::EMPTY_BOUNDARY_IDS),
64 : GeometricSearchInterface(this),
65 : BlockRestrictable(this),
66 2048 : _current_elem(_assembly.elem()),
67 1024 : _coord_sys(_assembly.coordSystem()),
68 1024 : _dirac_kernel_info(_subproblem.diracKernelInfo()),
69 1024 : _q_point(_assembly.qPoints()),
70 1024 : _physical_point(_assembly.physicalPoints()),
71 1024 : _qrule(_assembly.qRule()),
72 1024 : _JxW(_assembly.JxW()),
73 1024 : _drop_duplicate_points(parameters.get<bool>("drop_duplicate_points")),
74 1024 : _point_not_found_behavior(
75 1024 : parameters.get<MooseEnum>("point_not_found_behavior").getEnum<PointNotFoundBehavior>()),
76 5120 : _allow_moving_sources(getParam<bool>("allow_moving_sources"))
77 : {
78 : // Stateful material properties are not allowed on DiracKernels
79 1024 : statefulPropertiesAllowed(false);
80 1024 : }
81 :
82 : void
83 408759 : DiracKernelBase::addPoint(const Elem * elem, Point p, unsigned /*id*/)
84 : {
85 408759 : if (!elem || !hasBlocks(elem->subdomain_id()))
86 : {
87 13164 : std::stringstream msg;
88 52656 : msg << "Point " << p << " not found in block(s) " << Moose::stringify(blockIDs(), ", ") << ".";
89 13164 : switch (_point_not_found_behavior)
90 : {
91 4 : case PointNotFoundBehavior::ERROR:
92 4 : mooseError(msg.str());
93 : break;
94 329 : case PointNotFoundBehavior::WARNING:
95 329 : mooseDoOnce(mooseWarning(msg.str()));
96 329 : break;
97 12831 : case PointNotFoundBehavior::IGNORE:
98 12831 : break;
99 0 : default:
100 0 : mooseError("Internal enum error.");
101 : }
102 13160 : return;
103 13160 : }
104 :
105 395595 : if (elem->processor_id() != processor_id())
106 0 : return;
107 :
108 395595 : _dirac_kernel_info.addPoint(elem, p);
109 395595 : _local_dirac_kernel_info.addPoint(elem, p);
110 : }
111 :
112 : const Elem *
113 82139 : DiracKernelBase::addPoint(Point p, unsigned id)
114 : {
115 : // Make sure that this method was called with the same id on all
116 : // processors. It's an extra communication, though, so let's only
117 : // do it in DEBUG mode.
118 : libmesh_assert(comm().verify(id));
119 :
120 82139 : if (id != libMesh::invalid_uint)
121 42722 : return addPointWithValidId(p, id);
122 :
123 : // If id == libMesh::invalid_uint (the default), the user is not
124 : // enabling caching when they add Dirac points. So all we can do is
125 : // the PointLocator lookup, and call the other addPoint() method.
126 39417 : const Elem * elem = _dirac_kernel_info.findPoint(p, _mesh, blockIDs());
127 39417 : addPoint(elem, p, id);
128 39417 : return elem;
129 : }
130 :
131 : const Elem *
132 42722 : DiracKernelBase::addPointWithValidId(Point p, unsigned id)
133 : {
134 : // The Elem we'll eventually return. We can't return early on some
135 : // processors, because we may need to call parallel_only() functions in
136 : // the remainder of this scope.
137 42722 : const Elem * return_elem = NULL;
138 :
139 : // May be set if the Elem is found in our cache, otherwise stays as NULL.
140 42722 : const Elem * cached_elem = NULL;
141 :
142 : // OK, the user gave us an ID, let's see if we already have it...
143 42722 : point_cache_t::iterator it = _point_cache.find(id);
144 :
145 : // Was the point found in a _point_cache on at least one processor?
146 42722 : unsigned int i_found_it = static_cast<unsigned int>(it != _point_cache.end());
147 42722 : unsigned int we_found_it = i_found_it;
148 42722 : comm().max(we_found_it);
149 :
150 : // If nobody found it in their local caches, it means we need to
151 : // do the PointLocator look-up and update the caches. This is
152 : // safe, because all processors have the same value of we_found_it.
153 42722 : if (!we_found_it)
154 : {
155 4190 : const Elem * elem = _dirac_kernel_info.findPoint(p, _mesh, blockIDs());
156 :
157 : // Only add the point to the cache on this processor if the Elem is local
158 4190 : if (elem && (elem->processor_id() == processor_id()))
159 : {
160 : // Add the point to the cache...
161 642 : _point_cache[id] = std::make_pair(elem, p);
162 :
163 : // ... and to the reverse cache.
164 642 : std::vector<std::pair<Point, unsigned>> & points = _reverse_point_cache[elem];
165 642 : points.push_back(std::make_pair(p, id));
166 : }
167 :
168 : // Call the other addPoint() method. This method ignores non-local
169 : // and NULL elements automatically.
170 4190 : addPoint(elem, p, id);
171 4186 : return_elem = elem;
172 : }
173 :
174 : // If the point was found in a cache, but not my cache, I'm not
175 : // responsible for it.
176 : //
177 : // We can't return early here: then we aren't allowed to call any more
178 : // parallel_only() functions in the remainder of this function!
179 42718 : if (we_found_it && !i_found_it)
180 10066 : return_elem = NULL;
181 :
182 : // This flag may be set by the processor that cached the Elem because it
183 : // needs to call findPoint() (due to moving mesh, etc.). If so, we will
184 : // call it at the end of the while loop below.
185 42718 : bool i_need_find_point = false;
186 :
187 : // Now that we only cache local data, some processors may enter
188 : // this if statement and some may not. Therefore we can't call
189 : // any parallel_only() functions inside this if statement.
190 42718 : while (i_found_it)
191 : {
192 : // We have something cached, now make sure it's actually the same Point.
193 : // TODO: we should probably use this same comparison in the DiracKernelInfo code!
194 28466 : Point cached_point = (it->second).second;
195 :
196 28466 : if (cached_point.relative_fuzzy_equals(p))
197 : {
198 : // Find the cached element associated to this point
199 28462 : cached_elem = (it->second).first;
200 :
201 : // If the cached element's processor ID doesn't match ours, we
202 : // are no longer responsible for caching it. This can happen
203 : // due to adaptivity...
204 28462 : if (cached_elem->processor_id() != processor_id())
205 : {
206 : // Update the caches, telling them to drop the cached Elem.
207 : // Analogously to the rest of the DiracKernel system, we
208 : // also return NULL because the Elem is non-local.
209 0 : updateCaches(cached_elem, NULL, p, id);
210 0 : return_elem = NULL;
211 28462 : break; // out of while loop
212 : }
213 :
214 28462 : bool active = cached_elem->active();
215 28462 : bool contains_point = cached_elem->contains_point(p);
216 :
217 : // If the cached Elem is active and the point is still
218 : // contained in it, call the other addPoint() method and
219 : // return its result.
220 28462 : if (active && contains_point)
221 : {
222 28366 : addPoint(cached_elem, p, id);
223 28366 : return_elem = cached_elem;
224 28366 : break; // out of while loop
225 : }
226 :
227 : // Is the Elem not active (been refined) but still contains the point?
228 : // Then search in its active children and update the caches.
229 96 : else if (!active && contains_point)
230 : {
231 : // Get the list of active children
232 0 : std::vector<const Elem *> active_children;
233 0 : cached_elem->active_family_tree(active_children);
234 :
235 : // Linear search through active children for the one that contains p
236 0 : for (unsigned c = 0; c < active_children.size(); ++c)
237 0 : if (active_children[c]->contains_point(p))
238 : {
239 0 : updateCaches(cached_elem, active_children[c], p, id);
240 0 : addPoint(active_children[c], p, id);
241 0 : return_elem = active_children[c];
242 0 : break; // out of for loop
243 : }
244 :
245 : // If we got here without setting return_elem, it means the Point was
246 : // found in the parent element, but not in any of the active
247 : // children... this is not possible under normal
248 : // circumstances, so something must have gone seriously
249 : // wrong!
250 0 : if (!return_elem)
251 0 : mooseError("Error, Point not found in any of the active children!");
252 :
253 0 : break; // out of while loop
254 0 : }
255 :
256 96 : else if (
257 : // Is the Elem active but the point is not contained in it any
258 : // longer? (For example, did the Mesh move out from under
259 : // it?) Then we fall back to the expensive Point Locator
260 : // lookup. TODO: we could try and do something more optimized
261 : // like checking if any of the active neighbors contains the
262 : // point. Update the caches.
263 96 : (active && !contains_point) ||
264 :
265 : // The Elem has been refined *and* the Mesh has moved out
266 : // from under it, we fall back to doing the expensive Point
267 : // Locator lookup. TODO: We could try and look in the
268 : // active children of this Elem's neighbors for the Point.
269 : // Update the caches.
270 0 : (!active && !contains_point))
271 : {
272 96 : i_need_find_point = true;
273 96 : break; // out of while loop
274 : }
275 :
276 : else
277 0 : mooseError("We'll never get here!");
278 : } // if (cached_point.relative_fuzzy_equals(p))
279 : else
280 4 : mooseError("Cached Dirac point ",
281 : cached_point,
282 : " already exists with ID: ",
283 : id,
284 : " and does not match point ",
285 : p,
286 : "If Dirac sources are moving, please set 'allow_moving_sources' to true");
287 :
288 : // We only want one iteration of this while loop at maximum.
289 : i_found_it = false;
290 : } // while (i_found_it)
291 :
292 : // We are back to all processors here because we do not return
293 : // early in the code above...
294 :
295 : // Does we need to call findPoint() on all processors.
296 42714 : unsigned int we_need_find_point = static_cast<unsigned int>(i_need_find_point);
297 42714 : comm().max(we_need_find_point);
298 :
299 42714 : if (we_need_find_point)
300 : {
301 : // findPoint() is a parallel-only function
302 126 : const Elem * elem = _dirac_kernel_info.findPoint(p, _mesh, blockIDs());
303 :
304 126 : updateCaches(cached_elem, elem, p, id);
305 126 : addPoint(elem, p, id);
306 126 : return_elem = elem;
307 : }
308 :
309 42714 : return return_elem;
310 : }
311 :
312 : bool
313 361735 : DiracKernelBase::hasPointsOnElem(const Elem * elem)
314 : {
315 361735 : return _local_dirac_kernel_info.getElements().count(_mesh.elemPtr(elem->id())) != 0;
316 : }
317 :
318 : bool
319 355937 : DiracKernelBase::isActiveAtPoint(const Elem * elem, const Point & p)
320 : {
321 355937 : return _local_dirac_kernel_info.hasPoint(elem, p);
322 : }
323 :
324 : void
325 55361 : DiracKernelBase::clearPoints()
326 : {
327 55361 : _local_dirac_kernel_info.clearPoints();
328 :
329 : // If points can move, we can't trust caches
330 55361 : if (_allow_moving_sources)
331 0 : clearPointsCaches();
332 55361 : }
333 :
334 : void
335 89 : DiracKernelBase::clearPointsCaches()
336 : {
337 89 : _point_cache.clear();
338 89 : _reverse_point_cache.clear();
339 89 : }
340 :
341 : void
342 126 : DiracKernelBase::updateCaches(const Elem * old_elem, const Elem * new_elem, Point p, unsigned id)
343 : {
344 : // Update the point cache. Remove old cached data, only cache
345 : // new_elem if it is non-NULL and local.
346 126 : _point_cache.erase(id);
347 126 : if (new_elem && (new_elem->processor_id() == processor_id()))
348 78 : _point_cache[id] = std::make_pair(new_elem, p);
349 :
350 : // Update the reverse cache
351 : //
352 : // First, remove the Point from the old_elem's vector
353 126 : reverse_cache_t::iterator it = _reverse_point_cache.find(old_elem);
354 126 : if (it != _reverse_point_cache.end())
355 : {
356 96 : reverse_cache_t::mapped_type & points = it->second;
357 : {
358 96 : reverse_cache_t::mapped_type::iterator points_it = points.begin(), points_end = points.end();
359 :
360 96 : for (; points_it != points_end; ++points_it)
361 : {
362 : // If the point matches, remove it from the vector of points
363 96 : if (p.relative_fuzzy_equals(points_it->first))
364 : {
365 : // Vector erasure. It can be slow but these vectors are
366 : // generally very short. It also invalidates existing
367 : // iterators, so we need to break out of the loop and
368 : // not use them any more.
369 96 : points.erase(points_it);
370 96 : break;
371 : }
372 : }
373 :
374 : // If the points vector is now empty, remove old_elem from the reverse cache entirely
375 96 : if (points.empty())
376 96 : _reverse_point_cache.erase(old_elem);
377 : }
378 : }
379 :
380 : // Next, if new_elem is not NULL and local, add the point to the new_elem's vector
381 126 : if (new_elem && (new_elem->processor_id() == processor_id()))
382 : {
383 78 : reverse_cache_t::mapped_type & points = _reverse_point_cache[new_elem];
384 78 : points.push_back(std::make_pair(p, id));
385 : }
386 126 : }
387 :
388 : unsigned
389 63232 : DiracKernelBase::currentPointCachedID()
390 : {
391 63232 : reverse_cache_t::iterator it = _reverse_point_cache.find(_current_elem);
392 :
393 : // If the current Elem is not in the cache, return invalid_uint
394 63232 : if (it == _reverse_point_cache.end())
395 0 : return libMesh::invalid_uint;
396 :
397 : // Do a linear search in the (hopefully small) vector of Points for this Elem
398 63232 : reverse_cache_t::mapped_type & points = it->second;
399 :
400 63232 : for (const auto & points_it : points)
401 : {
402 : // If the current_point equals the cached point, return the associated id
403 63232 : if (_current_point.relative_fuzzy_equals(points_it.first))
404 63232 : return points_it.second;
405 : }
406 :
407 : // If we made it here, we didn't find the cached point, so return invalid_uint
408 0 : return libMesh::invalid_uint;
409 : }
|
