libMesh
augment_sparsity_on_interface.C
Go to the documentation of this file.
1 // local includes
3 
4 // libMesh includes
5 #include "libmesh/elem.h"
6 
7 using namespace libMesh;
8 
10  boundary_id_type crack_boundary_lower,
11  boundary_id_type crack_boundary_upper) :
12  _mesh(mesh),
13  _crack_boundary_lower(crack_boundary_lower),
14  _crack_boundary_upper(crack_boundary_upper),
15  _initialized(false)
16 {
17  this->mesh_reinit();
18 }
19 
21 {
23  return _lower_to_upper;
24 }
25 
27 {
28  this->_initialized = true;
29 
30  // Loop over all elements (not just local elements) to make sure we find
31  // "neighbor" elements on opposite sides of the crack.
32 
33  // Map from (elem, side) to centroid
34  std::map<std::pair<const Elem *, unsigned char>, Point> lower_centroids;
35  std::map<std::pair<const Elem *, unsigned char>, Point> upper_centroids;
36 
37  for (const auto & elem : _mesh.active_element_ptr_range())
38  for (auto side : elem->side_index_range())
39  if (elem->neighbor_ptr(side) == libmesh_nullptr)
40  {
42  {
43  UniquePtr<const Elem> side_elem = elem->build_side_ptr(side);
44 
45  lower_centroids[std::make_pair(elem, side)] = side_elem->centroid();
46  }
47 
49  {
50  UniquePtr<const Elem> side_elem = elem->build_side_ptr(side);
51 
52  upper_centroids[std::make_pair(elem, side)] = side_elem->centroid();
53  }
54  }
55 
56  // If we're doing a reinit on a distributed mesh then we may not see
57  // all the centroids, or even a matching number of centroids.
58  // std::size_t n_lower_centroids = lower_centroids.size();
59  // std::size_t n_upper_centroids = upper_centroids.size();
60  // libmesh_assert(n_lower_centroids == n_upper_centroids);
61 
62  // Clear _lower_to_upper. This map will be used for matrix assembly later on.
63  _lower_to_upper.clear();
64 
65  // Clear _upper_to_lower. This map will be used for element ghosting
66  // on distributed meshes, communication send_list construction in
67  // parallel, and sparsity calculations
68  _upper_to_lower.clear();
69 
70  // We do an N^2 search to find elements with matching centroids. This could be optimized,
71  // e.g. by first sorting the centroids based on their (x,y,z) location.
72  {
73  std::map<std::pair<const Elem *, unsigned char>, Point>::iterator it = lower_centroids.begin();
74  std::map<std::pair<const Elem *, unsigned char>, Point>::iterator it_end = lower_centroids.end();
75  for ( ; it != it_end; ++it)
76  {
77  Point lower_centroid = it->second;
78 
79  // find closest centroid in upper_centroids
80  Real min_distance = std::numeric_limits<Real>::max();
81 
82  std::map<std::pair<const Elem *, unsigned char>, Point>::iterator inner_it = upper_centroids.begin();
83  std::map<std::pair<const Elem *, unsigned char>, Point>::iterator inner_it_end = upper_centroids.end();
84 
85  for ( ; inner_it != inner_it_end; ++inner_it)
86  {
87  Point upper_centroid = inner_it->second;
88 
89  Real distance = (upper_centroid - lower_centroid).norm();
90  if (distance < min_distance)
91  {
92  min_distance = distance;
93  _lower_to_upper[it->first] = inner_it->first.first;
94  }
95  }
96 
97  // For pairs with local elements, we should have found a
98  // matching pair by now.
99  const Elem * elem = it->first.first;
100  const Elem * neighbor = _lower_to_upper[it->first];
101  if (min_distance < TOLERANCE)
102  {
103  // fill up the inverse map
104  _upper_to_lower[neighbor] = elem;
105  }
106  else
107  {
108  libmesh_assert_not_equal_to(elem->processor_id(), _mesh.processor_id());
109  // This must have a false positive; a remote element would
110  // have been closer.
111  _lower_to_upper.erase(it->first);
112  }
113  }
114 
115  // Let's make sure we didn't miss any upper elements either
116 #ifndef NDEBUG
117  std::map<std::pair<const Elem *, unsigned char>, Point>::iterator inner_it = upper_centroids.begin();
118  std::map<std::pair<const Elem *, unsigned char>, Point>::iterator inner_it_end = upper_centroids.end();
119 
120  for ( ; inner_it != inner_it_end; ++inner_it)
121  {
122  const Elem * neighbor = inner_it->first.first;
123  if (neighbor->processor_id() != _mesh.processor_id())
124  continue;
125  ElementMap::const_iterator utl_it =
126  _upper_to_lower.find(neighbor);
127  libmesh_assert(utl_it != _upper_to_lower.end());
128  }
129 #endif
130  }
131 }
132 
133 void AugmentSparsityOnInterface::operator()
134  (const MeshBase::const_element_iterator & range_begin,
135  const MeshBase::const_element_iterator & range_end,
137  map_type & coupled_elements)
138 {
140 
141  const CouplingMatrix * const null_mat = libmesh_nullptr;
142 
143  for (MeshBase::const_element_iterator elem_it = range_begin;
144  elem_it != range_end; ++elem_it)
145  {
146  const Elem * const elem = *elem_it;
147 
148  if (elem->processor_id() != p)
149  coupled_elements.insert (std::make_pair(elem, null_mat));
150 
151  for (auto side : elem->side_index_range())
152  if (elem->neighbor_ptr(side) == libmesh_nullptr)
153  {
154  ElementSideMap::const_iterator ltu_it =
155  _lower_to_upper.find(std::make_pair(elem, side));
156  if (ltu_it != _lower_to_upper.end())
157  {
158  const Elem * neighbor = ltu_it->second;
159  if (neighbor->processor_id() != p)
160  coupled_elements.insert (std::make_pair(neighbor, null_mat));
161  }
162  }
163 
164  ElementMap::const_iterator utl_it =
165  _upper_to_lower.find(elem);
166  if (utl_it != _upper_to_lower.end())
167  {
168  const Elem * neighbor = utl_it->second;
169  if (neighbor->processor_id() != p)
170  coupled_elements.insert (std::make_pair(neighbor, null_mat));
171  }
172 
173  }
174 }
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
The definition of the const_element_iterator struct.
Definition: mesh_base.h:1494
IntRange< unsigned short > side_index_range() const
Definition: elem.h:2083
unsigned short int side
Definition: xdr_io.C:49
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
MeshBase & mesh
uint8_t processor_id_type
Definition: id_types.h:99
const class libmesh_nullptr_t libmesh_nullptr
bool _initialized
Make sure we&#39;ve been initialized before use.
const ElementSideMap & get_lower_to_upper() const
static const Real TOLERANCE
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
The libMesh namespace provides an interface to certain functionality in the library.
long double max(long double a, double b)
Real distance(const Point &p)
This is the MeshBase class.
Definition: mesh_base.h:68
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
boundary_id_type _crack_boundary_lower
Boundary IDs for the lower and upper faces of the "crack" in the mesh.
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:1967
AugmentSparsityOnInterface(MeshBase &mesh, boundary_id_type crack_boundary_lower, boundary_id_type crack_boundary_upper)
Constructor.
int8_t boundary_id_type
Definition: id_types.h:51
std::unordered_map< const Elem *, const CouplingMatrix * > map_type
What elements do we care about and what variables do we care about on each element?
ElementMap _upper_to_lower
The inverse (ignoring sides) of the above map.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
bool has_boundary_id(const Node *const node, const boundary_id_type id) const
ElementSideMap _lower_to_upper
A map from (lower element ID, side ID) to matching upper element ID.
virtual void mesh_reinit() libmesh_override
Rebuild the cached _lower_to_upper map whenever our Mesh has changed.
std::map< std::pair< const Elem *, unsigned char >, const Elem * > ElementSideMap
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
MeshBase & _mesh
The Mesh we&#39;re calculating on.
processor_id_type processor_id() const
processor_id_type processor_id() const
Definition: dof_object.h:694
This class defines a coupling matrix.