https://mooseframework.inl.gov
DiracKernelBase.C
Go to the documentation of this file.
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 
21 {
26 
27  params.addParam<bool>("use_displaced_mesh",
28  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  params.addParam<bool>(
34  "drop_duplicate_points",
35  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  MooseEnum point_not_found_behavior("ERROR WARNING IGNORE", "IGNORE");
41  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  params.addParam<bool>("allow_moving_sources",
49  false,
50  "If true, allow Dirac sources to move, even if the mesh does not move, "
51  "during the simulation.");
52 
53  params.addParamNamesToGroup("use_displaced_mesh drop_duplicate_points", "Advanced");
54  params.declareControllable("enable");
55  params.registerSystemAttributeName("DiracKernel");
56 
57  return params;
58 }
59 
61  : ResidualObject(parameters),
65  BlockRestrictable(this),
66  _current_elem(_assembly.elem()),
67  _coord_sys(_assembly.coordSystem()),
68  _dirac_kernel_info(_subproblem.diracKernelInfo()),
69  _q_point(_assembly.qPoints()),
70  _physical_point(_assembly.physicalPoints()),
71  _qrule(_assembly.qRule()),
72  _JxW(_assembly.JxW()),
73  _drop_duplicate_points(parameters.get<bool>("drop_duplicate_points")),
74  _point_not_found_behavior(
75  parameters.get<MooseEnum>("point_not_found_behavior").getEnum<PointNotFoundBehavior>()),
76  _allow_moving_sources(getParam<bool>("allow_moving_sources"))
77 {
78  // Stateful material properties are not allowed on DiracKernels
80 }
81 
82 Real
84 {
85  return 0;
86 }
87 
88 void
89 DiracKernelBase::addPoint(const Elem * elem, Point p, unsigned /*id*/)
90 {
91  if (!elem || !hasBlocks(elem->subdomain_id()))
92  {
93  std::stringstream msg;
94  msg << "Point " << p << " not found in block(s) " << Moose::stringify(blockIDs(), ", ") << ".";
96  {
98  mooseError(msg.str());
99  break;
101  mooseDoOnce(mooseWarning(msg.str()));
102  break;
104  break;
105  default:
106  mooseError("Internal enum error.");
107  }
108  return;
109  }
110 
111  if (elem->processor_id() != processor_id())
112  return;
113 
114  _dirac_kernel_info.addPoint(elem, p);
116 }
117 
118 const Elem *
119 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  if (id != libMesh::invalid_uint)
127  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  const Elem * elem = _dirac_kernel_info.findPoint(p, _mesh, blockIDs());
133  addPoint(elem, p, id);
134  return elem;
135 }
136 
137 const Elem *
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  const Elem * return_elem = NULL;
144 
145  // May be set if the Elem is found in our cache, otherwise stays as NULL.
146  const Elem * cached_elem = NULL;
147 
148  // OK, the user gave us an ID, let's see if we already have it...
149  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  unsigned int i_found_it = static_cast<unsigned int>(it != _point_cache.end());
153  unsigned int we_found_it = i_found_it;
154  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  if (!we_found_it)
160  {
161  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  if (elem && (elem->processor_id() == processor_id()))
165  {
166  // Add the point to the cache...
167  _point_cache[id] = std::make_pair(elem, p);
168 
169  // ... and to the reverse cache.
170  std::vector<std::pair<Point, unsigned>> & points = _reverse_point_cache[elem];
171  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  addPoint(elem, p, id);
177  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  if (we_found_it && !i_found_it)
186  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  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  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  Point cached_point = (it->second).second;
201 
202  if (cached_point.relative_fuzzy_equals(p))
203  {
204  // Find the cached element associated to this point
205  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  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  updateCaches(cached_elem, NULL, p, id);
216  return_elem = NULL;
217  break; // out of while loop
218  }
219 
220  bool active = cached_elem->active();
221  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  if (active && contains_point)
227  {
228  addPoint(cached_elem, p, id);
229  return_elem = cached_elem;
230  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  else if (!active && contains_point)
236  {
237  // Get the list of active children
238  std::vector<const Elem *> active_children;
239  cached_elem->active_family_tree(active_children);
240 
241  // Linear search through active children for the one that contains p
242  for (unsigned c = 0; c < active_children.size(); ++c)
243  if (active_children[c]->contains_point(p))
244  {
245  updateCaches(cached_elem, active_children[c], p, id);
246  addPoint(active_children[c], p, id);
247  return_elem = active_children[c];
248  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  if (!return_elem)
257  mooseError("Error, Point not found in any of the active children!");
258 
259  break; // out of while loop
260  }
261 
262  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  (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  (!active && !contains_point))
277  {
278  i_need_find_point = true;
279  break; // out of while loop
280  }
281 
282  else
283  mooseError("We'll never get here!");
284  } // if (cached_point.relative_fuzzy_equals(p))
285  else
286  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  unsigned int we_need_find_point = static_cast<unsigned int>(i_need_find_point);
303  comm().max(we_need_find_point);
304 
305  if (we_need_find_point)
306  {
307  // findPoint() is a parallel-only function
308  const Elem * elem = _dirac_kernel_info.findPoint(p, _mesh, blockIDs());
309 
310  updateCaches(cached_elem, elem, p, id);
311  addPoint(elem, p, id);
312  return_elem = elem;
313  }
314 
315  return return_elem;
316 }
317 
318 bool
320 {
321  return _local_dirac_kernel_info.getElements().count(_mesh.elemPtr(elem->id())) != 0;
322 }
323 
324 bool
325 DiracKernelBase::isActiveAtPoint(const Elem * elem, const Point & p)
326 {
327  return _local_dirac_kernel_info.hasPoint(elem, p);
328 }
329 
330 void
332 {
334 
335  // If points can move, we can't trust caches
338 }
339 
340 void
342 {
343  _point_cache.clear();
344  _reverse_point_cache.clear();
345 }
346 
347 void
348 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  _point_cache.erase(id);
353  if (new_elem && (new_elem->processor_id() == processor_id()))
354  _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  reverse_cache_t::iterator it = _reverse_point_cache.find(old_elem);
360  if (it != _reverse_point_cache.end())
361  {
362  reverse_cache_t::mapped_type & points = it->second;
363  {
364  reverse_cache_t::mapped_type::iterator points_it = points.begin(), points_end = points.end();
365 
366  for (; points_it != points_end; ++points_it)
367  {
368  // If the point matches, remove it from the vector of points
369  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  points.erase(points_it);
376  break;
377  }
378  }
379 
380  // If the points vector is now empty, remove old_elem from the reverse cache entirely
381  if (points.empty())
382  _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  if (new_elem && (new_elem->processor_id() == processor_id()))
388  {
389  reverse_cache_t::mapped_type & points = _reverse_point_cache[new_elem];
390  points.push_back(std::make_pair(p, id));
391  }
392 }
393 
394 unsigned
396 {
397  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  if (it == _reverse_point_cache.end())
401  return libMesh::invalid_uint;
402 
403  // Do a linear search in the (hopefully small) vector of Points for this Elem
404  reverse_cache_t::mapped_type & points = it->second;
405 
406  for (const auto & points_it : points)
407  {
408  // If the current_point equals the cached point, return the associated id
409  if (_current_point.relative_fuzzy_equals(points_it.first))
410  return points_it.second;
411  }
412 
413  // If we made it here, we didn't find the cached point, so return invalid_uint
414  return libMesh::invalid_uint;
415 }
MooseMesh & _mesh
Reference to this Kernel&#39;s mesh object.
const bool _allow_moving_sources
Whether Dirac sources can move during the simulation.
bool hasPointsOnElem(const Elem *elem)
Whether or not this DiracKernel has something to distribute on this element.
const PointNotFoundBehavior _point_not_found_behavior
bool hasPoint(const Elem *elem, const Point &p)
Return true if we have Point &#39;p&#39; in Element &#39;elem&#39;.
const unsigned int invalid_uint
const Elem * findPoint(const Point &p, const MooseMesh &mesh, const std::set< SubdomainID > &blocks)
Used by client DiracKernel classes to determine the Elem in which the Point p resides.
virtual Elem * elemPtr(const dof_id_type i)
Definition: MooseMesh.C:3153
T * get(const std::unique_ptr< T > &u)
The MooseUtils::get() specializations are used to support making forwards-compatible code changes fro...
Definition: MooseUtils.h:1133
void clearPoints()
Remove all of the current points and elements.
DiracKernelInfo & _dirac_kernel_info
Place for storing Point/Elem information shared across all DiracKernel objects.
point_cache_t _point_cache
void addPoint(const Elem *elem, Point p, unsigned id=libMesh::invalid_uint)
Add the physical x,y,z point located in the element "elem" to the list of points this DiracKernel wil...
void registerSystemAttributeName(const std::string &value)
This method is used to define the MOOSE system name that is used by the TheWarehouse object for stori...
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
const Parallel::Communicator & comm() const
unsigned currentPointCachedID()
Returns the user-assigned ID of the current Dirac point if it exits, and libMesh::invalid_uint otherw...
static InputParameters validParams()
virtual const std::set< SubdomainID > & blockIDs() const
Return the block subdomain ids for this object Note, if this is not block restricted, this function returns all mesh subdomain ids.
static InputParameters validParams()
void addPoint(const Elem *elem, const Point &p)
Adds a point source.
virtual Real computeQpOffDiagJacobian(unsigned int jvar)
This gets called by computeOffDiagJacobian() at each quadrature point.
DiracKernelBase(const InputParameters &parameters)
static InputParameters validParams()
bool isActiveAtPoint(const Elem *elem, const Point &p)
Whether or not this DiracKernel has something to distribute at this Point.
libmesh_assert(ctx)
This is a "smart" enum class intended to replace many of the shortcomings in the C++ enum type It sho...
Definition: MooseEnum.h:33
static InputParameters validParams()
DiracKernelInfo _local_dirac_kernel_info
Place for storing Point/Elem information only for this DiracKernel.
void clearPoints()
Remove all of the current points and elements.
std::string stringify(const T &t)
conversion to string
Definition: Conversion.h:64
const Elem *const & _current_elem
reverse_cache_t _reverse_point_cache
void statefulPropertiesAllowed(bool)
Derived classes can declare whether or not they work with stateful material properties.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
This is the common base class for objects that give residual contributions.
const std::set< SubdomainID > EMPTY_BLOCK_IDS
Definition: MooseTypes.h:684
void max(const T &r, T &o, Request &req) const
An interface for accessing Materials.
const std::set< BoundaryID > EMPTY_BOUNDARY_IDS
Definition: MooseTypes.h:685
Intermediate base class that ties together all the interfaces for getting MooseVariableFEBases with t...
void clearPointsCaches()
Clear the cache of points because the points may have moved.
An interface that restricts an object to subdomains via the &#39;blocks&#39; input parameter.
void mooseWarning(Args &&... args) const
Emits a warning prefixed with object name and type.
Definition: MooseBase.h:299
void mooseError(Args &&... args) const
Emits an error prefixed with object name and type and optionally a file path to the top-level block p...
Definition: MooseBase.h:271
void updateCaches(const Elem *old_elem, const Elem *new_elem, Point p, unsigned id)
This function is used internally when the Elem for a locally-cached point needs to be updated...
void addParam(const std::string &name, const S &value, const std::string &doc_string)
These methods add an optional parameter and a documentation string to the InputParameters object...
std::set< const Elem * > & getElements()
Returns a writeable reference to the _elements container.
MOOSE now contains C++17 code, so give a reasonable error message stating what the user can do to add...
bool hasBlocks(const SubdomainName &name) const
Test if the supplied block name is valid for this object.
static InputParameters validParams()
processor_id_type processor_id() const
void declareControllable(const std::string &name, std::set< ExecFlagType > execute_flags={})
Declare the given parameters as controllable.
const Elem * addPointWithValidId(Point p, unsigned id)
A helper function for addPoint(Point, id) for when id != invalid_uint.
Point _current_point
The current point.
void addParamNamesToGroup(const std::string &space_delim_names, const std::string group_name)
This method takes a space delimited list of parameter names and adds them to the specified group name...