Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #include "mesh.hpp"
5 :
6 : #include "fem/coefficient.hpp"
7 : #include "fem/fespace.hpp"
8 : #include "fem/libceed/integrator.hpp"
9 :
10 : namespace palace
11 : {
12 :
13 : namespace
14 : {
15 :
16 11627 : const auto &GetParentMesh(const mfem::ParMesh &mesh)
17 : {
18 : // Get the parent mesh if the mesh is a boundary submesh (no submesh of submesh
19 : // capabilities, for now).
20 11627 : const auto *submesh = dynamic_cast<const mfem::ParSubMesh *>(&mesh);
21 11627 : if (submesh && submesh->GetFrom() == mfem::SubMesh::From::Boundary)
22 : {
23 0 : return *submesh->GetParent();
24 : }
25 : return mesh;
26 : }
27 :
28 : auto &GetParentMesh(mfem::ParMesh &mesh)
29 : {
30 : return const_cast<mfem::ParMesh &>(
31 11627 : GetParentMesh(const_cast<const mfem::ParMesh &>(mesh)));
32 : }
33 :
34 1187162 : auto GetBdrNeighborAttribute(int i, const mfem::ParMesh &mesh,
35 : mfem::FaceElementTransformations &FET,
36 : mfem::IsoparametricTransformation &T1,
37 : mfem::IsoparametricTransformation &T2)
38 : {
39 : // For internal boundaries, use the element which corresponds to the domain with lower
40 : // attribute number (ensures all boundary elements are aligned).
41 1187162 : BdrGridFunctionCoefficient::GetBdrElementNeighborTransformations(i, mesh, FET, T1, T2);
42 1187162 : return (FET.Elem2 && FET.Elem2->Attribute < FET.Elem1->Attribute) ? FET.Elem2->Attribute
43 1187162 : : FET.Elem1->Attribute;
44 : }
45 :
46 11627 : auto BuildCeedAttributes(const mfem::ParMesh &mesh)
47 : {
48 : // Set up sparse map from global domain attributes to local ones on this process.
49 : // Include ghost elements for all shared faces so we have their material properties
50 : // stored locally. New attributes for libCEED are contiguous and 1-based.
51 : std::unordered_map<int, int> loc_attr;
52 11627 : mfem::FaceElementTransformations FET;
53 11627 : mfem::IsoparametricTransformation T1, T2;
54 : int count = 0;
55 1372557 : for (int i = 0; i < mesh.GetNE(); i++)
56 : {
57 1360930 : const int attr = mesh.GetAttribute(i);
58 1360930 : if (loc_attr.find(attr) == loc_attr.end())
59 : {
60 139091 : loc_attr[attr] = ++count;
61 : }
62 : }
63 196027 : for (int i = 0; i < mesh.GetNSharedFaces(); i++)
64 : {
65 184400 : mesh.GetSharedFaceTransformations(i, FET, T1, T2);
66 184400 : int attr = FET.Elem1->Attribute;
67 184400 : if (loc_attr.find(attr) == loc_attr.end())
68 : {
69 0 : loc_attr[attr] = ++count;
70 : }
71 184400 : attr = FET.Elem2->Attribute;
72 184400 : if (loc_attr.find(attr) == loc_attr.end())
73 : {
74 7082 : loc_attr[attr] = ++count;
75 : }
76 : }
77 11627 : return loc_attr;
78 11627 : }
79 :
80 11627 : auto BuildCeedBdrAttributes(const mfem::ParMesh &mesh)
81 : {
82 : // Set up sparse map from global boundary attributes to local ones on this process. Each
83 : // original global boundary attribute maps to a key-value pairing of global domain
84 : // attributes which neighbor the given boundary and local boundary attributes. New
85 : // attributes for libCEED are contiguous and 1-based.
86 : std::unordered_map<int, std::unordered_map<int, int>> loc_bdr_attr;
87 11627 : mfem::FaceElementTransformations FET;
88 11627 : mfem::IsoparametricTransformation T1, T2;
89 : int count = 0;
90 618521 : for (int i = 0; i < mesh.GetNBE(); i++)
91 : {
92 606894 : const int attr = mesh.GetBdrAttribute(i);
93 606894 : const int nbr_attr = GetBdrNeighborAttribute(i, mesh, FET, T1, T2);
94 : auto &bdr_attr_map = loc_bdr_attr[attr];
95 606894 : if (bdr_attr_map.find(nbr_attr) == bdr_attr_map.end())
96 : {
97 200828 : bdr_attr_map[nbr_attr] = ++count;
98 : }
99 : }
100 11627 : return loc_bdr_attr;
101 11627 : }
102 :
103 45108 : auto GetElementIndices(const mfem::ParMesh &mesh, bool use_bdr, int start, int stop)
104 : {
105 : // Count the number of elements of each type in the local mesh.
106 : std::unordered_map<mfem::Geometry::Type, int> counts;
107 1855164 : for (int i = start; i < stop; i++)
108 : {
109 3620112 : const auto geom = use_bdr ? mesh.GetBdrElementGeometry(i) : mesh.GetElementGeometry(i);
110 : auto it = counts.find(geom);
111 1810056 : if (it == counts.end())
112 : {
113 54674 : counts[geom] = 1;
114 : }
115 : else
116 : {
117 1755382 : it->second++;
118 : }
119 : }
120 :
121 : // Populate the indices arrays for each element geometry.
122 : std::unordered_map<mfem::Geometry::Type, int> offsets;
123 : std::unordered_map<mfem::Geometry::Type, std::vector<int>> element_indices;
124 99782 : for (auto it = counts.begin(); it != counts.end(); ++it)
125 : {
126 54674 : offsets[it->first] = 0;
127 54674 : element_indices[it->first].resize(it->second);
128 : }
129 1855164 : for (int i = start; i < stop; i++)
130 : {
131 3620112 : const auto geom = use_bdr ? mesh.GetBdrElementGeometry(i) : mesh.GetElementGeometry(i);
132 : auto &offset = offsets[geom];
133 : auto &indices = element_indices[geom];
134 1810056 : indices[offset++] = i;
135 : }
136 :
137 45108 : return element_indices;
138 : }
139 :
140 54674 : auto AssembleGeometryData(Ceed ceed, mfem::Geometry::Type geom, std::vector<int> &indices,
141 : const mfem::GridFunction &mesh_nodes, const Vector &elem_attr)
142 : {
143 : const mfem::FiniteElementSpace &mesh_fespace = *mesh_nodes.FESpace();
144 : const mfem::Mesh &mesh = *mesh_fespace.GetMesh();
145 :
146 : ceed::CeedGeomFactorData data;
147 54674 : data.dim = mfem::Geometry::Dimension[geom];
148 54674 : data.space_dim = mesh.SpaceDimension();
149 54674 : data.indices = std::move(indices);
150 : const std::size_t num_elem = data.indices.size();
151 :
152 : // Construct mesh node element restriction and basis.
153 : CeedElemRestriction mesh_restr =
154 54674 : FiniteElementSpace::BuildCeedElemRestriction(mesh_fespace, ceed, geom, data.indices);
155 54674 : CeedBasis mesh_basis = FiniteElementSpace::BuildCeedBasis(mesh_fespace, ceed, geom);
156 : CeedVector mesh_nodes_vec;
157 54674 : ceed::InitCeedVector(mesh_nodes, ceed, &mesh_nodes_vec);
158 : CeedInt num_qpts;
159 54674 : PalaceCeedCall(ceed, CeedBasisGetNumQuadraturePoints(mesh_basis, &num_qpts));
160 :
161 : // Construct element attribute element restriction and basis.
162 : CeedElemRestriction attr_restr;
163 : CeedBasis attr_basis;
164 54674 : PalaceCeedCall(ceed, CeedElemRestrictionCreateStrided(ceed, num_elem, 1, 1, num_elem,
165 : CEED_STRIDES_BACKEND, &attr_restr));
166 : {
167 : // Note: ceed::GetCeedTopology(CEED_TOPOLOGY_LINE) == 1.
168 218696 : mfem::Vector Bt(num_qpts), Gt(num_qpts), qX(num_qpts), qW(num_qpts);
169 54674 : Bt = 1.0;
170 54674 : Gt = 0.0;
171 54674 : qX = 0.0;
172 54674 : qW = 0.0;
173 54674 : PalaceCeedCall(ceed, CeedBasisCreateH1(ceed, CEED_TOPOLOGY_LINE, 1, 1, num_qpts,
174 : Bt.GetData(), Gt.GetData(), qX.GetData(),
175 : qW.GetData(), &attr_basis));
176 : }
177 : CeedVector elem_attr_vec;
178 54674 : ceed::InitCeedVector(elem_attr, ceed, &elem_attr_vec);
179 :
180 : // Allocate storage for geometry factor data (stored as attribute + quadrature weight +
181 : // Jacobian, column-major).
182 54674 : CeedInt geom_data_size = 2 + data.space_dim * data.dim;
183 54674 : PalaceCeedCall(ceed,
184 : CeedVectorCreate(ceed, (CeedSize)num_elem * num_qpts * geom_data_size,
185 : &data.geom_data));
186 54674 : PalaceCeedCall(
187 : ceed, CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, geom_data_size,
188 : (CeedSize)num_elem * num_qpts * geom_data_size,
189 : CEED_STRIDES_BACKEND, &data.geom_data_restr));
190 :
191 : // Compute the required geometry factors at quadrature points.
192 54674 : ceed::AssembleCeedGeometryData(ceed, mesh_restr, mesh_basis, mesh_nodes_vec, attr_restr,
193 : attr_basis, elem_attr_vec, data.geom_data,
194 : data.geom_data_restr);
195 54674 : PalaceCeedCall(ceed, CeedVectorDestroy(&mesh_nodes_vec));
196 54674 : PalaceCeedCall(ceed, CeedElemRestrictionDestroy(&mesh_restr));
197 54674 : PalaceCeedCall(ceed, CeedBasisDestroy(&mesh_basis));
198 54674 : PalaceCeedCall(ceed, CeedVectorDestroy(&elem_attr_vec));
199 54674 : PalaceCeedCall(ceed, CeedElemRestrictionDestroy(&attr_restr));
200 54674 : PalaceCeedCall(ceed, CeedBasisDestroy(&attr_basis));
201 :
202 54674 : return data;
203 : }
204 :
205 22554 : auto BuildCeedGeomFactorData(
206 : const mfem::ParMesh &mesh, const std::unordered_map<int, int> &loc_attr,
207 : const std::unordered_map<int, std::unordered_map<int, int>> &loc_bdr_attr, Ceed ceed)
208 : {
209 : // Create a list of the element indices in the mesh corresponding to a given thread and
210 : // element geometry type and corresponding geometry factor data. libCEED operators will be
211 : // constructed in parallel over threads, where each thread builds a composite operator
212 : // with sub-operators for each geometry.
213 22554 : const std::size_t nt = ceed::internal::NumCeeds();
214 22554 : auto it = std::find(ceed::internal::GetCeedObjects().begin(),
215 22554 : ceed::internal::GetCeedObjects().end(), ceed);
216 22554 : MFEM_VERIFY(it != ceed::internal::GetCeedObjects().end(),
217 : "Unable to find matching Ceed context in BuildCeedGeomFactorData!");
218 22554 : std::size_t i = std::distance(ceed::internal::GetCeedObjects().begin(), it);
219 22554 : mfem::FaceElementTransformations FET;
220 22554 : mfem::IsoparametricTransformation T1, T2;
221 : ceed::GeometryObjectMap<ceed::CeedGeomFactorData> geom_data_map;
222 :
223 : // First domain elements.
224 : {
225 22554 : const int num_elem = mesh.GetNE();
226 22554 : const int stride = (num_elem + nt - 1) / nt;
227 22554 : const int start = i * stride;
228 22554 : const int stop = std::min(start + stride, num_elem);
229 : constexpr bool use_bdr = false;
230 22554 : auto element_indices = GetElementIndices(mesh, use_bdr, start, stop);
231 22554 : auto GetCeedAttribute = [&]() -> std::function<int(int)>
232 : {
233 22554 : if (const auto *submesh = dynamic_cast<const mfem::ParSubMesh *>(&mesh))
234 : {
235 0 : MFEM_VERIFY(submesh->GetFrom() == mfem::SubMesh::From::Boundary,
236 : "Unexpected non-SubMesh object for BuildCeedGeomFactorData with Mesh "
237 : "with (dim, space_dim) = ("
238 : << mesh.Dimension() << ", " << mesh.SpaceDimension() << ")!");
239 0 : return [&, submesh](int i)
240 : {
241 : // Mesh is actually a boundary submesh, so we use the boundary attribute mappings
242 : // from the parent mesh.
243 0 : const int attr = mesh.GetAttribute(i);
244 0 : const int nbr_attr = GetBdrNeighborAttribute(submesh->GetParentElementIDMap()[i],
245 0 : *submesh->GetParent(), FET, T1, T2);
246 : MFEM_ASSERT(loc_bdr_attr.find(attr) != loc_bdr_attr.end() &&
247 : loc_bdr_attr.at(attr).find(nbr_attr) !=
248 : loc_bdr_attr.at(attr).end(),
249 : "Missing libCEED boundary attribute for attribute " << attr << "!");
250 0 : return loc_bdr_attr.at(attr).at(nbr_attr);
251 0 : };
252 : }
253 : else
254 : {
255 22554 : return [&](int i)
256 : {
257 1229788 : const int attr = mesh.GetAttribute(i);
258 : MFEM_ASSERT(loc_attr.find(attr) != loc_attr.end(),
259 : "Missing libCEED domain attribute for attribute " << attr << "!");
260 1229788 : return loc_attr.at(attr);
261 : };
262 : }
263 22554 : }();
264 51021 : for (auto &[geom, indices] : element_indices)
265 : {
266 28467 : Vector elem_attr(indices.size());
267 1258255 : for (std::size_t k = 0; k < indices.size(); k++)
268 : {
269 2459576 : elem_attr[k] = GetCeedAttribute(indices[k]);
270 : }
271 : geom_data_map.emplace(
272 56934 : geom, AssembleGeometryData(ceed, geom, indices, *mesh.GetNodes(), elem_attr));
273 : }
274 : }
275 :
276 : // Then boundary elements (no support for boundary integrators on meshes embedded in
277 : // higher dimensional space for now).
278 22554 : if (mesh.Dimension() == mesh.SpaceDimension())
279 : {
280 22554 : const int nbe = mesh.GetNBE();
281 22554 : const int stride = (nbe + nt - 1) / nt;
282 22554 : const int start = i * stride;
283 22554 : const int stop = std::min(start + stride, nbe);
284 : constexpr bool use_bdr = true;
285 22554 : auto element_indices = GetElementIndices(mesh, use_bdr, start, stop);
286 580268 : auto GetCeedAttribute = [&](int i)
287 : {
288 580268 : const int attr = mesh.GetBdrAttribute(i);
289 580268 : const int nbr_attr = GetBdrNeighborAttribute(i, mesh, FET, T1, T2);
290 : MFEM_ASSERT(loc_bdr_attr.find(attr) != loc_bdr_attr.end() &&
291 : loc_bdr_attr.at(attr).find(nbr_attr) != loc_bdr_attr.at(attr).end(),
292 : "Missing libCEED boundary attribute for attribute " << attr << "!");
293 580268 : return loc_bdr_attr.at(attr).at(nbr_attr);
294 22554 : };
295 48761 : for (auto &[geom, indices] : element_indices)
296 : {
297 26207 : Vector elem_attr(indices.size());
298 606475 : for (std::size_t k = 0; k < indices.size(); k++)
299 : {
300 580268 : elem_attr[k] = GetCeedAttribute(indices[k]);
301 : }
302 : geom_data_map.emplace(
303 52414 : geom, AssembleGeometryData(ceed, geom, indices, *mesh.GetNodes(), elem_attr));
304 : }
305 : }
306 :
307 22554 : return geom_data_map;
308 22554 : }
309 :
310 : } // namespace
311 :
312 : const ceed::GeometryObjectMap<ceed::CeedGeomFactorData> &
313 22692 : Mesh::GetCeedGeomFactorData(Ceed ceed) const
314 : {
315 22692 : MFEM_VERIFY(!loc_attr.empty(),
316 : "Mesh attribute mappings have not been built for GetCeedGeomFactorData!");
317 : auto it = geom_data.find(ceed);
318 : MFEM_ASSERT(it != geom_data.end(), "Unknown Ceed context in GetCeedGeomFactorData!");
319 22692 : auto &geom_data_map = it->second;
320 22692 : if (geom_data_map.empty())
321 : {
322 45108 : geom_data_map = BuildCeedGeomFactorData(*mesh, loc_attr, loc_bdr_attr, ceed);
323 : }
324 22692 : return geom_data_map;
325 : }
326 :
327 25846 : void Mesh::ResetCeedObjects()
328 : {
329 54284 : for (auto &[ceed, geom_data_map] : geom_data)
330 : {
331 83112 : for (auto &[key, val] : geom_data_map)
332 : {
333 54674 : PalaceCeedCall(ceed, CeedVectorDestroy(&val.geom_data));
334 54674 : PalaceCeedCall(ceed, CeedElemRestrictionDestroy(&val.geom_data_restr));
335 : }
336 : }
337 : geom_data.clear();
338 77538 : for (std::size_t i = 0; i < ceed::internal::GetCeedObjects().size(); i++)
339 : {
340 51692 : Ceed ceed = ceed::internal::GetCeedObjects()[i];
341 103384 : geom_data.emplace(ceed, ceed::GeometryObjectMap<ceed::CeedGeomFactorData>());
342 : }
343 25846 : }
344 :
345 11627 : void Mesh::Update()
346 : {
347 : // Attribute mappings, etc. are always constructed for the parent mesh (use boundary
348 : // attribute maps for the domain attributes of a boundary submesh, for example).
349 : auto &parent_mesh = GetParentMesh(*mesh);
350 11627 : parent_mesh.ExchangeFaceNbrData();
351 : loc_attr.clear();
352 : loc_bdr_attr.clear();
353 11627 : loc_attr = BuildCeedAttributes(parent_mesh);
354 11627 : loc_bdr_attr = BuildCeedBdrAttributes(parent_mesh);
355 11627 : ResetCeedObjects();
356 11627 : }
357 :
358 : } // namespace palace
|
