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