Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
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 "KokkosMesh.h"
11 :
12 : #include "Assembly.h"
13 : #include "MooseMesh.h"
14 :
15 : #include "libmesh/elem_side_builder.h"
16 : #include "libmesh/reference_elem.h"
17 :
18 : namespace Moose::Kokkos
19 : {
20 :
21 : void
22 2622 : Mesh::update()
23 : {
24 2622 : initMap();
25 2622 : initElement();
26 :
27 2622 : _initialized = true;
28 :
29 2622 : if (_needs_element_geometry)
30 181 : initElementGeometry();
31 2622 : if (_needs_element_side_geometry)
32 181 : initElementSideGeometry();
33 2622 : }
34 :
35 : void
36 2622 : Mesh::initMap()
37 : {
38 2622 : if (!_maps)
39 2622 : _maps = std::make_shared<MeshMap>();
40 :
41 2622 : if (_elem_id_integer == libMesh::invalid_uint)
42 7866 : _elem_id_integer = _mesh.getMesh().add_elem_integer("kokkos_contiguous_elem_id");
43 :
44 2622 : if (_node_id_integer == libMesh::invalid_uint)
45 7866 : _node_id_integer = _mesh.getMesh().add_node_integer("kokkos_contiguous_node_id");
46 :
47 2622 : _num_local_elems = 0;
48 2622 : _num_ghost_elems = 0;
49 2622 : _num_local_nodes = 0;
50 :
51 2622 : _maps->subdomain_id_mapping.clear();
52 2622 : _maps->boundary_id_mapping.clear();
53 2622 : _maps->elem_type_id_mapping.clear();
54 2622 : _maps->local_nodes.clear();
55 2622 : _maps->ghost_node_id_mapping.clear();
56 2622 : _maps->ghost_elem_id_mapping.clear();
57 2622 : _maps->subdomain_elem_id_ranges.clear();
58 2622 : _maps->subdomain_node_ids.clear();
59 2622 : _maps->boundary_node_ids.clear();
60 :
61 2622 : std::unordered_set<ElemType> elem_types;
62 2622 : std::unordered_set<Node *> ghost_nodes;
63 :
64 5974 : for (const auto subdomain : _mesh.meshSubdomains())
65 : {
66 3352 : _maps->subdomain_id_mapping[subdomain] = _maps->subdomain_id_mapping.size();
67 3352 : _maps->subdomain_node_ids[subdomain];
68 :
69 3352 : dof_id_type begin = libMesh::DofObject::invalid_id;
70 3352 : dof_id_type eid = libMesh::DofObject::invalid_id;
71 :
72 681952 : for (auto elem : _mesh.getMesh().active_local_subdomain_elements_ptr_range(subdomain))
73 : {
74 678600 : elem_types.insert(elem->type());
75 :
76 678600 : eid = _num_local_elems;
77 678600 : elem->set_extra_integer(_elem_id_integer, _num_local_elems++);
78 :
79 678600 : if (begin == libMesh::DofObject::invalid_id)
80 3156 : begin = eid;
81 :
82 3463740 : for (auto & node : elem->node_ref_range())
83 2785140 : if (node.processor_id() != _mesh.processor_id())
84 14870 : ghost_nodes.insert(&node);
85 3352 : }
86 :
87 3352 : _maps->subdomain_elem_id_ranges[subdomain] = eid != libMesh::DofObject::invalid_id
88 3616 : ? std::make_pair(begin, eid + 1)
89 1906 : : std::make_pair(begin, eid);
90 : }
91 :
92 12869 : for (const auto boundary : _mesh.meshBoundaryIds())
93 : {
94 10247 : _maps->boundary_id_mapping[boundary] = _maps->boundary_id_mapping.size();
95 10247 : _maps->boundary_node_ids[boundary];
96 : }
97 :
98 : // Enumerate one off-process neighbor layer for element-side geometry data.
99 2622 : if (_needs_element_side_geometry)
100 21554 : for (const auto elem : _mesh.getMesh().active_local_element_ptr_range())
101 102487 : for (const auto side : elem->side_index_range())
102 : {
103 81114 : const auto neighbor = elem->neighbor_ptr(side);
104 79330 : if (neighbor && neighbor != libMesh::remote_elem &&
105 79408 : neighbor->processor_id() != _mesh.processor_id() &&
106 16 : !_maps->ghost_elem_id_mapping.count(neighbor))
107 : {
108 32 : _maps->ghost_elem_id_mapping[neighbor] = _num_local_elems + _num_ghost_elems++;
109 32 : elem_types.insert(neighbor->type());
110 : }
111 181 : }
112 :
113 : // Setup node maps
114 :
115 5236 : for (const auto type : elem_types)
116 2614 : _maps->elem_type_id_mapping[type] = _maps->elem_type_id_mapping.size();
117 :
118 4036 : _maps->local_nodes.insert(
119 2828 : _maps->local_nodes.end(), _mesh.localNodesBegin(), _mesh.localNodesEnd());
120 2622 : _maps->local_nodes.insert(_maps->local_nodes.end(), ghost_nodes.begin(), ghost_nodes.end());
121 :
122 776312 : for (auto node : _maps->local_nodes)
123 : {
124 773690 : if (node->processor_id() == _mesh.processor_id())
125 : {
126 765666 : node->set_extra_integer(_node_id_integer, _num_local_nodes);
127 :
128 1617937 : for (auto subdomain : _mesh.getNodeBlockIds(*node))
129 852271 : _maps->subdomain_node_ids[subdomain].push_back(_num_local_nodes);
130 : }
131 : else
132 8024 : _maps->ghost_node_id_mapping[node] = _num_local_nodes;
133 :
134 773690 : ++_num_local_nodes;
135 : }
136 :
137 134069 : for (const auto bnd_node : as_range(_mesh.bndNodesBegin(), _mesh.bndNodesEnd()))
138 131447 : if (bnd_node->_node->processor_id() == _mesh.processor_id())
139 112217 : _maps->boundary_node_ids[bnd_node->_bnd_id].insert(getContiguousNodeID(bnd_node->_node));
140 2622 : }
141 :
142 : void
143 2622 : Mesh::initElement()
144 : {
145 : // Cache reference element data
146 :
147 2622 : const auto num_elem_types = getNumLocalElementTypes();
148 :
149 2622 : _num_sides.create(num_elem_types);
150 2622 : _num_nodes.create(num_elem_types);
151 2622 : _num_side_nodes.create(num_elem_types);
152 2622 : _local_side_node.create(num_elem_types);
153 :
154 5236 : for (auto & [elem_type, elem_type_id] : getElementTypeMap())
155 : {
156 2614 : auto elem = &libMesh::ReferenceElem::get(elem_type);
157 :
158 2614 : _num_sides[elem_type_id] = elem->n_sides();
159 2614 : _num_nodes[elem_type_id] = elem->n_nodes();
160 :
161 2614 : _num_side_nodes[elem_type_id].create(elem->n_sides());
162 2614 : _local_side_node[elem_type_id].create(elem->n_nodes(), elem->n_sides());
163 :
164 12472 : for (const auto side : elem->side_index_range())
165 : {
166 9858 : const auto side_ptr = elem->build_side_ptr(side);
167 :
168 9858 : _num_side_nodes[elem_type_id][side] = side_ptr->n_nodes();
169 :
170 30334 : for (const auto node : side_ptr->node_index_range())
171 20476 : _local_side_node[elem_type_id](node, side) = elem->local_side_node(side, node);
172 9858 : }
173 :
174 2614 : _num_side_nodes[elem_type_id].moveToDevice();
175 2614 : _local_side_node[elem_type_id].moveToDevice();
176 : }
177 :
178 2622 : _num_sides.moveToDevice();
179 2622 : _num_nodes.moveToDevice();
180 2622 : _num_side_nodes.copyToDevice();
181 2622 : _local_side_node.copyToDevice();
182 :
183 : // Cache element data
184 :
185 1208 : const auto num_local_plus_one_neighbor_layer_elems =
186 1414 : getNumLocalAndPossiblyOneNeighborLayerGhostElements();
187 2622 : const auto num_elems = getNumLocalElements();
188 2622 : const auto num_subdomains = getNumSubdomains();
189 :
190 2622 : _elem_info.create(num_local_plus_one_neighbor_layer_elems);
191 2622 : _elem_neighbor.create(_mesh.getMaxSidesPerElem(), num_elems);
192 2622 : _extra_elem_ids.create(num_elems, _mesh.getMesh().n_elem_integers());
193 2622 : _starting_elem_id.create(num_subdomains);
194 :
195 2622 : _elem_neighbor = libMesh::DofObject::invalid_id;
196 2622 : _extra_elem_ids = libMesh::DofObject::invalid_id;
197 :
198 5974 : for (const auto subdomain : _mesh.meshSubdomains())
199 : {
200 3352 : const auto sid = getContiguousSubdomainID(subdomain);
201 :
202 1360552 : for (const auto elem : _mesh.getMesh().active_local_subdomain_elements_ptr_range(subdomain))
203 : {
204 678600 : const auto elem_type = getElementTypeID(elem);
205 678600 : const auto eid = getContiguousElementID(elem);
206 :
207 678600 : _elem_info[eid].type = elem_type;
208 678600 : _elem_info[eid].id = eid;
209 678600 : _elem_info[eid].subdomain = sid;
210 :
211 2365716 : for (const auto i : make_range(_mesh.getMesh().n_elem_integers()))
212 1687116 : _extra_elem_ids(eid, i) = elem->get_extra_integer(i);
213 3352 : }
214 :
215 3352 : _starting_elem_id[sid] = *getSubdomainContiguousElementIDRange(subdomain).begin();
216 : }
217 :
218 : // Cache ElementInfo for one-layer ghost neighbors so side-based computations can identify
219 : // off-process neighbor element subdomains.
220 2654 : for (const auto & [ghost_elem, ghost_eid] : _maps->ghost_elem_id_mapping)
221 : {
222 32 : _elem_info[ghost_eid].type = getElementTypeID(ghost_elem);
223 32 : _elem_info[ghost_eid].id = ghost_eid;
224 32 : _elem_info[ghost_eid].subdomain = getContiguousSubdomainID(ghost_elem->subdomain_id());
225 : }
226 :
227 2622 : _elem_info.moveToDevice();
228 2622 : _extra_elem_ids.moveToDevice();
229 2622 : _starting_elem_id.moveToDevice();
230 :
231 : // Cache node data
232 :
233 2622 : const auto num_nodes = getNumLocalNodes();
234 :
235 2622 : _points.create(num_nodes);
236 2622 : _nodes.create(_mesh.getMaxNodesPerElem(), num_local_plus_one_neighbor_layer_elems);
237 2622 : _boundary_nodes.create(_mesh.meshBoundaryIds().size());
238 :
239 776312 : for (const auto node : getLocalNodes())
240 773690 : _points[getContiguousNodeID(node)] = *node;
241 :
242 1359822 : for (const auto elem : _mesh.getMesh().active_local_element_ptr_range())
243 : {
244 678600 : const auto eid = getContiguousElementID(elem);
245 :
246 3463740 : for (const auto node : elem->node_index_range())
247 2785140 : _nodes(node, eid) = getContiguousNodeID(elem->node_ptr(node));
248 2622 : }
249 :
250 12869 : for (const auto boundary : _mesh.meshBoundaryIds())
251 10247 : _boundary_nodes[getContiguousBoundaryID(boundary)].copySet<false, true>(
252 5563 : _maps->boundary_node_ids[boundary]);
253 :
254 2622 : _points.moveToDevice();
255 2622 : _nodes.moveToDevice();
256 2622 : _boundary_nodes.copyToDevice();
257 2622 : }
258 :
259 : ContiguousSubdomainID
260 1092290 : Mesh::getContiguousSubdomainID(const SubdomainID subdomain) const
261 : {
262 1092290 : return libmesh_map_find(_maps->subdomain_id_mapping, subdomain);
263 : }
264 :
265 : ContiguousBoundaryID
266 10487 : Mesh::getContiguousBoundaryID(const BoundaryID boundary) const
267 : {
268 10487 : return libmesh_map_find(_maps->boundary_id_mapping, boundary);
269 : }
270 :
271 : unsigned int
272 678632 : Mesh::getElementTypeID(const Elem * elem) const
273 : {
274 : mooseAssert(elem, "Element pointer is null");
275 :
276 678632 : return libmesh_map_find(_maps->elem_type_id_mapping, elem->type());
277 : }
278 :
279 : ContiguousElementID
280 4131192 : Mesh::getContiguousElementID(const Elem * elem) const
281 : {
282 : mooseAssert(elem, "Element pointer is null");
283 :
284 4131192 : if (elem->processor_id() == _mesh.processor_id())
285 4131160 : return elem->get_extra_integer(_elem_id_integer);
286 :
287 32 : return libmesh_map_find(_maps->ghost_elem_id_mapping, elem);
288 : }
289 :
290 : ContiguousNodeID
291 5209386 : Mesh::getContiguousNodeID(const Node * node) const
292 : {
293 : mooseAssert(node, "Node pointer is null");
294 :
295 5209386 : if (node->processor_id() == _mesh.processor_id())
296 5186492 : return node->get_extra_integer(_node_id_integer);
297 : else
298 22894 : return libmesh_map_find(_maps->ghost_node_id_mapping, node);
299 : }
300 :
301 : void
302 181 : Mesh::initElementGeometry()
303 : {
304 181 : const auto num_local_elems = getNumLocalElements();
305 :
306 181 : _elem_volume.create(num_local_elems);
307 181 : _elem_centroid.create(num_local_elems);
308 :
309 181 : _elem_volume = 0;
310 181 : _elem_centroid = Real3();
311 :
312 42927 : for (const auto elem : _mesh.getMesh().active_local_element_ptr_range())
313 : {
314 21373 : const auto eid = getContiguousElementID(elem);
315 21373 : const auto elem_centroid = elem->vertex_average();
316 : Real coord_factor;
317 21373 : ::coordTransformFactor(_mesh, elem->subdomain_id(), elem_centroid, coord_factor);
318 :
319 21373 : _elem_volume[eid] = elem->volume() * coord_factor;
320 21373 : _elem_centroid[eid] = elem_centroid;
321 181 : }
322 :
323 181 : _elem_volume.moveToDevice();
324 181 : _elem_centroid.moveToDevice();
325 :
326 181 : _element_geometry_initialized = true;
327 181 : }
328 :
329 : void
330 181 : Mesh::initElementSideGeometry()
331 : {
332 : mooseAssert(_element_geometry_initialized,
333 : "initElementGeometry() must be called before initElementSideGeometry().");
334 :
335 181 : const auto num_elems = getNumLocalElements();
336 181 : const auto max_sides = _mesh.getMaxSidesPerElem();
337 :
338 181 : _side_area.create(max_sides, num_elems);
339 181 : _side_centroid.create(max_sides, num_elems);
340 181 : _side_normal.create(max_sides, num_elems);
341 181 : _elem_centroid_to_side_centroid.create(max_sides, num_elems);
342 181 : _elem_centroid_to_side_centroid_distance.create(max_sides, num_elems);
343 181 : _elem_centroid_to_neighbor_centroid.create(max_sides, num_elems);
344 181 : _elem_centroid_to_neighbor_centroid_distance.create(max_sides, num_elems);
345 181 : _side_boundary_id.create(max_sides, num_elems);
346 :
347 181 : _side_area = 0;
348 181 : _side_centroid = Real3();
349 181 : _side_normal = Real3();
350 181 : _elem_centroid_to_side_centroid = Real3();
351 181 : _elem_centroid_to_side_centroid_distance = 0;
352 181 : _elem_centroid_to_neighbor_centroid = Real3();
353 181 : _elem_centroid_to_neighbor_centroid_distance = 0;
354 181 : _side_boundary_id = Moose::INVALID_BOUNDARY_ID;
355 :
356 181 : libMesh::ElemSideBuilder side_builder;
357 :
358 42927 : for (const auto elem : _mesh.getMesh().active_local_element_ptr_range())
359 : {
360 21373 : const auto eid = getContiguousElementID(elem);
361 21373 : const auto elem_centroid = elem->vertex_average();
362 :
363 102487 : for (const auto side : elem->side_index_range())
364 : {
365 81114 : auto & face = side_builder(*elem, side);
366 81114 : const auto face_centroid = face.vertex_average();
367 81114 : const auto * const neighbor = elem->neighbor_ptr(side);
368 : Real coord_factor;
369 : mooseAssert(!neighbor || neighbor != libMesh::remote_elem,
370 : "First layer neighbors should be ghosted");
371 81114 : ::coordTransformFactor(_mesh,
372 40588 : elem->subdomain_id(),
373 : face_centroid,
374 : coord_factor,
375 38804 : neighbor ? neighbor->subdomain_id()
376 : : libMesh::Elem::invalid_subdomain_id);
377 :
378 81114 : ContiguousElementID neighbor_eid = libMesh::DofObject::invalid_id;
379 81114 : if (neighbor)
380 : {
381 77562 : neighbor_eid = getContiguousElementID(neighbor);
382 : mooseAssert(
383 : neighbor_eid != libMesh::DofObject::invalid_id,
384 : "We should have determined a contiguous element ID for our local element neighbors");
385 77562 : _elem_neighbor(side, eid) = neighbor_eid;
386 : }
387 :
388 81114 : const auto elem_centroid_to_side_centroid = face_centroid - elem_centroid;
389 40526 : const auto elem_centroid_to_neighbor_centroid =
390 40588 : neighbor ? neighbor->vertex_average() - elem_centroid : Real3();
391 :
392 81114 : _side_area(side, eid) = face.volume() * coord_factor;
393 81114 : _side_centroid(side, eid) = face_centroid;
394 81114 : _side_normal(side, eid) = elem->side_vertex_average_normal(side);
395 81114 : _elem_centroid_to_side_centroid(side, eid) = elem_centroid_to_side_centroid;
396 81114 : _elem_centroid_to_side_centroid_distance(side, eid) = elem_centroid_to_side_centroid.norm();
397 81114 : _elem_centroid_to_neighbor_centroid(side, eid) = elem_centroid_to_neighbor_centroid;
398 81114 : _elem_centroid_to_neighbor_centroid_distance(side, eid) =
399 40588 : neighbor ? elem_centroid_to_neighbor_centroid.norm() : Real(0);
400 :
401 81114 : const auto boundary_ids = _mesh.getBoundaryIDs(elem, side);
402 81114 : if (boundary_ids.size() > 1)
403 0 : mooseError("Kokkos element-side geometry does not support multiple boundary IDs on a "
404 : "single side. "
405 : "Element ",
406 : eid,
407 : " side ",
408 : side,
409 : " has ",
410 0 : boundary_ids.size(),
411 : " boundary IDs.");
412 81114 : if (!boundary_ids.empty())
413 4048 : _side_boundary_id(side, eid) = boundary_ids.front();
414 81114 : }
415 181 : }
416 :
417 181 : _elem_neighbor.moveToDevice();
418 181 : _side_area.moveToDevice();
419 181 : _side_centroid.moveToDevice();
420 181 : _side_normal.moveToDevice();
421 181 : _elem_centroid_to_side_centroid.moveToDevice();
422 181 : _elem_centroid_to_side_centroid_distance.moveToDevice();
423 181 : _elem_centroid_to_neighbor_centroid.moveToDevice();
424 181 : _elem_centroid_to_neighbor_centroid_distance.moveToDevice();
425 181 : _side_boundary_id.moveToDevice();
426 :
427 181 : _element_side_geometry_initialized = true;
428 181 : }
429 :
430 : } // namespace Moose::Kokkos
|
