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