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 "fespace.hpp"
5 :
6 : #include "fem/bilinearform.hpp"
7 : #include "fem/integrator.hpp"
8 : #include "fem/libceed/basis.hpp"
9 : #include "fem/libceed/restriction.hpp"
10 : #include "linalg/rap.hpp"
11 :
12 : namespace palace
13 : {
14 :
15 30742 : CeedBasis FiniteElementSpace::GetCeedBasis(Ceed ceed, mfem::Geometry::Type geom) const
16 : {
17 : auto it = basis.find(ceed);
18 : MFEM_ASSERT(it != basis.end(), "Unknown Ceed context in GetCeedBasis!");
19 : auto &basis_map = it->second;
20 : auto basis_it = basis_map.find(geom);
21 30742 : if (basis_it != basis_map.end())
22 : {
23 9215 : return basis_it->second;
24 : }
25 21527 : return basis_map.emplace(geom, BuildCeedBasis(*this, ceed, geom)).first->second;
26 : }
27 :
28 : CeedElemRestriction
29 32717 : FiniteElementSpace::GetCeedElemRestriction(Ceed ceed, mfem::Geometry::Type geom,
30 : const std::vector<int> &indices) const
31 : {
32 : auto it = restr.find(ceed);
33 : MFEM_ASSERT(it != restr.end(), "Unknown Ceed context in GetCeedElemRestriction!");
34 : auto &restr_map = it->second;
35 : auto restr_it = restr_map.find(geom);
36 32717 : if (restr_it != restr_map.end())
37 : {
38 9239 : return restr_it->second;
39 : }
40 23478 : return restr_map.emplace(geom, BuildCeedElemRestriction(*this, ceed, geom, indices))
41 23478 : .first->second;
42 : }
43 :
44 : CeedElemRestriction
45 2374 : FiniteElementSpace::GetInterpCeedElemRestriction(Ceed ceed, mfem::Geometry::Type geom,
46 : const std::vector<int> &indices) const
47 : {
48 2374 : const mfem::FiniteElement &fe = *GetFEColl().FiniteElementForGeometry(geom);
49 2374 : if (!HasUniqueInterpRestriction(fe))
50 : {
51 1975 : return GetCeedElemRestriction(ceed, geom, indices);
52 : }
53 : auto it = interp_restr.find(ceed);
54 : MFEM_ASSERT(it != interp_restr.end(),
55 : "Unknown Ceed context in GetInterpCeedElemRestriction!");
56 : auto &restr_map = it->second;
57 : auto restr_it = restr_map.find(geom);
58 399 : if (restr_it != restr_map.end())
59 : {
60 6 : return restr_it->second;
61 : }
62 : return restr_map
63 393 : .emplace(geom, BuildCeedElemRestriction(*this, ceed, geom, indices, true, false))
64 393 : .first->second;
65 : }
66 :
67 : CeedElemRestriction
68 1254 : FiniteElementSpace::GetInterpRangeCeedElemRestriction(Ceed ceed, mfem::Geometry::Type geom,
69 : const std::vector<int> &indices) const
70 : {
71 1254 : const mfem::FiniteElement &fe = *GetFEColl().FiniteElementForGeometry(geom);
72 1254 : if (!HasUniqueInterpRangeRestriction(fe))
73 : {
74 1120 : return GetInterpCeedElemRestriction(ceed, geom, indices);
75 : }
76 : auto it = interp_range_restr.find(ceed);
77 : MFEM_ASSERT(it != interp_range_restr.end(),
78 : "Unknown Ceed context in GetInterpRangeCeedElemRestriction!");
79 : auto &restr_map = it->second;
80 : auto restr_it = restr_map.find(geom);
81 134 : if (restr_it != restr_map.end())
82 : {
83 6 : return restr_it->second;
84 : }
85 : return restr_map
86 128 : .emplace(geom, BuildCeedElemRestriction(*this, ceed, geom, indices, true, true))
87 128 : .first->second;
88 : }
89 :
90 41276 : void FiniteElementSpace::ResetCeedObjects()
91 : {
92 82552 : for (auto &[ceed, basis_map] : basis)
93 : {
94 62803 : for (auto &[key, val] : basis_map)
95 : {
96 21527 : PalaceCeedCall(ceed, CeedBasisDestroy(&val));
97 : }
98 : }
99 82552 : for (auto &[ceed, restr_map] : restr)
100 : {
101 64754 : for (auto &[key, val] : restr_map)
102 : {
103 23478 : PalaceCeedCall(ceed, CeedElemRestrictionDestroy(&val));
104 : }
105 : }
106 82552 : for (auto &[ceed, restr_map] : interp_restr)
107 : {
108 41669 : for (auto &[key, val] : restr_map)
109 : {
110 393 : PalaceCeedCall(ceed, CeedElemRestrictionDestroy(&val));
111 : }
112 : }
113 82552 : for (auto &[ceed, restr_map] : interp_range_restr)
114 : {
115 41404 : for (auto &[key, val] : restr_map)
116 : {
117 128 : PalaceCeedCall(ceed, CeedElemRestrictionDestroy(&val));
118 : }
119 : }
120 : basis.clear();
121 : restr.clear();
122 : interp_restr.clear();
123 : interp_range_restr.clear();
124 123828 : for (std::size_t i = 0; i < ceed::internal::GetCeedObjects().size(); i++)
125 : {
126 82552 : Ceed ceed = ceed::internal::GetCeedObjects()[i];
127 82552 : basis.emplace(ceed, ceed::GeometryObjectMap<CeedBasis>());
128 82552 : restr.emplace(ceed, ceed::GeometryObjectMap<CeedElemRestriction>());
129 82552 : interp_restr.emplace(ceed, ceed::GeometryObjectMap<CeedElemRestriction>());
130 165104 : interp_range_restr.emplace(ceed, ceed::GeometryObjectMap<CeedElemRestriction>());
131 : }
132 41276 : }
133 :
134 76201 : CeedBasis FiniteElementSpace::BuildCeedBasis(const mfem::FiniteElementSpace &fespace,
135 : Ceed ceed, mfem::Geometry::Type geom)
136 : {
137 : // Find the appropriate integration rule for the element.
138 76201 : mfem::IsoparametricTransformation T;
139 : const mfem::FiniteElement *fe_nodal =
140 76201 : fespace.GetMesh()->GetNodalFESpace()->FEColl()->FiniteElementForGeometry(geom);
141 76201 : if (!fe_nodal)
142 : {
143 : fe_nodal =
144 0 : fespace.GetMesh()->GetNodalFESpace()->FEColl()->TraceFiniteElementForGeometry(geom);
145 : }
146 : T.SetFE(fe_nodal);
147 76201 : const int q_order = fem::DefaultIntegrationOrder::Get(T);
148 76201 : const mfem::IntegrationRule &ir = mfem::IntRules.Get(geom, q_order);
149 :
150 : // Build the libCEED basis.
151 : CeedBasis val;
152 76201 : const mfem::FiniteElement *fe = fespace.FEColl()->FiniteElementForGeometry(geom);
153 76201 : if (!fe)
154 : {
155 0 : fe = fespace.FEColl()->TraceFiniteElementForGeometry(geom);
156 : }
157 : const int vdim = fespace.GetVDim();
158 76201 : ceed::InitBasis(*fe, ir, vdim, ceed, &val);
159 76201 : return val;
160 76201 : }
161 :
162 78673 : CeedElemRestriction FiniteElementSpace::BuildCeedElemRestriction(
163 : const mfem::FiniteElementSpace &fespace, Ceed ceed, mfem::Geometry::Type geom,
164 : const std::vector<int> &indices, bool is_interp, bool is_interp_range)
165 : {
166 : // Construct the libCEED element restriction for this element type.
167 : CeedElemRestriction val;
168 78673 : const bool use_bdr = (mfem::Geometry::Dimension[geom] != fespace.GetMesh()->Dimension());
169 78673 : ceed::InitRestriction(fespace, indices, use_bdr, is_interp, is_interp_range, ceed, &val);
170 78673 : return val;
171 : }
172 :
173 6 : const Operator &FiniteElementSpace::BuildDiscreteInterpolator() const
174 : {
175 : // Allow finite element spaces to be swapped in their order (intended as deriv(aux) ->
176 : // primal). G is always partially assembled.
177 : const int dim = Dimension();
178 : const bool swap =
179 6 : (aux_fespace->GetFEColl().GetMapType(dim) == GetFEColl().GetDerivMapType(dim));
180 6 : MFEM_VERIFY(!swap, "Incorrect order for primal/auxiliary (test/trial) spaces in discrete "
181 : "interpolator construction!");
182 6 : MFEM_VERIFY(
183 : GetFEColl().GetMapType(dim) == aux_fespace->GetFEColl().GetDerivMapType(dim),
184 : "Unsupported trial/test FE spaces for FiniteElementSpace discrete interpolator!");
185 6 : const FiniteElementSpace &trial_fespace = !swap ? *aux_fespace : *this;
186 : const FiniteElementSpace &test_fespace = !swap ? *this : *aux_fespace;
187 6 : const auto aux_map_type = trial_fespace.GetFEColl().GetMapType(dim);
188 6 : const auto primal_map_type = test_fespace.GetFEColl().GetMapType(dim);
189 6 : if (aux_map_type == mfem::FiniteElement::VALUE &&
190 6 : primal_map_type == mfem::FiniteElement::H_CURL)
191 : {
192 : // Discrete gradient interpolator.
193 : DiscreteLinearOperator interp(trial_fespace, test_fespace);
194 3 : interp.AddDomainInterpolator<GradientInterpolator>();
195 3 : G = std::make_unique<ParOperator>(interp.PartialAssemble(), trial_fespace, test_fespace,
196 6 : true);
197 : }
198 3 : else if (aux_map_type == mfem::FiniteElement::H_CURL &&
199 3 : primal_map_type == mfem::FiniteElement::H_DIV)
200 : {
201 : // Discrete curl interpolator.
202 : DiscreteLinearOperator interp(trial_fespace, test_fespace);
203 3 : interp.AddDomainInterpolator<CurlInterpolator>();
204 3 : G = std::make_unique<ParOperator>(interp.PartialAssemble(), trial_fespace, test_fespace,
205 6 : true);
206 : }
207 0 : else if (aux_map_type == mfem::FiniteElement::H_DIV &&
208 0 : primal_map_type == mfem::FiniteElement::INTEGRAL)
209 : {
210 : // Discrete divergence interpolator.
211 : DiscreteLinearOperator interp(trial_fespace, test_fespace);
212 0 : interp.AddDomainInterpolator<DivergenceInterpolator>();
213 0 : G = std::make_unique<ParOperator>(interp.PartialAssemble(), trial_fespace, test_fespace,
214 0 : true);
215 : }
216 : else
217 : {
218 0 : MFEM_ABORT(
219 : "Unsupported trial/test FE spaces for FiniteElementSpace discrete interpolator!");
220 : }
221 :
222 6 : return *G;
223 : }
224 :
225 0 : const Operator &FiniteElementSpaceHierarchy::BuildProlongationAtLevel(std::size_t l) const
226 : {
227 : // P is always partially assembled.
228 0 : MFEM_VERIFY(l + 1 < GetNumLevels(),
229 : "Can only construct a finite element space prolongation with more than one "
230 : "space in the hierarchy!");
231 0 : if (&fespaces[l]->GetMesh() != &fespaces[l + 1]->GetMesh())
232 : {
233 0 : P[l] = std::make_unique<ParOperator>(
234 0 : std::make_unique<mfem::TransferOperator>(*fespaces[l], *fespaces[l + 1]),
235 0 : *fespaces[l], *fespaces[l + 1], true);
236 : }
237 : else
238 : {
239 : DiscreteLinearOperator p(*fespaces[l], *fespaces[l + 1]);
240 0 : p.AddDomainInterpolator<IdentityInterpolator>();
241 0 : P[l] = std::make_unique<ParOperator>(p.PartialAssemble(), *fespaces[l],
242 0 : *fespaces[l + 1], true);
243 : }
244 :
245 0 : return *P[l];
246 : }
247 :
248 : } // namespace palace
|
