Line data Source code
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 243 : MeshAlignment::MeshAlignment(const MooseMesh & mesh, bool require_same_translation)
17 : : MeshAlignmentBase(mesh),
18 243 : _require_same_translation(require_same_translation),
19 243 : _meshes_are_coincident(false),
20 243 : _meshes_have_same_translation(false)
21 : {
22 243 : }
23 :
24 : void
25 155 : MeshAlignment::initialize(
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 155 : _primary_elem_ids = primary_elem_ids;
30 155 : extractFrom1DElements(
31 155 : _primary_elem_ids, _primary_elem_points, _primary_node_ids, _primary_node_points);
32 :
33 155 : extractFromBoundaryInfo(secondary_boundary_info,
34 155 : _secondary_elem_ids,
35 155 : _secondary_side_ids,
36 155 : _secondary_elem_points,
37 155 : _secondary_node_ids,
38 155 : _secondary_node_points);
39 :
40 155 : buildMapping();
41 155 : }
42 :
43 : void
44 76 : MeshAlignment::initialize(
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 76 : extractFromBoundaryInfo(primary_boundary_info,
49 76 : _primary_elem_ids,
50 76 : _primary_side_ids,
51 76 : _primary_elem_points,
52 76 : _primary_node_ids,
53 76 : _primary_node_points);
54 :
55 76 : extractFromBoundaryInfo(secondary_boundary_info,
56 76 : _secondary_elem_ids,
57 76 : _secondary_side_ids,
58 76 : _secondary_elem_points,
59 76 : _secondary_node_ids,
60 76 : _secondary_node_points);
61 :
62 76 : buildMapping();
63 76 : }
64 :
65 : void
66 231 : MeshAlignment::buildMapping()
67 : {
68 231 : _meshes_are_coincident = true;
69 231 : std::vector<unsigned int> primary_elem_pairing_count(_primary_elem_ids.size(), 0);
70 :
71 : // Build the element mapping
72 231 : if (_primary_elem_points.size() > 0 && _secondary_elem_points.size() > 0)
73 : {
74 : // find the primary elements that are nearest to the secondary elements
75 231 : KDTree kd_tree(_primary_elem_points, _mesh.getMaxLeafSize());
76 2934 : for (std::size_t i_secondary = 0; i_secondary < _secondary_elem_points.size(); i_secondary++)
77 : {
78 : unsigned int patch_size = 1;
79 2703 : std::vector<std::size_t> return_index(patch_size);
80 2703 : kd_tree.neighborSearch(_secondary_elem_points[i_secondary], patch_size, return_index);
81 2703 : const std::size_t i_primary = return_index[0];
82 :
83 : // Flip flag if any pair of points are not coincident
84 2703 : if (!_secondary_elem_points[i_secondary].absolute_fuzzy_equals(
85 : _primary_elem_points[i_primary]))
86 844 : _meshes_are_coincident = false;
87 :
88 2703 : const auto primary_elem_id = _primary_elem_ids[i_primary];
89 2703 : const auto secondary_elem_id = _secondary_elem_ids[i_secondary];
90 :
91 2703 : _coupled_elem_ids.insert({primary_elem_id, secondary_elem_id});
92 2703 : _coupled_elem_ids.insert({secondary_elem_id, primary_elem_id});
93 :
94 2703 : primary_elem_pairing_count[i_primary]++;
95 2703 : }
96 : }
97 :
98 : // Build the node mapping
99 : Point reference_translation_vector;
100 231 : _meshes_have_same_translation = true;
101 231 : if (_primary_node_points.size() > 0 && _secondary_node_points.size() > 0)
102 : {
103 : // find the primary nodes that are nearest to the secondary nodes
104 231 : KDTree kd_tree(_primary_node_points, _mesh.getMaxLeafSize());
105 3165 : for (std::size_t i_secondary = 0; i_secondary < _secondary_node_points.size(); i_secondary++)
106 : {
107 : unsigned int patch_size = 1;
108 2934 : std::vector<std::size_t> return_index(patch_size);
109 2934 : kd_tree.neighborSearch(_secondary_node_points[i_secondary], patch_size, return_index);
110 2934 : 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 2934 : if (i_secondary == 0)
116 231 : reference_translation_vector = translation_vector;
117 : else
118 : {
119 2703 : if (!MooseUtils::absoluteFuzzyEqual(
120 2703 : (reference_translation_vector - translation_vector).norm(), 0))
121 66 : _meshes_have_same_translation = false;
122 : }
123 :
124 2934 : const auto primary_node_id = _primary_node_ids[i_primary];
125 2934 : const auto secondary_node_id = _secondary_node_ids[i_secondary];
126 :
127 2934 : _coupled_node_ids.insert({primary_node_id, secondary_node_id});
128 2934 : _coupled_node_ids.insert({secondary_node_id, primary_node_id});
129 2934 : }
130 :
131 454 : if (_meshes_have_same_translation &&
132 223 : MooseUtils::absoluteFuzzyEqual(reference_translation_vector.norm(), 0))
133 133 : _meshes_are_coincident = true;
134 : else
135 98 : _meshes_are_coincident = false;
136 :
137 : // Check if meshes are aligned
138 231 : if (_require_same_translation)
139 155 : _meshes_are_aligned = _meshes_have_same_translation;
140 : else
141 : {
142 : // Else require that all primary boundary elements have exactly one pairing
143 76 : _meshes_are_aligned = true;
144 756 : for (std::size_t i_primary = 0; i_primary < _primary_elem_ids.size(); i_primary++)
145 680 : if (primary_elem_pairing_count[i_primary] != 1)
146 24 : _meshes_are_aligned = false;
147 : }
148 : }
149 231 : }
150 :
151 : void
152 173 : MeshAlignment::buildCoupledElemQpIndexMap(Assembly & assembly)
153 : {
154 : // get local quadrature point maps on each side
155 : std::map<dof_id_type, std::vector<Point>> primary_qp_map =
156 173 : getLocalQuadraturePointMap(assembly, _primary_elem_ids, _primary_side_ids);
157 : std::map<dof_id_type, std::vector<Point>> secondary_qp_map =
158 173 : getLocalQuadraturePointMap(assembly, _secondary_elem_ids, _secondary_side_ids);
159 :
160 : // build mapping
161 2513 : for (const auto & primary_elem_id : _primary_elem_ids)
162 : {
163 2340 : 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 2340 : if (primary_elem && primary_elem->processor_id() == _mesh.processor_id())
167 : {
168 1726 : std::vector<Point> & primary_qps = primary_qp_map[primary_elem_id];
169 :
170 1726 : const dof_id_type secondary_elem_id = getCoupledElemID(primary_elem_id);
171 1726 : 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 1726 : _coupled_elem_qp_indices[primary_elem_id].resize(primary_qps.size());
177 1726 : KDTree kd_tree_qp(secondary_qps, _mesh.getMaxLeafSize());
178 5178 : for (std::size_t i = 0; i < primary_qps.size(); i++)
179 : {
180 : unsigned int patch_size = 1;
181 3452 : std::vector<std::size_t> return_index(patch_size);
182 3452 : kd_tree_qp.neighborSearch(primary_qps[i], patch_size, return_index);
183 3452 : _coupled_elem_qp_indices[primary_elem_id][i] = return_index[0];
184 3452 : }
185 : }
186 : }
187 173 : }
188 :
189 : std::map<dof_id_type, std::vector<Point>>
190 346 : MeshAlignment::getLocalQuadraturePointMap(Assembly & assembly,
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 5026 : for (unsigned int i = 0; i < elem_ids.size(); i++)
196 : {
197 4680 : const auto elem_id = elem_ids[i];
198 4680 : const Elem * elem = _mesh.queryElemPtr(elem_id);
199 4680 : if (elem)
200 : {
201 : assembly.setCurrentSubdomainID(elem->subdomain_id());
202 :
203 : MooseArray<Point> q_points;
204 4349 : if (side_ids.size() > 0)
205 : {
206 2191 : assembly.reinit(elem, side_ids[i]);
207 2191 : q_points = assembly.qPointsFace();
208 : }
209 : else
210 : {
211 2158 : assembly.reinit(elem);
212 2158 : q_points = assembly.qPoints();
213 : }
214 :
215 13047 : for (std::size_t i = 0; i < q_points.size(); i++)
216 8698 : elem_id_to_qps[elem_id].push_back(q_points[i]);
217 : }
218 : }
219 :
220 346 : return elem_id_to_qps;
221 : }
222 :
223 : bool
224 2637 : MeshAlignment::hasCoupledElemID(const dof_id_type & elem_id) const
225 : {
226 2637 : return _coupled_elem_ids.find(elem_id) != _coupled_elem_ids.end();
227 : }
228 :
229 : const dof_id_type &
230 60866 : MeshAlignment::getCoupledElemID(const dof_id_type & elem_id) const
231 : {
232 : mooseAssert(hasCoupledElemID(elem_id), "The element ID has no coupled element.");
233 60866 : return _coupled_elem_ids.find(elem_id)->second;
234 : }
235 :
236 : bool
237 960000 : MeshAlignment::hasCoupledNodeID(const dof_id_type & node_id) const
238 : {
239 960000 : return _coupled_node_ids.find(node_id) != _coupled_node_ids.end();
240 : }
241 :
242 : const dof_id_type &
243 1108690 : MeshAlignment::getCoupledNodeID(const dof_id_type & node_id) const
244 : {
245 : mooseAssert(hasCoupledNodeID(node_id), "The node ID has no coupled node.");
246 1108690 : return _coupled_node_ids.find(node_id)->second;
247 : }
248 :
249 : unsigned int
250 113006 : 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 113006 : return it->second[qp];
257 : }
|
