www.mooseframework.org
DiracKernel.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 "DiracKernel.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 template <>
22 {
26  params.addRequiredParam<NonlinearVariableName>(
27  "variable", "The name of the variable that this kernel operates on");
28 
29  params.addParam<bool>("use_displaced_mesh",
30  false,
31  "Whether or not this object should use the displaced mesh for computation. "
32  "Note that in the case this is true but no displacements are provided in "
33  "the Mesh block the undisplaced mesh will still be used.");
34 
35  params.addParam<bool>(
36  "drop_duplicate_points",
37  true,
38  "By default points added to a DiracKernel are dropped if a point at the same location"
39  "has been added before. If this option is set to false duplicate points are retained"
40  "and contribute to residual and Jacobian.");
41 
42  params.addParamNamesToGroup("use_displaced_mesh drop_duplicate_points", "Advanced");
43 
44  params.declareControllable("enable");
45  params.registerBase("DiracKernel");
46 
47  return params;
48 }
49 
51  : MooseObject(parameters),
52  SetupInterface(this),
54  MooseVariableInterface<Real>(this,
55  false,
56  "variable",
59  FunctionInterface(this),
60  UserObjectInterface(this),
61  TransientInterface(this),
65  Restartable(this, "DiracKernels"),
66  MeshChangedInterface(parameters),
67  TaggingInterface(this),
68  _subproblem(*getCheckedPointerParam<SubProblem *>("_subproblem")),
69  _sys(*getCheckedPointerParam<SystemBase *>("_sys")),
70  _tid(parameters.get<THREAD_ID>("_tid")),
71  _assembly(_subproblem.assembly(_tid)),
72  _var(*mooseVariable()),
73  _mesh(_subproblem.mesh()),
74  _coord_sys(_assembly.coordSystem()),
75  _dirac_kernel_info(_subproblem.diracKernelInfo()),
76  _current_elem(_var.currentElem()),
77  _q_point(_assembly.qPoints()),
78  _physical_point(_assembly.physicalPoints()),
79  _qrule(_assembly.qRule()),
80  _JxW(_assembly.JxW()),
81  _phi(_assembly.phi(_var)),
82  _grad_phi(_assembly.gradPhi(_var)),
83  _test(_var.phi()),
84  _grad_test(_var.gradPhi()),
85  _u(_var.sln()),
86  _grad_u(_var.gradSln()),
87  _drop_duplicate_points(parameters.get<bool>("drop_duplicate_points"))
88 {
90 
91  // Stateful material properties are not allowed on DiracKernels
93 }
94 
95 void
97 {
99 
100  const std::vector<unsigned int> * multiplicities =
102  unsigned int local_qp = 0;
103  Real multiplicity = 1.0;
104 
105  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
106  {
109  {
111  multiplicity = (*multiplicities)[local_qp++];
112 
113  for (_i = 0; _i < _test.size(); _i++)
114  _local_re(_i) += multiplicity * computeQpResidual();
115  }
116  }
117 
119 }
120 
121 void
123 {
125 
126  const std::vector<unsigned int> * multiplicities =
128  unsigned int local_qp = 0;
129  Real multiplicity = 1.0;
130 
131  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
132  {
135  {
137  multiplicity = (*multiplicities)[local_qp++];
138 
139  for (_i = 0; _i < _test.size(); _i++)
140  for (_j = 0; _j < _phi.size(); _j++)
141  _local_ke(_i, _j) += multiplicity * computeQpJacobian();
142  }
143  }
144 
146 }
147 
148 void
150 {
151  if (jvar == _var.number())
152  {
153  computeJacobian();
154  }
155  else
156  {
158 
159  const std::vector<unsigned int> * multiplicities =
161  unsigned int local_qp = 0;
162  Real multiplicity = 1.0;
163 
164  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
165  {
168  {
170  multiplicity = (*multiplicities)[local_qp++];
171 
172  for (_i = 0; _i < _test.size(); _i++)
173  for (_j = 0; _j < _phi.size(); _j++)
174  _local_ke(_i, _j) += multiplicity * computeQpOffDiagJacobian(jvar);
175  }
176  }
177 
179  }
180 }
181 
182 Real
184 {
185  return 0;
186 }
187 
188 Real
190 {
191  return 0;
192 }
193 
194 void
195 DiracKernel::addPoint(const Elem * elem, Point p, unsigned /*id*/)
196 {
197  if (!elem || (elem->processor_id() != processor_id()))
198  return;
199 
200  _dirac_kernel_info.addPoint(elem, p);
202 }
203 
204 const Elem *
205 DiracKernel::addPoint(Point p, unsigned id)
206 {
207  // Make sure that this method was called with the same id on all
208  // processors. It's an extra communication, though, so let's only
209  // do it in DEBUG mode.
210  libmesh_assert(comm().verify(id));
211 
212  if (id != libMesh::invalid_uint)
213  return addPointWithValidId(p, id);
214 
215  // If id == libMesh::invalid_uint (the default), the user is not
216  // enabling caching when they add Dirac points. So all we can do is
217  // the PointLocator lookup, and call the other addPoint() method.
218  const Elem * elem = _dirac_kernel_info.findPoint(p, _mesh);
219  addPoint(elem, p, id);
220  return elem;
221 }
222 
223 const Elem *
224 DiracKernel::addPointWithValidId(Point p, unsigned id)
225 {
226  // The Elem we'll eventually return. We can't return early on some
227  // processors, because we may need to call parallel_only() functions in
228  // the remainder of this scope.
229  const Elem * return_elem = NULL;
230 
231  // May be set if the Elem is found in our cache, otherwise stays as NULL.
232  const Elem * cached_elem = NULL;
233 
234  // OK, the user gave us an ID, let's see if we already have it...
235  point_cache_t::iterator it = _point_cache.find(id);
236 
237  // Was the point found in a _point_cache on at least one processor?
238  unsigned int i_found_it = static_cast<unsigned int>(it != _point_cache.end());
239  unsigned int we_found_it = i_found_it;
240  comm().max(we_found_it);
241 
242  // If nobody found it in their local caches, it means we need to
243  // do the PointLocator look-up and update the caches. This is
244  // safe, because all processors have the same value of we_found_it.
245  if (!we_found_it)
246  {
247  const Elem * elem = _dirac_kernel_info.findPoint(p, _mesh);
248 
249  // Only add the point to the cache on this processor if the Elem is local
250  if (elem && (elem->processor_id() == processor_id()))
251  {
252  // Add the point to the cache...
253  _point_cache[id] = std::make_pair(elem, p);
254 
255  // ... and to the reverse cache.
256  std::vector<std::pair<Point, unsigned>> & points = _reverse_point_cache[elem];
257  points.push_back(std::make_pair(p, id));
258  }
259 
260  // Call the other addPoint() method. This method ignores non-local
261  // and NULL elements automatically.
262  addPoint(elem, p, id);
263  return_elem = elem;
264  }
265 
266  // If the point was found in a cache, but not my cache, I'm not
267  // responsible for it.
268  //
269  // We can't return early here: then we aren't allowed to call any more
270  // parallel_only() functions in the remainder of this function!
271  if (we_found_it && !i_found_it)
272  return_elem = NULL;
273 
274  // This flag may be set by the processor that cached the Elem because it
275  // needs to call findPoint() (due to moving mesh, etc.). If so, we will
276  // call it at the end of the while loop below.
277  bool i_need_find_point = false;
278 
279  // Now that we only cache local data, some processors may enter
280  // this if statement and some may not. Therefore we can't call
281  // any parallel_only() functions inside this if statement.
282  while (i_found_it)
283  {
284  // We have something cached, now make sure it's actually the same Point.
285  // TODO: we should probably use this same comparison in the DiracKernelInfo code!
286  Point cached_point = (it->second).second;
287 
288  if (cached_point.relative_fuzzy_equals(p))
289  {
290  // Find the cached element associated to this point
291  cached_elem = (it->second).first;
292 
293  // If the cached element's processor ID doesn't match ours, we
294  // are no longer responsible for caching it. This can happen
295  // due to adaptivity...
296  if (cached_elem->processor_id() != processor_id())
297  {
298  // Update the caches, telling them to drop the cached Elem.
299  // Analogously to the rest of the DiracKernel system, we
300  // also return NULL because the Elem is non-local.
301  updateCaches(cached_elem, NULL, p, id);
302  return_elem = NULL;
303  break; // out of while loop
304  }
305 
306  bool active = cached_elem->active();
307  bool contains_point = cached_elem->contains_point(p);
308 
309  // If the cached Elem is active and the point is still
310  // contained in it, call the other addPoint() method and
311  // return its result.
312  if (active && contains_point)
313  {
314  addPoint(cached_elem, p, id);
315  return_elem = cached_elem;
316  break; // out of while loop
317  }
318 
319  // Is the Elem not active (been refined) but still contains the point?
320  // Then search in its active children and update the caches.
321  else if (!active && contains_point)
322  {
323  // Get the list of active children
324  std::vector<const Elem *> active_children;
325  cached_elem->active_family_tree(active_children);
326 
327  // Linear search through active children for the one that contains p
328  for (unsigned c = 0; c < active_children.size(); ++c)
329  if (active_children[c]->contains_point(p))
330  {
331  updateCaches(cached_elem, active_children[c], p, id);
332  addPoint(active_children[c], p, id);
333  return_elem = active_children[c];
334  break; // out of for loop
335  }
336 
337  // If we got here without setting return_elem, it means the Point was
338  // found in the parent element, but not in any of the active
339  // children... this is not possible under normal
340  // circumstances, so something must have gone seriously
341  // wrong!
342  if (!return_elem)
343  mooseError("Error, Point not found in any of the active children!");
344 
345  break; // out of while loop
346  }
347 
348  else if (
349  // Is the Elem active but the point is not contained in it any
350  // longer? (For example, did the Mesh move out from under
351  // it?) Then we fall back to the expensive Point Locator
352  // lookup. TODO: we could try and do something more optimized
353  // like checking if any of the active neighbors contains the
354  // point. Update the caches.
355  (active && !contains_point) ||
356 
357  // The Elem has been refined *and* the Mesh has moved out
358  // from under it, we fall back to doing the expensive Point
359  // Locator lookup. TODO: We could try and look in the
360  // active children of this Elem's neighbors for the Point.
361  // Update the caches.
362  (!active && !contains_point))
363  {
364  i_need_find_point = true;
365  break; // out of while loop
366  }
367 
368  else
369  mooseError("We'll never get here!");
370  } // if (cached_point.relative_fuzzy_equals(p))
371  else
372  mooseError("Cached Dirac point ",
373  cached_point,
374  " already exists with ID: ",
375  id,
376  " and does not match point ",
377  p);
378 
379  // We only want one iteration of this while loop at maximum.
380  i_found_it = false;
381  } // while (i_found_it)
382 
383  // We are back to all processors here because we do not return
384  // early in the code above...
385 
386  // Does we need to call findPoint() on all processors.
387  unsigned int we_need_find_point = static_cast<unsigned int>(i_need_find_point);
388  comm().max(we_need_find_point);
389 
390  if (we_need_find_point)
391  {
392  // findPoint() is a parallel-only function
393  const Elem * elem = _dirac_kernel_info.findPoint(p, _mesh);
394 
395  updateCaches(cached_elem, elem, p, id);
396  addPoint(elem, p, id);
397  return_elem = elem;
398  }
399 
400  return return_elem;
401 }
402 
403 unsigned
405 {
406  reverse_cache_t::iterator it = _reverse_point_cache.find(_current_elem);
407 
408  // If the current Elem is not in the cache, return invalid_uint
409  if (it == _reverse_point_cache.end())
410  return libMesh::invalid_uint;
411 
412  // Do a linear search in the (hopefully small) vector of Points for this Elem
413  reverse_cache_t::mapped_type & points = it->second;
414 
415  for (const auto & points_it : points)
416  {
417  // If the current_point equals the cached point, return the associated id
418  if (_current_point.relative_fuzzy_equals(points_it.first))
419  return points_it.second;
420  }
421 
422  // If we made it here, we didn't find the cached point, so return invalid_uint
423  return libMesh::invalid_uint;
424 }
425 
426 bool
428 {
429  return _local_dirac_kernel_info.getElements().count(_mesh.elemPtr(elem->id())) != 0;
430 }
431 
432 bool
433 DiracKernel::isActiveAtPoint(const Elem * elem, const Point & p)
434 {
435  return _local_dirac_kernel_info.hasPoint(elem, p);
436 }
437 
438 void
440 {
442 }
443 
444 void
446 {
447  _point_cache.clear();
448  _reverse_point_cache.clear();
449 }
450 
453 {
454  return _var;
455 }
456 
457 SubProblem &
459 {
460  return _subproblem;
461 }
462 
463 void
464 DiracKernel::updateCaches(const Elem * old_elem, const Elem * new_elem, Point p, unsigned id)
465 {
466  // Update the point cache. Remove old cached data, only cache
467  // new_elem if it is non-NULL and local.
468  _point_cache.erase(id);
469  if (new_elem && (new_elem->processor_id() == processor_id()))
470  _point_cache[id] = std::make_pair(new_elem, p);
471 
472  // Update the reverse cache
473  //
474  // First, remove the Point from the old_elem's vector
475  reverse_cache_t::iterator it = _reverse_point_cache.find(old_elem);
476  if (it != _reverse_point_cache.end())
477  {
478  reverse_cache_t::mapped_type & points = it->second;
479  {
480  reverse_cache_t::mapped_type::iterator points_it = points.begin(), points_end = points.end();
481 
482  for (; points_it != points_end; ++points_it)
483  {
484  // If the point matches, remove it from the vector of points
485  if (p.relative_fuzzy_equals(points_it->first))
486  {
487  // Vector erasure. It can be slow but these vectors are
488  // generally very short. It also invalidates existing
489  // iterators, so we need to break out of the loop and
490  // not use them any more.
491  points.erase(points_it);
492  break;
493  }
494  }
495 
496  // If the points vector is now empty, remove old_elem from the reverse cache entirely
497  if (points.empty())
498  _reverse_point_cache.erase(old_elem);
499  }
500  }
501 
502  // Next, if new_elem is not NULL and local, add the point to the new_elem's vector
503  if (new_elem && (new_elem->processor_id() == processor_id()))
504  {
505  reverse_cache_t::mapped_type & points = _reverse_point_cache[new_elem];
506  points.push_back(std::make_pair(p, id));
507  }
508 }
InputParameters validParams< MaterialPropertyInterface >()
VarFieldType
Definition: MooseTypes.h:488
unsigned int _qp
Quadrature point index.
Definition: DiracKernel.h:184
A class for creating restricted objects.
Definition: Restartable.h:29
unsigned currentPointCachedID()
Returns the user-assigned ID of the current Dirac point if it exits, and libMesh::invalid_uint otherw...
Definition: DiracKernel.C:404
reverse_cache_t _reverse_point_cache
Definition: DiracKernel.h:229
virtual Real computeQpOffDiagJacobian(unsigned int jvar)
This gets called by computeOffDiagJacobian() at each quadrature point.
Definition: DiracKernel.C:189
void accumulateTaggedLocalResidual()
Local residual blocks will be appended by adding the current local kernel residual.
virtual Real computeQpResidual()=0
This is the virtual that derived classes should override for computing the residual.
virtual Elem * elemPtr(const dof_id_type i)
Definition: MooseMesh.C:2267
SubProblem & subProblem()
Return a reference to the subproblem.
Definition: DiracKernel.C:458
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...
Definition: DiracKernel.C:195
MooseMesh & _mesh
Mesh this kernels acts on.
Definition: DiracKernel.h:165
unsigned int number() const
Get variable number coming from libMesh.
const bool _drop_duplicate_points
drop duplicate points or consider them in residual and Jacobian
Definition: DiracKernel.h:217
const Elem * addPointWithValidId(Point p, unsigned id)
A helper function for addPoint(Point, id) for when id != invalid_uint.
Definition: DiracKernel.C:224
virtual void computeResidual()
Computes the residual for the current element.
Definition: DiracKernel.C:96
bool isActiveAtPoint(const Elem *elem, const Point &p)
Whether or not this DiracKernel has something to distribute at this Point.
Definition: DiracKernel.C:433
MooseVariable & _var
Variable this kernel acts on.
Definition: DiracKernel.h:162
point_cache_t _point_cache
Definition: DiracKernel.h:223
DiracKernel(const InputParameters &parameters)
Definition: DiracKernel.C:50
The main MOOSE class responsible for handling user-defined parameters in almost every MOOSE system...
Assembly & _assembly
Definition: DiracKernel.h:159
void mooseError(Args &&... args) const
Definition: MooseObject.h:147
const Elem *const & _current_elem
Definition: DiracKernel.h:181
Base class for a system (of equations)
Definition: SystemBase.h:92
DiracKernelInfo & _dirac_kernel_info
Place for storing Point/Elem information shared across all DiracKernel objects.
Definition: DiracKernel.h:172
DenseMatrix< Number > _local_ke
Holds residual entries as they are accumulated by this Kernel.
bool hasPointsOnElem(const Elem *elem)
Whether or not this DiracKernel has something to distribute on this element.
Definition: DiracKernel.C:427
void addRequiredParam(const std::string &name, const std::string &doc_string)
This method adds a parameter and documentation string to the InputParameters object that will be extr...
void addMooseVariableDependency(MooseVariableFEBase *var)
Call this function to add the passed in MooseVariableFEBase as a variable that this object depends on...
virtual void computeOffDiagJacobian(unsigned int jvar)
Computes the off-diagonal Jacobian for variable jvar.
Definition: DiracKernel.C:149
virtual Real computeQpJacobian()
This is the virtual that derived classes should override for computing the Jacobian.
Definition: DiracKernel.C:183
MooseVariableFE< Real > * mooseVariable() const
Get the variable that this object is using.
Interface for objects that needs transient capabilities.
unsigned int _j
Definition: DiracKernel.h:195
unsigned int size() const
The number of elements that can currently be stored in the array.
Definition: MooseArray.h:259
MultiPointMap & getPoints()
Returns a writeable reference to the _points container.
Interface for notifications that the mesh has changed.
Every object that can be built by the factory should be derived from this class.
Definition: MooseObject.h:42
VarKindType
Framework-wide stuff.
Definition: MooseTypes.h:481
DiracKernelInfo _local_dirac_kernel_info
Place for storing Point/Elem information only for this DiracKernel.
Definition: DiracKernel.h:175
const MooseArray< Point > & _physical_point
Physical points.
Definition: DiracKernel.h:188
const Elem * findPoint(Point p, const MooseMesh &mesh)
Used by client DiracKernel classes to determine the Elem in which the Point p resides.
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...
Definition: DiracKernel.C:464
void addPoint(const Elem *elem, Point p)
Adds a point source.
const QBase *const & _qrule
Quadrature rule.
Definition: DiracKernel.h:190
Interface for objects that need to use UserObjects.
void clearPoints()
Remove all of the current points and elements.
void accumulateTaggedLocalMatrix()
Local Jacobian blocks will be appended by adding the current local kernel Jacobian.
bool hasPoint(const Elem *elem, Point p)
Return true if we have Point &#39;p&#39; in Element &#39;elem&#39;.
MooseVariable & variable()
The variable number that this kernel operates on.
Definition: DiracKernel.C:452
InputParameters validParams< MooseObject >()
Definition: MooseObject.C:25
Point _current_point
The current point.
Definition: DiracKernel.h:178
const VariableTestValue & _test
Values of test functions at QPs.
Definition: DiracKernel.h:207
unsigned int _i
i-th, j-th index for enumerating shape and test functions
Definition: DiracKernel.h:195
virtual void meshChanged() override
Clear point cache when the mesh changes, so that element coarsening, element deletion, and distributed mesh repartitioning don&#39;t leave this with an invalid cache.
Definition: DiracKernel.C:445
void statefulPropertiesAllowed(bool)
Derived classes can declare whether or not they work with stateful material properties.
const std::set< SubdomainID > EMPTY_BLOCK_IDS
Definition: MooseTypes.h:451
Generic class for solving transient nonlinear problems.
Definition: SubProblem.h:59
An interface for accessing Materials.
const std::set< BoundaryID > EMPTY_BOUNDARY_IDS
Definition: MooseTypes.h:452
InputParameters validParams< TaggingInterface >()
const VariablePhiValue & _phi
Values of shape functions at QPs.
Definition: DiracKernel.h:200
DenseVector< Number > _local_re
Holds residual entries as they are accumulated by this Kernel.
Intermediate base class that ties together all the interfaces for getting MooseVariableFEBases with t...
Interface for objects that need to get values of MooseVariables.
InputParameters validParams< DiracKernel >()
Definition: DiracKernel.C:21
MPI_Comm comm
virtual const OutputTools< Real >::VariableSecond & second()
The second derivative of the variable this object is operating on.
std::set< const Elem * > & getElements()
Returns a writeable reference to the _elements container.
Definition: Moose.h:112
void prepareVectorTag(Assembly &assembly, unsigned int ivar)
Prepare data for computing element residual the according to active tags.
void prepareMatrixTag(Assembly &assembly, unsigned int ivar, unsigned int jvar)
Prepare data for computing element jacobian according to the ative tags.
virtual void computeJacobian()
Computes the jacobian for the current element.
Definition: DiracKernel.C:122
Interface for objects that need to use functions.
SubProblem & _subproblem
Definition: DiracKernel.h:154
Interface class for classes which interact with Postprocessors.
void clearPoints()
Remove all of the current points and elements.
Definition: DiracKernel.C:439
unsigned int THREAD_ID
Definition: MooseTypes.h:161