https://mooseframework.inl.gov
MeshAlignment.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 "MeshAlignment.h"
11 #include "Assembly.h"
12 #include "KDTree.h"
13 
14 #include "libmesh/elem.h"
15 
16 MeshAlignment::MeshAlignment(const MooseMesh & mesh, bool require_same_translation)
18  _require_same_translation(require_same_translation),
19  _meshes_are_coincident(false),
20  _meshes_have_same_translation(false)
21 {
22 }
23 
24 void
26  const std::vector<dof_id_type> & primary_elem_ids,
27  const std::vector<std::tuple<dof_id_type, unsigned short int>> & secondary_boundary_info)
28 {
29  _primary_elem_ids = primary_elem_ids;
32 
33  extractFromBoundaryInfo(secondary_boundary_info,
39 
40  buildMapping();
41 }
42 
43 void
45  const std::vector<std::tuple<dof_id_type, unsigned short int>> & primary_boundary_info,
46  const std::vector<std::tuple<dof_id_type, unsigned short int>> & secondary_boundary_info)
47 {
48  extractFromBoundaryInfo(primary_boundary_info,
54 
55  extractFromBoundaryInfo(secondary_boundary_info,
61 
62  buildMapping();
63 }
64 
65 void
67 {
69  std::vector<unsigned int> primary_elem_pairing_count(_primary_elem_ids.size(), 0);
70 
71  // Build the element mapping
72  if (_primary_elem_points.size() > 0 && _secondary_elem_points.size() > 0)
73  {
74  // find the primary elements that are nearest to the secondary elements
76  for (std::size_t i_secondary = 0; i_secondary < _secondary_elem_points.size(); i_secondary++)
77  {
78  unsigned int patch_size = 1;
79  std::vector<std::size_t> return_index(patch_size);
80  kd_tree.neighborSearch(_secondary_elem_points[i_secondary], patch_size, return_index);
81  const std::size_t i_primary = return_index[0];
82 
83  // Flip flag if any pair of points are not coincident
85  _primary_elem_points[i_primary]))
86  _meshes_are_coincident = false;
87 
88  const auto primary_elem_id = _primary_elem_ids[i_primary];
89  const auto secondary_elem_id = _secondary_elem_ids[i_secondary];
90 
91  _coupled_elem_ids.insert({primary_elem_id, secondary_elem_id});
92  _coupled_elem_ids.insert({secondary_elem_id, primary_elem_id});
93 
94  primary_elem_pairing_count[i_primary]++;
95  }
96  }
97 
98  // Build the node mapping
99  Point reference_translation_vector;
101  if (_primary_node_points.size() > 0 && _secondary_node_points.size() > 0)
102  {
103  // find the primary nodes that are nearest to the secondary nodes
105  for (std::size_t i_secondary = 0; i_secondary < _secondary_node_points.size(); i_secondary++)
106  {
107  unsigned int patch_size = 1;
108  std::vector<std::size_t> return_index(patch_size);
109  kd_tree.neighborSearch(_secondary_node_points[i_secondary], patch_size, return_index);
110  const std::size_t i_primary = return_index[0];
111 
112  const Point translation_vector =
113  _primary_node_points[i_primary] - _secondary_node_points[i_secondary];
114 
115  if (i_secondary == 0)
116  reference_translation_vector = translation_vector;
117  else
118  {
120  (reference_translation_vector - translation_vector).norm(), 0))
122  }
123 
124  const auto primary_node_id = _primary_node_ids[i_primary];
125  const auto secondary_node_id = _secondary_node_ids[i_secondary];
126 
127  _coupled_node_ids.insert({primary_node_id, secondary_node_id});
128  _coupled_node_ids.insert({secondary_node_id, primary_node_id});
129  }
130 
132  MooseUtils::absoluteFuzzyEqual(reference_translation_vector.norm(), 0))
133  _meshes_are_coincident = true;
134  else
135  _meshes_are_coincident = false;
136 
137  // Check if meshes are aligned
140  else
141  {
142  // Else require that all primary boundary elements have exactly one pairing
143  _meshes_are_aligned = true;
144  for (std::size_t i_primary = 0; i_primary < _primary_elem_ids.size(); i_primary++)
145  if (primary_elem_pairing_count[i_primary] != 1)
146  _meshes_are_aligned = false;
147  }
148  }
149 }
150 
151 void
153 {
154  // get local quadrature point maps on each side
155  std::map<dof_id_type, std::vector<Point>> primary_qp_map =
157  std::map<dof_id_type, std::vector<Point>> secondary_qp_map =
159 
160  // build mapping
161  for (const auto & primary_elem_id : _primary_elem_ids)
162  {
163  const Elem * primary_elem = _mesh.queryElemPtr(primary_elem_id);
164  // The PID check is needed to exclude ghost elements, since those don't
165  // necessarily have a coupled element on this processor:
166  if (primary_elem && primary_elem->processor_id() == _mesh.processor_id())
167  {
168  std::vector<Point> & primary_qps = primary_qp_map[primary_elem_id];
169 
170  const dof_id_type secondary_elem_id = getCoupledElemID(primary_elem_id);
171  std::vector<Point> & secondary_qps = secondary_qp_map[secondary_elem_id];
172 
173  mooseAssert(primary_qps.size() == secondary_qps.size(),
174  "The numbers of quadrature points for each element must be the same");
175 
176  _coupled_elem_qp_indices[primary_elem_id].resize(primary_qps.size());
177  KDTree kd_tree_qp(secondary_qps, _mesh.getMaxLeafSize());
178  for (std::size_t i = 0; i < primary_qps.size(); i++)
179  {
180  unsigned int patch_size = 1;
181  std::vector<std::size_t> return_index(patch_size);
182  kd_tree_qp.neighborSearch(primary_qps[i], patch_size, return_index);
183  _coupled_elem_qp_indices[primary_elem_id][i] = return_index[0];
184  }
185  }
186  }
187 }
188 
189 std::map<dof_id_type, std::vector<Point>>
191  const std::vector<dof_id_type> & elem_ids,
192  const std::vector<unsigned short int> & side_ids) const
193 {
194  std::map<dof_id_type, std::vector<Point>> elem_id_to_qps;
195  for (unsigned int i = 0; i < elem_ids.size(); i++)
196  {
197  const auto elem_id = elem_ids[i];
198  const Elem * elem = _mesh.queryElemPtr(elem_id);
199  if (elem)
200  {
201  assembly.setCurrentSubdomainID(elem->subdomain_id());
202 
203  MooseArray<Point> q_points;
204  if (side_ids.size() > 0)
205  {
206  assembly.reinit(elem, side_ids[i]);
207  q_points = assembly.qPointsFace();
208  }
209  else
210  {
211  assembly.reinit(elem);
212  q_points = assembly.qPoints();
213  }
214 
215  for (std::size_t i = 0; i < q_points.size(); i++)
216  elem_id_to_qps[elem_id].push_back(q_points[i]);
217  }
218  }
219 
220  return elem_id_to_qps;
221 }
222 
223 bool
225 {
226  return _coupled_elem_ids.find(elem_id) != _coupled_elem_ids.end();
227 }
228 
229 const dof_id_type &
231 {
232  mooseAssert(hasCoupledElemID(elem_id), "The element ID has no coupled element.");
233  return _coupled_elem_ids.find(elem_id)->second;
234 }
235 
236 bool
238 {
239  return _coupled_node_ids.find(node_id) != _coupled_node_ids.end();
240 }
241 
242 const dof_id_type &
244 {
245  mooseAssert(hasCoupledNodeID(node_id), "The node ID has no coupled node.");
246  return _coupled_node_ids.find(node_id)->second;
247 }
248 
249 unsigned int
250 MeshAlignment::getCoupledElemQpIndex(const dof_id_type & elem_id, const unsigned int & qp) const
251 {
252  auto it = _coupled_elem_qp_indices.find(elem_id);
253  mooseAssert(it != _coupled_elem_qp_indices.end(),
254  "The element ID has no coupled quadrature point indices.");
255  mooseAssert(qp < it->second.size(), "The quadrature index does not exist in the map.");
256  return it->second[qp];
257 }
void extractFrom1DElements(const std::vector< dof_id_type > &elem_ids, std::vector< Point > &elem_points, std::vector< dof_id_type > &node_ids, std::vector< Point > &node_points) const
Extracts mesh information from 1D elements.
std::map< dof_id_type, std::vector< Point > > getLocalQuadraturePointMap(Assembly &assembly, const std::vector< dof_id_type > &elem_ids, const std::vector< unsigned short int > &side_ids) const
Gets the local quadrature point map for the primary or secondary side.
bool _meshes_have_same_translation
Flag that each pairing has the same translation vector.
void buildMapping()
Builds the mapping using the extracted mesh information.
Definition: MeshAlignment.C:66
std::vector< Point > _secondary_elem_points
List of secondary element points.
bool absoluteFuzzyEqual(const T &var1, const T2 &var2, const T3 &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
Builds mapping between two aligned subdomains/boundaries.
std::map< dof_id_type, dof_id_type > _coupled_node_ids
Map of node ID to coupled node ID.
bool hasCoupledNodeID(const dof_id_type &node_id) const
Returns true if the node ID has a coupled node ID.
const MooseMesh & _mesh
Mesh.
bool _meshes_are_aligned
Flag that meshes are aligned.
MeshBase & mesh
void initialize(const std::vector< dof_id_type > &primary_elem_ids, const std::vector< std::tuple< dof_id_type, unsigned short int >> &secondary_boundary_info)
Extracts mesh information and builds the mapping.
Definition: MeshAlignment.C:25
std::vector< dof_id_type > _secondary_elem_ids
List of secondary element IDs.
const MooseArray< Point > & qPoints() const
std::vector< dof_id_type > _secondary_node_ids
List of secondary node IDs.
std::vector< unsigned short int > _secondary_side_ids
List of secondary side IDs (if any)
void reinit(const Elem *elem)
virtual Elem * queryElemPtr(const dof_id_type i)
const dof_id_type & getCoupledNodeID(const dof_id_type &node_id) const
Gets the coupled node ID for a given node ID.
const dof_id_type & getCoupledElemID(const dof_id_type &elem_id) const
Gets the coupled element ID for a given element ID.
unsigned int size() const
const bool _require_same_translation
Flag that alignment requires that pairings have same translation vector.
unsigned int getMaxLeafSize() const
MeshAlignment(const MooseMesh &mesh, bool require_same_translation=false)
Constructor.
Definition: MeshAlignment.C:16
std::vector< dof_id_type > _primary_node_ids
List of primary node IDs.
auto norm(const T &a) -> decltype(std::abs(a))
std::vector< Point > _primary_elem_points
List of primary element points.
void buildCoupledElemQpIndexMap(Assembly &assembly)
Builds the map used for getting the coupled quadrature point index.
std::vector< Point > _primary_node_points
List of primary node points.
bool absolute_fuzzy_equals(const T &var1, const T2 &var2, const Real tol=TOLERANCE *TOLERANCE)
void setCurrentSubdomainID(SubdomainID i)
const MooseArray< Point > & qPointsFace() const
std::map< dof_id_type, std::vector< unsigned int > > _coupled_elem_qp_indices
Map of element ID to vector of coupled quadrature points.
std::vector< unsigned short int > _primary_side_ids
List of primary side IDs (if any)
void extractFromBoundaryInfo(const std::vector< std::tuple< dof_id_type, unsigned short int >> &boundary_info, std::vector< dof_id_type > &elem_ids, std::vector< unsigned short int > &side_ids, std::vector< Point > &side_points, std::vector< dof_id_type > &node_ids, std::vector< Point > &node_points) const
Extracts mesh information from boundary info.
std::vector< dof_id_type > _primary_elem_ids
List of primary element IDs.
std::map< dof_id_type, dof_id_type > _coupled_elem_ids
Map of element ID to coupled element ID.
processor_id_type processor_id() const
unsigned int getCoupledElemQpIndex(const dof_id_type &elem_id, const unsigned int &qp) const
Gets the quadrature point index on the coupled element corresponding to the quadrature point index on...
bool _meshes_are_coincident
Flag that meshes are coincident.
uint8_t dof_id_type
bool hasCoupledElemID(const dof_id_type &elem_id) const
Returns true if the element ID has a coupled element ID.
std::vector< Point > _secondary_node_points
List of secondary node points.