https://mooseframework.inl.gov
DiracKernelInfo.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 #include "DiracKernelInfo.h"
11 #include "MooseMesh.h"
12 
13 // LibMesh
14 #include "libmesh/point_locator_base.h"
15 #include "libmesh/elem.h"
16 #include "libmesh/enum_point_locator_type.h"
17 #include "libmesh/point.h"
18 
19 using namespace libMesh;
20 
22  : _point_locator(), _point_equal_distance_sq(libMesh::TOLERANCE * libMesh::TOLERANCE)
23 {
24 }
25 
27 
28 void
29 DiracKernelInfo::addPoint(const Elem * elem, const Point & p)
30 {
31  _elements.insert(elem);
32 
33  std::pair<std::vector<Point>, std::vector<unsigned int>> & multi_point_list = _points[elem];
34 
35  const unsigned int npoint = multi_point_list.first.size();
36  mooseAssert(npoint == multi_point_list.second.size(),
37  "Different sizes for location and multiplicity data");
38 
39  for (unsigned int i = 0; i < npoint; ++i)
40  if (pointsFuzzyEqual(multi_point_list.first[i], p))
41  {
42  // a point at the same (within a tolerance) location as p exists, increase its multiplicity
43  multi_point_list.second[i]++;
44  return;
45  }
46 
47  // no prior point found at this location, add it with a multiplicity of one
48  multi_point_list.first.push_back(p);
49  multi_point_list.second.push_back(1);
50 }
51 
52 void
54 {
55  _elements.clear();
56  _points.clear();
57 }
58 
59 bool
60 DiracKernelInfo::hasPoint(const Elem * elem, const Point & p)
61 {
62  std::vector<Point> & point_list = _points[elem].first;
63 
64  for (const auto & pt : point_list)
65  if (pointsFuzzyEqual(pt, p))
66  return true;
67 
68  // If we haven't found it, we don't have it.
69  return false;
70 }
71 
72 void
74 {
75  // Note: we could update the PointLocator *every* time we call this
76  // function, but that may introduce an unacceptable overhead in
77  // problems which don't need a PointLocator at all. This issue will
78  // most likely become a moot point when we eventually add a shared
79  // "CachingPointLocator" to MOOSE.
80  // _point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, mesh);
81 
82  // Construct the PointLocator object if *any* processors have Dirac
83  // points. Note: building a PointLocator object is a parallel_only()
84  // function, so this is an all-or-nothing thing.
85  unsigned pl_needs_rebuild = _elements.size();
86  mesh.comm().max(pl_needs_rebuild);
87 
88  if (pl_needs_rebuild)
89  {
90  // PointLocatorBase::build() is a parallel_only function! So we
91  // can't skip building it just becuase our local _elements is
92  // empty, it might be non-empty on some other processor!
93  _point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, mesh);
94 
95  // We may be querying for points which are not in the semilocal
96  // part of a distributed mesh.
97  _point_locator->enable_out_of_mesh_mode();
98  }
99  else
100  {
101  // There are no elements with Dirac points, but we have been
102  // requested to update the PointLocator so we have to assume the
103  // old one is invalid. Therefore we reset it to NULL... however
104  // adding this line causes the code to hang because it triggers
105  // the PointLocator to be rebuilt in a non-parallel-only segment
106  // of the code later... so it's commented out for now even though
107  // it's probably the right behavior.
108  // _point_locator.reset(NULL);
109  }
110 }
111 
112 const Elem *
114  const MooseMesh & mesh,
115  const std::set<SubdomainID> & blocks)
116 {
117  // If the PointLocator has never been created, do so now. NOTE - WE
118  // CAN'T DO THIS if findPoint() is only called on some processors,
119  // PointLocatorBase::build() is a 'parallel_only' method!
120  if (_point_locator.get() == NULL)
121  {
122  _point_locator = PointLocatorBase::build(TREE_LOCAL_ELEMENTS, mesh);
123  _point_locator->enable_out_of_mesh_mode();
124  }
125 
126  // Check that the PointLocator is ready to start locating points.
127  // So far I do not have any tests that trip this...
128  if (_point_locator->initialized() == false)
129  mooseError("Error, PointLocator is not initialized!");
130 
131  // Note: The PointLocator object returns NULL when the Point is not
132  // found within the Mesh. This is not considered to be an error as
133  // far as the DiracKernels are concerned: sometimes the Mesh moves
134  // out from the Dirac point entirely and in that case the Point just
135  // gets "deactivated".
136  const Elem * elem = (*_point_locator)(p, &blocks);
137 
138  // The processors may not agree on which Elem the point is in. This
139  // can happen if a Dirac point lies on the processor boundary, and
140  // two or more neighboring processors think the point is in the Elem
141  // on *their* side.
142  dof_id_type elem_id = elem ? elem->id() : DofObject::invalid_id;
143 
144  // We are going to let the element with the smallest ID "win", all other
145  // procs will return NULL.
146  dof_id_type min_elem_id = elem_id;
147  mesh.comm().min(min_elem_id);
148 
149  return min_elem_id == elem_id ? elem : NULL;
150 }
151 
152 bool
154 {
155  const Real dist_sq = (a - b).norm_sq();
156  return dist_sq < _point_equal_distance_sq;
157 }
bool hasPoint(const Elem *elem, const Point &p)
Return true if we have Point &#39;p&#39; in Element &#39;elem&#39;.
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.
std::set< const Elem * > _elements
The list of elements that need distributions.
char ** blocks
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:302
static constexpr Real TOLERANCE
MeshBase & mesh
const Parallel::Communicator & comm() const
MultiPointMap _points
The list of physical xyz Points that need to be evaluated in each element.
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
void addPoint(const Elem *elem, const Point &p)
Adds a point source.
std::unique_ptr< libMesh::PointLocatorBase > _point_locator
The DiracKernelInfo object manages a PointLocator object which is used by all DiracKernels to find Po...
dof_id_type id() const
void min(const T &r, T &o, Request &req) const
MooseMesh wraps a libMesh::Mesh object and enhances its capabilities by caching additional data and s...
Definition: MooseMesh.h:88
void clearPoints()
Remove all of the current points and elements.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void max(const T &r, T &o, Request &req) const
auto norm_sq(const T &a) -> decltype(std::norm(a))
const Real _point_equal_distance_sq
threshold distance squared below which two points are considered identical
TREE_LOCAL_ELEMENTS
void updatePointLocator(const MooseMesh &mesh)
Called during FEProblemBase::meshChanged() to update the PointLocator object used by the DiracKernels...
uint8_t dof_id_type
virtual ~DiracKernelInfo()
bool pointsFuzzyEqual(const Point &, const Point &)
Check if two points are equal with respect to a tolerance.