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