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 "interpolator.hpp"
5 :
6 : #include <algorithm>
7 : #include "fem/fespace.hpp"
8 : #include "fem/gridfunction.hpp"
9 : #include "utils/communication.hpp"
10 : #include "utils/iodata.hpp"
11 :
12 : namespace palace
13 : {
14 :
15 : namespace
16 : {
17 :
18 : constexpr auto GSLIB_BB_TOL = 0.01; // MFEM defaults, slightly reduced bounding box
19 : constexpr auto GSLIB_NEWTON_TOL = 1.0e-12;
20 :
21 : } // namespace
22 15 : InterpolationOperator::InterpolationOperator(const IoData &iodata,
23 15 : FiniteElementSpace &nd_space)
24 : #if defined(MFEM_USE_GSLIB)
25 15 : : op(nd_space.GetParMesh().GetComm()), v_dim_fes(nd_space.Get().GetVectorDim())
26 : {
27 : auto &mesh = nd_space.GetParMesh();
28 : // Set up probes interpolation. All processes search for all points.
29 15 : if (iodata.domains.postpro.probe.empty())
30 : {
31 15 : return;
32 : }
33 : const int dim = mesh.SpaceDimension();
34 0 : MFEM_VERIFY(
35 : mesh.Dimension() == dim,
36 : "Probe postprocessing functionality requires mesh dimension == space dimension!");
37 0 : const int npts = static_cast<int>(iodata.domains.postpro.probe.size());
38 0 : mfem::Vector xyz(npts * dim);
39 0 : op_idx.resize(npts);
40 : int i = 0;
41 0 : for (const auto &[idx, data] : iodata.domains.postpro.probe)
42 : {
43 0 : for (int d = 0; d < dim; d++)
44 : {
45 : // Use default ordering byNODES.
46 0 : xyz(d * npts + i) = data.center[d];
47 : }
48 0 : op_idx[i++] = idx;
49 : }
50 0 : op.Setup(mesh, GSLIB_BB_TOL, GSLIB_NEWTON_TOL, npts);
51 0 : op.FindPoints(xyz, mfem::Ordering::byNODES);
52 : op.SetDefaultInterpolationValue(0.0);
53 : i = 0;
54 0 : for (const auto &[idx, data] : iodata.domains.postpro.probe)
55 : {
56 0 : if (op.GetCode()[i++] == 2)
57 : {
58 0 : Mpi::Warning(
59 : "Probe {:d} at ({:.3e}) m could not be found!\n Using default value 0.0!\n", idx,
60 0 : fmt::join(iodata.units.Dimensionalize<Units::ValueType::LENGTH>(data.center),
61 : ", "));
62 : }
63 : }
64 0 : }
65 : #else
66 : {
67 : MFEM_CONTRACT_VAR(GSLIB_BB_TOL);
68 : MFEM_CONTRACT_VAR(GSLIB_NEWTON_TOL);
69 : MFEM_VERIFY(iodata.domains.postpro.probe.empty(),
70 : "InterpolationOperator class requires MFEM_USE_GSLIB!");
71 : }
72 : #endif
73 :
74 0 : std::vector<double> InterpolationOperator::ProbeField(const mfem::ParGridFunction &U)
75 : {
76 : #if defined(MFEM_USE_GSLIB)
77 : // Interpolated vector values are returned from GSLIB interpolator with the same ordering
78 : // as the source grid function, which we transform to byVDIM for output.
79 : const int npts = op.GetCode().Size();
80 0 : const int vdim = U.VectorDim();
81 0 : std::vector<double> vals(npts * vdim);
82 0 : if (U.FESpace()->GetOrdering() == mfem::Ordering::byVDIM)
83 : {
84 : mfem::Vector v(vals.data(), npts * vdim);
85 0 : op.Interpolate(U, v);
86 : }
87 : else
88 : {
89 : mfem::Vector v(npts * vdim);
90 0 : op.Interpolate(U, v);
91 0 : for (int d = 0; d < vdim; d++)
92 : {
93 0 : for (int i = 0; i < npts; i++)
94 : {
95 0 : vals[i * vdim + d] = v(d * npts + i);
96 : }
97 : }
98 : }
99 0 : return vals;
100 : #else
101 : MFEM_ABORT("InterpolationOperator class requires MFEM_USE_GSLIB!");
102 : return {};
103 : #endif
104 : }
105 :
106 0 : std::vector<std::complex<double>> InterpolationOperator::ProbeField(const GridFunction &U)
107 : {
108 0 : std::vector<double> vr = ProbeField(U.Real());
109 0 : if (U.HasImag())
110 : {
111 0 : std::vector<double> vi = ProbeField(U.Imag());
112 0 : std::vector<std::complex<double>> vals(vr.size());
113 : std::transform(vr.begin(), vr.end(), vi.begin(), vals.begin(),
114 : [](double xr, double xi) { return std::complex<double>(xr, xi); });
115 : return vals;
116 : }
117 : else
118 : {
119 0 : return {vr.begin(), vr.end()};
120 : }
121 : }
122 :
123 : namespace fem
124 : {
125 :
126 0 : void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V)
127 : {
128 : #if defined(MFEM_USE_GSLIB)
129 : // Generate list of points where the grid function will be evaluated. If the grid function
130 : // to interpolate is an H1 space of the same order as the mesh nodes, we can use the
131 : // mesh node points directly. Otherwise, for a different basis order or type, we generate
132 : // the interpolation points in the physical space manually.
133 : auto &dest_mesh = *V.FESpace()->GetMesh();
134 0 : MFEM_VERIFY(dest_mesh.GetNodes(), "Destination mesh has no nodal FE space!");
135 : const int dim = dest_mesh.SpaceDimension();
136 : mfem::Vector xyz;
137 : mfem::Ordering::Type ordering;
138 : const auto *dest_fec_h1 =
139 0 : dynamic_cast<const mfem::H1_FECollection *>(V.FESpace()->FEColl());
140 0 : const auto *dest_nodes_h1 = dynamic_cast<const mfem::H1_FECollection *>(
141 : dest_mesh.GetNodes()->FESpace()->FEColl());
142 0 : int dest_fespace_order = V.FESpace()->GetMaxElementOrder();
143 0 : int dest_nodes_order = dest_mesh.GetNodes()->FESpace()->GetMaxElementOrder();
144 0 : dest_mesh.GetNodes()->HostRead();
145 0 : if (dest_fec_h1 && dest_nodes_h1 && dest_fespace_order == dest_nodes_order)
146 : {
147 : xyz.MakeRef(*dest_mesh.GetNodes(), 0, dest_mesh.GetNodes()->Size());
148 0 : ordering = dest_mesh.GetNodes()->FESpace()->GetOrdering();
149 : }
150 : else
151 : {
152 : int npts = 0, offset = 0;
153 0 : for (int i = 0; i < dest_mesh.GetNE(); i++)
154 : {
155 0 : npts += V.FESpace()->GetFE(i)->GetNodes().GetNPoints();
156 : }
157 0 : xyz.SetSize(npts * dim);
158 0 : mfem::DenseMatrix pointmat;
159 0 : for (int i = 0; i < dest_mesh.GetNE(); i++)
160 : {
161 0 : const mfem::FiniteElement &fe = *V.FESpace()->GetFE(i);
162 0 : mfem::ElementTransformation &T = *dest_mesh.GetElementTransformation(i);
163 0 : T.Transform(fe.GetNodes(), pointmat);
164 0 : for (int j = 0; j < pointmat.Width(); j++)
165 : {
166 0 : for (int d = 0; d < dim; d++)
167 : {
168 : // Use default ordering byNODES.
169 0 : xyz(d * npts + offset + j) = pointmat(d, j);
170 : }
171 : }
172 0 : offset += pointmat.Width();
173 : }
174 : ordering = mfem::Ordering::byNODES;
175 0 : }
176 0 : const int npts = xyz.Size() / dim;
177 :
178 : // Set up the interpolator.
179 : auto &src_mesh = *U.FESpace()->GetMesh();
180 0 : MFEM_VERIFY(src_mesh.GetNodes(), "Source mesh has no nodal FE space!");
181 0 : auto *src_pmesh = dynamic_cast<mfem::ParMesh *>(&src_mesh);
182 0 : MPI_Comm comm = (src_pmesh) ? src_pmesh->GetComm() : MPI_COMM_SELF;
183 0 : mfem::FindPointsGSLIB op(comm);
184 0 : op.Setup(src_mesh, GSLIB_BB_TOL, GSLIB_NEWTON_TOL, npts);
185 :
186 : // Perform the interpolation and fill the target GridFunction (see MFEM's field-interp
187 : // miniapp).
188 0 : const int vdim = U.VectorDim();
189 0 : mfem::Vector vals(npts * vdim);
190 : op.SetDefaultInterpolationValue(0.0);
191 : op.SetL2AvgType(mfem::FindPointsGSLIB::NONE);
192 0 : op.Interpolate(xyz, U, vals, ordering);
193 : const auto *dest_fec_l2 =
194 0 : dynamic_cast<const mfem::L2_FECollection *>(V.FESpace()->FEColl());
195 0 : if (dest_fec_h1 || dest_fec_l2)
196 : {
197 0 : if (dest_fec_h1 && dest_fespace_order != dest_nodes_order)
198 : {
199 : // H1 with order != mesh order needs to handle duplicated interpolation points.
200 : mfem::Vector elem_vals;
201 : mfem::Array<int> vdofs;
202 : int offset = 0;
203 0 : for (int i = 0; i < dest_mesh.GetNE(); i++)
204 : {
205 0 : const mfem::FiniteElement &fe = *V.FESpace()->GetFE(i);
206 : const int elem_npts = fe.GetNodes().GetNPoints();
207 0 : elem_vals.SetSize(elem_npts * vdim);
208 0 : for (int d = 0; d < vdim; d++)
209 : {
210 0 : for (int j = 0; j < elem_npts; j++)
211 : {
212 : // Arrange element values byNODES to align with GetElementVDofs.
213 : int idx = (U.FESpace()->GetOrdering() == mfem::Ordering::byNODES)
214 0 : ? d * npts + offset + j
215 0 : : (offset + j) * vdim + d;
216 0 : elem_vals(d * elem_npts + j) = vals(idx);
217 : }
218 : }
219 0 : const auto *dof_trans = V.FESpace()->GetElementVDofs(i, vdofs);
220 0 : if (dof_trans)
221 : {
222 : dof_trans->TransformPrimal(elem_vals);
223 : }
224 0 : V.SetSubVector(vdofs, elem_vals);
225 0 : offset += elem_npts;
226 : }
227 : }
228 : else
229 : {
230 : // Otherwise, H1 and L2 copy interpolated values to vdofs.
231 : MFEM_ASSERT(V.Size() == vals.Size(),
232 : "Unexpected size mismatch for interpolated values and grid function!");
233 0 : V = vals;
234 : }
235 : }
236 : else
237 : {
238 : // H(div) or H(curl) use ProjectFromNodes.
239 : mfem::Vector elem_vals, v;
240 : mfem::Array<int> vdofs;
241 : int offset = 0;
242 0 : for (int i = 0; i < dest_mesh.GetNE(); i++)
243 : {
244 0 : const mfem::FiniteElement &fe = *V.FESpace()->GetFE(i);
245 0 : mfem::ElementTransformation &T = *dest_mesh.GetElementTransformation(i);
246 : const int elem_npts = fe.GetNodes().GetNPoints();
247 0 : elem_vals.SetSize(elem_npts * vdim);
248 0 : for (int d = 0; d < vdim; d++)
249 : {
250 0 : for (int j = 0; j < elem_npts; j++)
251 : {
252 : // Arrange element values byVDIM for ProjectFromNodes.
253 : int idx = (U.FESpace()->GetOrdering() == mfem::Ordering::byNODES)
254 0 : ? d * npts + offset + j
255 0 : : (offset + j) * vdim + d;
256 0 : elem_vals(j * vdim + d) = vals(idx);
257 : }
258 : }
259 0 : const auto *dof_trans = V.FESpace()->GetElementVDofs(i, vdofs);
260 0 : v.SetSize(vdofs.Size());
261 0 : fe.ProjectFromNodes(elem_vals, T, v);
262 0 : if (dof_trans)
263 : {
264 : dof_trans->TransformPrimal(v);
265 : }
266 0 : V.SetSubVector(vdofs, v);
267 0 : offset += elem_npts;
268 : }
269 : }
270 : #else
271 : MFEM_ABORT("InterpolateFunction requires MFEM_USE_GSLIB!");
272 : #endif
273 0 : }
274 :
275 0 : void InterpolateFunction(const mfem::Vector &xyz, const mfem::GridFunction &U,
276 : mfem::Vector &vals, mfem::Ordering::Type ordering)
277 : {
278 : #if defined(MFEM_USE_GSLIB)
279 : // Set up the interpolator.
280 : auto &src_mesh = *U.FESpace()->GetMesh();
281 0 : MFEM_VERIFY(src_mesh.GetNodes(), "Source mesh has no nodal FE space!");
282 : const int dim = src_mesh.SpaceDimension();
283 0 : const int npts = xyz.Size() / dim;
284 0 : auto *src_pmesh = dynamic_cast<mfem::ParMesh *>(&src_mesh);
285 0 : MPI_Comm comm = (src_pmesh) ? src_pmesh->GetComm() : MPI_COMM_SELF;
286 0 : mfem::FindPointsGSLIB op(comm);
287 0 : op.Setup(src_mesh, GSLIB_BB_TOL, GSLIB_NEWTON_TOL, npts);
288 :
289 : // Perform the interpolation, with the ordering of the returned values matching the
290 : // ordering of the source grid function.
291 0 : const int vdim = U.VectorDim();
292 0 : MFEM_VERIFY(vals.Size() == npts * vdim, "Incorrect size for interpolated values vector!");
293 : op.SetDefaultInterpolationValue(0.0);
294 : op.SetL2AvgType(mfem::FindPointsGSLIB::NONE);
295 0 : op.Interpolate(xyz, U, vals, ordering);
296 : #else
297 : MFEM_ABORT("InterpolateFunction requires MFEM_USE_GSLIB!");
298 : #endif
299 0 : }
300 :
301 : } // namespace fem
302 :
303 : } // namespace palace
|
