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