Line data Source code
1 : // Copyright Amazon.com, Inc. or its affiliates. All Rights Reserved.
2 : // SPDX-License-Identifier: Apache-2.0
3 :
4 : #ifndef PALACE_FEM_FESPACE_HPP
5 : #define PALACE_FEM_FESPACE_HPP
6 :
7 : #include <memory>
8 : #include <vector>
9 : #include <mfem.hpp>
10 : #include "fem/libceed/ceed.hpp"
11 : #include "fem/mesh.hpp"
12 : #include "linalg/operator.hpp"
13 : #include "linalg/vector.hpp"
14 :
15 : namespace palace
16 : {
17 :
18 : //
19 : // Wrapper for MFEM's ParFiniteElementSpace class, with extensions for Palace.
20 : //
21 : class FiniteElementSpace
22 : {
23 : private:
24 : // Underlying MFEM object.
25 : mfem::ParFiniteElementSpace fespace;
26 :
27 : // Reference to the underlying mesh object (not owned).
28 : Mesh &mesh;
29 :
30 : // Members for constructing libCEED operators.
31 : mutable ceed::CeedObjectMap<CeedBasis> basis;
32 : mutable ceed::CeedObjectMap<CeedElemRestriction> restr, interp_restr, interp_range_restr;
33 :
34 : // Temporary storage for operator applications.
35 : mutable ComplexVector tx, lx, ly;
36 :
37 : // Members for discrete interpolators from an auxiliary space to a primal space.
38 : mutable const FiniteElementSpace *aux_fespace;
39 : mutable std::unique_ptr<Operator> G;
40 :
41 2374 : bool HasUniqueInterpRestriction(const mfem::FiniteElement &fe) const
42 : {
43 : // For interpolation operators and tensor-product elements, we need native (not
44 : // lexicographic) ordering.
45 : const mfem::TensorBasisElement *tfe =
46 2374 : dynamic_cast<const mfem::TensorBasisElement *>(&fe);
47 2374 : return (tfe && tfe->GetDofMap().Size() > 0 &&
48 2374 : fe.GetRangeType() != mfem::FiniteElement::VECTOR);
49 : }
50 :
51 1254 : bool HasUniqueInterpRangeRestriction(const mfem::FiniteElement &fe) const
52 : {
53 : // The range restriction for interpolation operators needs to use a special
54 : // DofTransformation (not equal to the transpose of the domain restriction).
55 1254 : if (mesh.Dimension() < 3)
56 : {
57 : return false;
58 : }
59 : const auto geom = fe.GetGeomType();
60 681 : const auto *dof_trans = fespace.FEColl()->DofTransformationForGeometry(geom);
61 681 : return (dof_trans && !dof_trans->IsIdentity());
62 : }
63 :
64 : const Operator &BuildDiscreteInterpolator() const;
65 :
66 : public:
67 : template <typename... T>
68 20638 : FiniteElementSpace(Mesh &mesh, T &&...args)
69 20638 : : fespace(&mesh.Get(), std::forward<T>(args)...), mesh(mesh), aux_fespace(nullptr)
70 : {
71 20638 : ResetCeedObjects();
72 20638 : tx.UseDevice(true);
73 20638 : lx.UseDevice(true);
74 20638 : ly.UseDevice(true);
75 20638 : }
76 41339 : virtual ~FiniteElementSpace() { ResetCeedObjects(); }
77 :
78 18 : const auto &Get() const { return fespace; }
79 189 : auto &Get() { return fespace; }
80 :
81 : operator const mfem::ParFiniteElementSpace &() const { return Get(); }
82 : operator mfem::ParFiniteElementSpace &() { return Get(); }
83 :
84 : const auto &GetFEColl() const { return *Get().FEColl(); }
85 : auto &GetFEColl() { return *Get().FEColl(); }
86 :
87 11366 : const auto &GetMesh() const { return mesh; }
88 0 : auto &GetMesh() { return mesh; }
89 :
90 0 : const auto &GetParMesh() const { return mesh.Get(); }
91 84 : auto &GetParMesh() { return mesh.Get(); }
92 :
93 : auto GetVDim() const { return Get().GetVDim(); }
94 : auto GetVSize() const { return Get().GetVSize(); }
95 : auto GlobalVSize() const { return Get().GlobalVSize(); }
96 54 : auto GetTrueVSize() const { return Get().GetTrueVSize(); }
97 : auto GlobalTrueVSize() const { return Get().GlobalTrueVSize(); }
98 5862 : auto Dimension() const { return mesh.Get().Dimension(); }
99 0 : auto SpaceDimension() const { return mesh.Get().SpaceDimension(); }
100 0 : auto GetMaxElementOrder() const { return Get().GetMaxElementOrder(); }
101 :
102 54 : const auto *GetProlongationMatrix() const { return Get().GetProlongationMatrix(); }
103 0 : const auto *GetRestrictionMatrix() const { return Get().GetRestrictionMatrix(); }
104 :
105 : // Return the discrete gradient, curl, or divergence matrix interpolating from the
106 : // auxiliary to the primal space, constructing it on the fly as necessary.
107 6 : const auto &GetDiscreteInterpolator(const FiniteElementSpace &aux_fespace_) const
108 : {
109 6 : if (&aux_fespace_ != aux_fespace)
110 : {
111 : G.reset();
112 6 : aux_fespace = &aux_fespace_;
113 : }
114 6 : return G ? *G : BuildDiscreteInterpolator();
115 : }
116 :
117 : // Return the basis object for elements of the given element geometry type.
118 : CeedBasis GetCeedBasis(Ceed ceed, mfem::Geometry::Type geom) const;
119 :
120 : // Return the element restriction object for the given element set (all with the same
121 : // geometry type).
122 : CeedElemRestriction GetCeedElemRestriction(Ceed ceed, mfem::Geometry::Type geom,
123 : const std::vector<int> &indices) const;
124 :
125 : // If the space has a special element restriction for discrete interpolators, return that.
126 : // Otherwise return the same restriction as given by GetCeedElemRestriction.
127 : CeedElemRestriction GetInterpCeedElemRestriction(Ceed ceed, mfem::Geometry::Type geom,
128 : const std::vector<int> &indices) const;
129 :
130 : // If the space has a special element restriction for the range space of discrete
131 : // interpolators, return that. Otherwise return the same restriction as given by
132 : // GetCeedElemRestriction.
133 : CeedElemRestriction
134 : GetInterpRangeCeedElemRestriction(Ceed ceed, mfem::Geometry::Type geom,
135 : const std::vector<int> &indices) const;
136 :
137 : // Clear the cached basis and element restriction objects owned by the finite element
138 : // space.
139 : void ResetCeedObjects();
140 :
141 : void Update() { ResetCeedObjects(); }
142 :
143 : static CeedBasis BuildCeedBasis(const mfem::FiniteElementSpace &fespace, Ceed ceed,
144 : mfem::Geometry::Type geom);
145 : static CeedElemRestriction
146 : BuildCeedElemRestriction(const mfem::FiniteElementSpace &fespace, Ceed ceed,
147 : mfem::Geometry::Type geom, const std::vector<int> &indices,
148 : bool is_interp = false, bool is_interp_range = false);
149 :
150 : template <typename VecType>
151 : auto &GetTVector() const
152 : {
153 12 : tx.SetSize(GetTrueVSize());
154 : if constexpr (std::is_same<VecType, ComplexVector>::value)
155 : {
156 : return tx;
157 : }
158 : else
159 : {
160 : return tx.Real();
161 : }
162 : }
163 :
164 : template <typename VecType>
165 : auto &GetLVector() const
166 : {
167 18 : lx.SetSize(GetVSize());
168 : if constexpr (std::is_same<VecType, ComplexVector>::value)
169 : {
170 : return lx;
171 : }
172 : else
173 : {
174 : return lx.Real();
175 : }
176 : }
177 :
178 : template <typename VecType>
179 : auto &GetLVector2() const
180 : {
181 18 : ly.SetSize(GetVSize());
182 : if constexpr (std::is_same<VecType, ComplexVector>::value)
183 : {
184 : return ly;
185 : }
186 : else
187 : {
188 : return ly.Real();
189 : }
190 : }
191 :
192 : // Get the associated MPI communicator.
193 : MPI_Comm GetComm() const { return fespace.GetComm(); }
194 : };
195 :
196 : //
197 : // A collection of FiniteElementSpace objects constructed on the same mesh with the ability
198 : // to construct the prolongation operators between them as needed.
199 : //
200 0 : class FiniteElementSpaceHierarchy
201 : {
202 : protected:
203 : std::vector<std::unique_ptr<FiniteElementSpace>> fespaces;
204 : mutable std::vector<std::unique_ptr<Operator>> P;
205 :
206 : const Operator &BuildProlongationAtLevel(std::size_t l) const;
207 :
208 : public:
209 : FiniteElementSpaceHierarchy() = default;
210 63 : FiniteElementSpaceHierarchy(std::unique_ptr<FiniteElementSpace> &&fespace)
211 63 : {
212 63 : AddLevel(std::move(fespace));
213 63 : }
214 :
215 : auto GetNumLevels() const { return fespaces.size(); }
216 :
217 63 : void AddLevel(std::unique_ptr<FiniteElementSpace> &&fespace)
218 : {
219 63 : fespaces.push_back(std::move(fespace));
220 63 : P.push_back(nullptr);
221 63 : }
222 :
223 : auto &GetFESpaceAtLevel(std::size_t l)
224 : {
225 : MFEM_ASSERT(l < GetNumLevels(),
226 : "Out of bounds request for finite element space at level " << l << "!");
227 : return *fespaces[l];
228 : }
229 : const auto &GetFESpaceAtLevel(std::size_t l) const
230 : {
231 : MFEM_ASSERT(l < GetNumLevels(),
232 : "Out of bounds request for finite element space at level " << l << "!");
233 : return *fespaces[l];
234 : }
235 :
236 : auto &GetFinestFESpace()
237 : {
238 : MFEM_ASSERT(GetNumLevels() > 0,
239 : "Out of bounds request for finite element space at level 0!");
240 : return *fespaces.back();
241 : }
242 : const auto &GetFinestFESpace() const
243 : {
244 : MFEM_ASSERT(GetNumLevels() > 0,
245 : "Out of bounds request for finite element space at level 0!");
246 : return *fespaces.back();
247 : }
248 :
249 : const auto &GetProlongationAtLevel(std::size_t l) const
250 : {
251 : MFEM_ASSERT(l + 1 < GetNumLevels(),
252 : "Out of bounds request for finite element space prolongation at level "
253 : << l << "!");
254 0 : return P[l] ? *P[l] : BuildProlongationAtLevel(l);
255 : }
256 :
257 0 : std::vector<const Operator *> GetProlongationOperators() const
258 : {
259 : MFEM_ASSERT(GetNumLevels() > 1,
260 : "Out of bounds request for finite element space prolongation at level 0!");
261 0 : std::vector<const Operator *> P_(GetNumLevels() - 1);
262 0 : for (std::size_t l = 0; l < P_.size(); l++)
263 : {
264 0 : P_[l] = &GetProlongationAtLevel(l);
265 : }
266 0 : return P_;
267 : }
268 :
269 : const auto &GetDiscreteInterpolatorAtLevel(std::size_t l,
270 : const FiniteElementSpace &aux_fespace) const
271 : {
272 0 : return GetFESpaceAtLevel(l).GetDiscreteInterpolator(aux_fespace);
273 : }
274 :
275 : std::vector<const Operator *>
276 0 : GetDiscreteInterpolators(const FiniteElementSpaceHierarchy &aux_fespaces) const
277 : {
278 0 : std::vector<const Operator *> G_(GetNumLevels());
279 0 : G_[0] = nullptr; // No discrete interpolator for coarsest level
280 0 : for (std::size_t l = 1; l < G_.size(); l++)
281 : {
282 0 : G_[l] = &GetDiscreteInterpolatorAtLevel(l, aux_fespaces.GetFESpaceAtLevel(l));
283 : }
284 0 : return G_;
285 : }
286 : };
287 :
288 : } // namespace palace
289 :
290 : #endif // PALACE_FEM_FESPACE_HPP
|
