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 {
25 
26  params.addParam<bool>("use_displaced_mesh",
27  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  params.addParam<bool>(
33  "drop_duplicate_points",
34  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  MooseEnum point_not_found_behavior("ERROR WARNING IGNORE", "IGNORE");
40  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  params.addParam<bool>("allow_moving_sources",
48  false,
49  "If true, allow Dirac sources to move, even if the mesh does not move, "
50  "during the simulation.");
51 
52  params.addParamNamesToGroup("use_displaced_mesh drop_duplicate_points", "Advanced");
53  params.declareControllable("enable");
54  params.registerSystemAttributeName("DiracKernel");
55 
56  return params;
57 }
58 
60  : ResidualObject(parameters),
64  BlockRestrictable(this),
65  _current_elem(_assembly.elem()),
66  _coord_sys(_assembly.coordSystem()),
67  _dirac_kernel_info(_subproblem.diracKernelInfo()),
68  _q_point(_assembly.qPoints()),
69  _physical_point(_assembly.physicalPoints()),
70  _qrule(_assembly.qRule()),
71  _JxW(_assembly.JxW()),
72  _drop_duplicate_points(parameters.get<bool>("drop_duplicate_points")),
73  _point_not_found_behavior(
74  parameters.get<MooseEnum>("point_not_found_behavior").getEnum<PointNotFoundBehavior>()),
75  _allow_moving_sources(getParam<bool>("allow_moving_sources"))
76 {
77  // Stateful material properties are not allowed on DiracKernels
79 }
80 
81 Real
83 {
84  return 0;
85 }
86 
87 void
88 DiracKernelBase::addPoint(const Elem * elem, Point p, unsigned /*id*/)
89 {
90  if (!elem || !hasBlocks(elem->subdomain_id()))
91  {
92  std::stringstream msg;
93  msg << "Point " << p << " not found in block(s) " << Moose::stringify(blockIDs(), ", ") << ".";
95  {
97  mooseError(msg.str());
98  break;
100  mooseDoOnce(mooseWarning(msg.str()));
101  break;
103  break;
104  default:
105  mooseError("Internal enum error.");
106  }
107  return;
108  }
109 
110  if (elem->processor_id() != processor_id())
111  return;
112 
113  _dirac_kernel_info.addPoint(elem, p);
115 }
116 
117 const Elem *
118 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  if (id != libMesh::invalid_uint)
126  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  const Elem * elem = _dirac_kernel_info.findPoint(p, _mesh, blockIDs());
132  addPoint(elem, p, id);
133  return elem;
134 }
135 
136 const Elem *
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  const Elem * return_elem = NULL;
143 
144  // May be set if the Elem is found in our cache, otherwise stays as NULL.
145  const Elem * cached_elem = NULL;
146 
147  // OK, the user gave us an ID, let's see if we already have it...
148  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  unsigned int i_found_it = static_cast<unsigned int>(it != _point_cache.end());
152  unsigned int we_found_it = i_found_it;
153  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  if (!we_found_it)
159  {
160  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  if (elem && (elem->processor_id() == processor_id()))
164  {
165  // Add the point to the cache...
166  _point_cache[id] = std::make_pair(elem, p);
167 
168  // ... and to the reverse cache.
169  std::vector<std::pair<Point, unsigned>> & points = _reverse_point_cache[elem];
170  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  addPoint(elem, p, id);
176  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  if (we_found_it && !i_found_it)
185  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  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  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  Point cached_point = (it->second).second;
200 
201  if (cached_point.relative_fuzzy_equals(p))
202  {
203  // Find the cached element associated to this point
204  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  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  updateCaches(cached_elem, NULL, p, id);
215  return_elem = NULL;
216  break; // out of while loop
217  }
218 
219  bool active = cached_elem->active();
220  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  if (active && contains_point)
226  {
227  addPoint(cached_elem, p, id);
228  return_elem = cached_elem;
229  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  else if (!active && contains_point)
235  {
236  // Get the list of active children
237  std::vector<const Elem *> active_children;
238  cached_elem->active_family_tree(active_children);
239 
240  // Linear search through active children for the one that contains p
241  for (unsigned c = 0; c < active_children.size(); ++c)
242  if (active_children[c]->contains_point(p))
243  {
244  updateCaches(cached_elem, active_children[c], p, id);
245  addPoint(active_children[c], p, id);
246  return_elem = active_children[c];
247  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  if (!return_elem)
256  mooseError("Error, Point not found in any of the active children!");
257 
258  break; // out of while loop
259  }
260 
261  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  (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  (!active && !contains_point))
276  {
277  i_need_find_point = true;
278  break; // out of while loop
279  }
280 
281  else
282  mooseError("We'll never get here!");
283  } // if (cached_point.relative_fuzzy_equals(p))
284  else
285  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  unsigned int we_need_find_point = static_cast<unsigned int>(i_need_find_point);
302  comm().max(we_need_find_point);
303 
304  if (we_need_find_point)
305  {
306  // findPoint() is a parallel-only function
307  const Elem * elem = _dirac_kernel_info.findPoint(p, _mesh, blockIDs());
308 
309  updateCaches(cached_elem, elem, p, id);
310  addPoint(elem, p, id);
311  return_elem = elem;
312  }
313 
314  return return_elem;
315 }
316 
317 bool
319 {
320  return _local_dirac_kernel_info.getElements().count(_mesh.elemPtr(elem->id())) != 0;
321 }
322 
323 bool
324 DiracKernelBase::isActiveAtPoint(const Elem * elem, const Point & p)
325 {
326  return _local_dirac_kernel_info.hasPoint(elem, p);
327 }
328 
329 void
331 {
333 
334  // If points can move, we can't trust caches
337 }
338 
339 void
341 {
342  _point_cache.clear();
343  _reverse_point_cache.clear();
344 }
345 
346 void
347 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  _point_cache.erase(id);
352  if (new_elem && (new_elem->processor_id() == processor_id()))
353  _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  reverse_cache_t::iterator it = _reverse_point_cache.find(old_elem);
359  if (it != _reverse_point_cache.end())
360  {
361  reverse_cache_t::mapped_type & points = it->second;
362  {
363  reverse_cache_t::mapped_type::iterator points_it = points.begin(), points_end = points.end();
364 
365  for (; points_it != points_end; ++points_it)
366  {
367  // If the point matches, remove it from the vector of points
368  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  points.erase(points_it);
375  break;
376  }
377  }
378 
379  // If the points vector is now empty, remove old_elem from the reverse cache entirely
380  if (points.empty())
381  _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  if (new_elem && (new_elem->processor_id() == processor_id()))
387  {
388  reverse_cache_t::mapped_type & points = _reverse_point_cache[new_elem];
389  points.push_back(std::make_pair(p, id));
390  }
391 }
392 
393 unsigned
395 {
396  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  if (it == _reverse_point_cache.end())
400  return libMesh::invalid_uint;
401 
402  // Do a linear search in the (hopefully small) vector of Points for this Elem
403  reverse_cache_t::mapped_type & points = it->second;
404 
405  for (const auto & points_it : points)
406  {
407  // If the current_point equals the cached point, return the associated id
408  if (_current_point.relative_fuzzy_equals(points_it.first))
409  return points_it.second;
410  }
411 
412  // If we made it here, we didn't find the cached point, so return invalid_uint
413  return libMesh::invalid_uint;
414 }
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:3108
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:1155
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.
void mooseWarning(Args &&... args) const
Emits a warning prefixed with object name and type.
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:683
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:684
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 mooseError(Args &&... args) const
Emits an error prefixed with object name and type.
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.
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...