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 "integrator.hpp"
5 :
6 : #include <string>
7 : #include <ceed/backend.h>
8 : #include <mfem.hpp>
9 : #include "utils/diagnostic.hpp"
10 :
11 : PalacePragmaDiagnosticPush
12 : PalacePragmaDiagnosticDisableUnused
13 :
14 : #include "fem/qfunctions/apply_qf.h"
15 : #include "fem/qfunctions/geom_qf.h"
16 :
17 : PalacePragmaDiagnosticPop
18 :
19 : namespace palace::ceed
20 : {
21 :
22 : namespace
23 : {
24 :
25 15407 : void AddQFunctionActiveInputs(unsigned int ops, Ceed ceed, CeedBasis basis,
26 : CeedQFunction qf, std::string name = "u")
27 : {
28 : // Add inputs or outputs with evaluation modes for the active vector of a QFunction.
29 : CeedInt num_comp;
30 15407 : PalaceCeedCall(ceed, CeedBasisGetNumComponents(basis, &num_comp));
31 15407 : if (ops & EvalMode::None)
32 : {
33 0 : PalaceCeedCall(ceed, CeedQFunctionAddInput(qf, name.c_str(), num_comp, CEED_EVAL_NONE));
34 : }
35 15407 : if (ops & EvalMode::Interp)
36 : {
37 : CeedInt q_comp;
38 8878 : PalaceCeedCall(ceed,
39 : CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp));
40 8878 : PalaceCeedCall(
41 : ceed, CeedQFunctionAddInput(qf, name.c_str(), num_comp * q_comp, CEED_EVAL_INTERP));
42 : }
43 15407 : if (ops & EvalMode::Grad)
44 : {
45 : CeedInt q_comp;
46 4226 : PalaceCeedCall(ceed,
47 : CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_GRAD, &q_comp));
48 12678 : PalaceCeedCall(ceed, CeedQFunctionAddInput(qf, (std::string("grad_") + name).c_str(),
49 : num_comp * q_comp, CEED_EVAL_GRAD));
50 : }
51 15407 : if (ops & EvalMode::Div)
52 : {
53 : CeedInt q_comp;
54 582 : PalaceCeedCall(ceed,
55 : CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_DIV, &q_comp));
56 1746 : PalaceCeedCall(ceed, CeedQFunctionAddInput(qf, (std::string("div_") + name).c_str(),
57 : num_comp * q_comp, CEED_EVAL_DIV));
58 : }
59 15407 : if (ops & EvalMode::Curl)
60 : {
61 : CeedInt q_comp;
62 1793 : PalaceCeedCall(ceed,
63 : CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_CURL, &q_comp));
64 5379 : PalaceCeedCall(ceed, CeedQFunctionAddInput(qf, (std::string("curl_") + name).c_str(),
65 : num_comp * q_comp, CEED_EVAL_CURL));
66 : }
67 15407 : }
68 :
69 15407 : void AddQFunctionActiveOutputs(unsigned int ops, Ceed ceed, CeedBasis basis,
70 : CeedQFunction qf, std::string name = "v")
71 : {
72 : // Add inputs or outputs with evaluation modes for the active vector of a QFunction.
73 : CeedInt num_comp;
74 15407 : PalaceCeedCall(ceed, CeedBasisGetNumComponents(basis, &num_comp));
75 15407 : if (ops & EvalMode::None)
76 : {
77 0 : PalaceCeedCall(ceed,
78 : CeedQFunctionAddOutput(qf, name.c_str(), num_comp, CEED_EVAL_NONE));
79 : }
80 15407 : if (ops & EvalMode::Interp)
81 : {
82 : CeedInt q_comp;
83 10244 : PalaceCeedCall(ceed,
84 : CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp));
85 10244 : PalaceCeedCall(ceed, CeedQFunctionAddOutput(qf, name.c_str(), num_comp * q_comp,
86 : CEED_EVAL_INTERP));
87 : }
88 15407 : if (ops & EvalMode::Grad)
89 : {
90 : CeedInt q_comp;
91 2860 : PalaceCeedCall(ceed,
92 : CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_GRAD, &q_comp));
93 8580 : PalaceCeedCall(ceed, CeedQFunctionAddOutput(qf, (std::string("grad_") + name).c_str(),
94 : num_comp * q_comp, CEED_EVAL_GRAD));
95 : }
96 15407 : if (ops & EvalMode::Div)
97 : {
98 : CeedInt q_comp;
99 582 : PalaceCeedCall(ceed,
100 : CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_DIV, &q_comp));
101 1746 : PalaceCeedCall(ceed, CeedQFunctionAddOutput(qf, (std::string("div_") + name).c_str(),
102 : num_comp * q_comp, CEED_EVAL_DIV));
103 : }
104 15407 : if (ops & EvalMode::Curl)
105 : {
106 : CeedInt q_comp;
107 1793 : PalaceCeedCall(ceed,
108 : CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_CURL, &q_comp));
109 5379 : PalaceCeedCall(ceed, CeedQFunctionAddOutput(qf, (std::string("curl_") + name).c_str(),
110 : num_comp * q_comp, CEED_EVAL_CURL));
111 : }
112 15407 : }
113 :
114 30814 : void AddOperatorActiveFields(unsigned int ops, Ceed ceed, CeedElemRestriction restr,
115 : CeedBasis basis, CeedOperator op, const std::string &name,
116 : CeedVector v)
117 : {
118 : // Set active input or output vector fields of an operator.
119 30814 : if (ops & EvalMode::None)
120 : {
121 0 : PalaceCeedCall(ceed, CeedOperatorSetField(op, name.c_str(), restr, CEED_BASIS_NONE, v));
122 : }
123 30814 : if (ops & EvalMode::Interp)
124 : {
125 19122 : PalaceCeedCall(ceed, CeedOperatorSetField(op, name.c_str(), restr, basis, v));
126 : }
127 30814 : if (ops & EvalMode::Grad)
128 : {
129 21258 : PalaceCeedCall(ceed, CeedOperatorSetField(op, (std::string("grad_") + name).c_str(),
130 : restr, basis, v));
131 : }
132 30814 : if (ops & EvalMode::Div)
133 : {
134 3492 : PalaceCeedCall(ceed, CeedOperatorSetField(op, (std::string("div_") + name).c_str(),
135 : restr, basis, v));
136 : }
137 30814 : if (ops & EvalMode::Curl)
138 : {
139 10758 : PalaceCeedCall(ceed, CeedOperatorSetField(op, (std::string("curl_") + name).c_str(),
140 : restr, basis, v));
141 : }
142 30814 : }
143 :
144 : void AddOperatorActiveInputFields(unsigned int ops, Ceed ceed, CeedElemRestriction restr,
145 : CeedBasis basis, CeedOperator op, std::string name = "u",
146 : CeedVector v = CEED_VECTOR_ACTIVE)
147 : {
148 15407 : AddOperatorActiveFields(ops, ceed, restr, basis, op, name, v);
149 15407 : }
150 :
151 : void AddOperatorActiveOutputFields(unsigned int ops, Ceed ceed, CeedElemRestriction restr,
152 : CeedBasis basis, CeedOperator op, std::string name = "v",
153 : CeedVector v = CEED_VECTOR_ACTIVE)
154 : {
155 15407 : AddOperatorActiveFields(ops, ceed, restr, basis, op, name, v);
156 15407 : }
157 :
158 0 : std::vector<CeedInt> QuadratureDataSetup(unsigned int ops, Ceed ceed,
159 : CeedElemRestriction restr, CeedBasis basis,
160 : CeedVector *q_data,
161 : CeedElemRestriction *q_data_restr)
162 : {
163 : // Operator application at each quadrature point should be square, so just use the inputs
164 : // and ignore the outputs.
165 : CeedInt num_comp;
166 0 : PalaceCeedCall(ceed, CeedBasisGetNumComponents(basis, &num_comp));
167 :
168 : std::vector<CeedInt> active_input_sizes;
169 0 : if (ops & EvalMode::None)
170 : {
171 : active_input_sizes.push_back(num_comp);
172 : }
173 0 : if (ops & EvalMode::Interp)
174 : {
175 : CeedInt q_comp;
176 0 : PalaceCeedCall(ceed,
177 : CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp));
178 0 : active_input_sizes.push_back(num_comp * q_comp);
179 : }
180 0 : if (ops & EvalMode::Grad)
181 : {
182 : CeedInt q_comp;
183 0 : PalaceCeedCall(ceed,
184 : CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_GRAD, &q_comp));
185 0 : active_input_sizes.push_back(num_comp * q_comp);
186 : }
187 0 : if (ops & EvalMode::Div)
188 : {
189 : CeedInt q_comp;
190 0 : PalaceCeedCall(ceed,
191 : CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_DIV, &q_comp));
192 0 : active_input_sizes.push_back(num_comp * q_comp);
193 : }
194 0 : if (ops & EvalMode::Curl)
195 : {
196 : CeedInt q_comp;
197 0 : PalaceCeedCall(ceed,
198 : CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_CURL, &q_comp));
199 0 : active_input_sizes.push_back(num_comp * q_comp);
200 : }
201 :
202 : CeedInt num_elem, num_qpts, q_data_size = 0;
203 0 : PalaceCeedCall(ceed, CeedElemRestrictionGetNumElements(restr, &num_elem));
204 0 : PalaceCeedCall(ceed, CeedBasisGetNumQuadraturePoints(basis, &num_qpts));
205 0 : for (auto size : active_input_sizes)
206 : {
207 0 : q_data_size += size * size;
208 : }
209 :
210 0 : PalaceCeedCall(
211 : ceed, CeedVectorCreate(ceed, (CeedSize)num_elem * num_qpts * q_data_size, q_data));
212 0 : PalaceCeedCall(
213 : ceed, CeedElemRestrictionCreateStrided(ceed, num_elem, num_qpts, q_data_size,
214 : (CeedSize)num_elem * num_qpts * q_data_size,
215 : CEED_STRIDES_BACKEND, q_data_restr));
216 :
217 0 : return active_input_sizes;
218 : }
219 :
220 0 : void QuadratureDataAssembly(const std::vector<CeedInt> &qf_active_sizes,
221 : const CeedQFunctionInfo &info, Ceed ceed,
222 : CeedElemRestriction trial_restr, CeedElemRestriction test_restr,
223 : CeedBasis trial_basis, CeedBasis test_basis, CeedVector q_data,
224 : CeedElemRestriction q_data_restr, CeedOperator *op)
225 : {
226 : // Assemble the quadrature data, destroy the operator, and create a new one for the
227 : // actual operator application.
228 0 : PalaceCeedCall(ceed,
229 : CeedOperatorApply(*op, CEED_VECTOR_NONE, q_data, CEED_REQUEST_IMMEDIATE));
230 0 : PalaceCeedCall(ceed, CeedOperatorDestroy(op));
231 :
232 0 : MFEM_VERIFY(!qf_active_sizes.empty() && qf_active_sizes.size() <= 2,
233 : "Invalid number of active QFunction input/output fields ("
234 : << qf_active_sizes.size() << ")!");
235 : CeedQFunction apply_qf;
236 0 : CeedInt qf_size_1 = qf_active_sizes[0],
237 0 : qf_size_2 = (qf_active_sizes.size() > 1) ? qf_active_sizes[1] : 0;
238 0 : switch (10 * qf_size_1 + qf_size_2)
239 : {
240 : case 1:
241 : case 10:
242 0 : PalaceCeedCall(ceed, CeedQFunctionCreateInterior(
243 : ceed, 1, f_apply_1,
244 : PalaceQFunctionRelativePath(f_apply_1_loc), &apply_qf));
245 : break;
246 : case 2:
247 : case 20:
248 0 : PalaceCeedCall(ceed, CeedQFunctionCreateInterior(
249 : ceed, 1, f_apply_2,
250 : PalaceQFunctionRelativePath(f_apply_2_loc), &apply_qf));
251 : break;
252 : case 3:
253 : case 30:
254 0 : PalaceCeedCall(ceed, CeedQFunctionCreateInterior(
255 : ceed, 1, f_apply_3,
256 : PalaceQFunctionRelativePath(f_apply_3_loc), &apply_qf));
257 : break;
258 : case 22:
259 0 : PalaceCeedCall(ceed, CeedQFunctionCreateInterior(
260 : ceed, 1, f_apply_22,
261 : PalaceQFunctionRelativePath(f_apply_22_loc), &apply_qf));
262 : break;
263 : case 33:
264 0 : PalaceCeedCall(ceed, CeedQFunctionCreateInterior(
265 : ceed, 1, f_apply_33,
266 : PalaceQFunctionRelativePath(f_apply_33_loc), &apply_qf));
267 : break;
268 : case 12:
269 0 : PalaceCeedCall(ceed, CeedQFunctionCreateInterior(
270 : ceed, 1, f_apply_12,
271 : PalaceQFunctionRelativePath(f_apply_12_loc), &apply_qf));
272 : break;
273 : case 13:
274 0 : PalaceCeedCall(ceed, CeedQFunctionCreateInterior(
275 : ceed, 1, f_apply_13,
276 : PalaceQFunctionRelativePath(f_apply_13_loc), &apply_qf));
277 : break;
278 : case 21:
279 0 : PalaceCeedCall(ceed, CeedQFunctionCreateInterior(
280 : ceed, 1, f_apply_21,
281 : PalaceQFunctionRelativePath(f_apply_21_loc), &apply_qf));
282 : break;
283 : case 31:
284 0 : PalaceCeedCall(ceed, CeedQFunctionCreateInterior(
285 : ceed, 1, f_apply_31,
286 : PalaceQFunctionRelativePath(f_apply_31_loc), &apply_qf));
287 : break;
288 0 : default:
289 0 : MFEM_ABORT("Invalid number of QFunction input/output components ("
290 : << qf_size_1 << ", " << qf_size_2 << ")!");
291 : apply_qf = nullptr; // Silence compiler warning
292 : }
293 :
294 : // Inputs/outputs.
295 : {
296 : CeedInt q_data_size;
297 0 : PalaceCeedCall(ceed, CeedElemRestrictionGetNumComponents(q_data_restr, &q_data_size));
298 0 : PalaceCeedCall(ceed,
299 : CeedQFunctionAddInput(apply_qf, "q_data", q_data_size, CEED_EVAL_NONE));
300 : }
301 0 : AddQFunctionActiveInputs(info.trial_ops, ceed, trial_basis, apply_qf);
302 0 : AddQFunctionActiveOutputs(info.test_ops, ceed, test_basis, apply_qf);
303 :
304 : // Create the actual operator.
305 0 : PalaceCeedCall(ceed, CeedOperatorCreate(ceed, apply_qf, nullptr, nullptr, op));
306 0 : PalaceCeedCall(ceed, CeedQFunctionDestroy(&apply_qf));
307 :
308 0 : PalaceCeedCall(
309 : ceed, CeedOperatorSetField(*op, "q_data", q_data_restr, CEED_BASIS_NONE, q_data));
310 0 : AddOperatorActiveInputFields(info.trial_ops, ceed, trial_restr, trial_basis, *op);
311 0 : AddOperatorActiveOutputFields(info.test_ops, ceed, test_restr, test_basis, *op);
312 :
313 0 : PalaceCeedCall(ceed, CeedOperatorCheckReady(*op));
314 0 : }
315 :
316 : } // namespace
317 :
318 12191 : int CeedGeometryDataGetSpaceDimension(CeedElemRestriction geom_data_restr, CeedInt dim,
319 : CeedInt *space_dim)
320 : {
321 12191 : if (space_dim)
322 : {
323 : Ceed ceed;
324 : CeedInt geom_data_size;
325 12191 : PalaceCeedCallBackend(CeedElemRestrictionGetCeed(geom_data_restr, &ceed));
326 12191 : PalaceCeedCall(ceed,
327 : CeedElemRestrictionGetNumComponents(geom_data_restr, &geom_data_size));
328 12191 : *space_dim = (geom_data_size - 2) / dim;
329 : MFEM_ASSERT(2 + (*space_dim) * dim == geom_data_size,
330 : "Invalid size for geometry quadrature data!");
331 : }
332 12191 : return CEED_ERROR_SUCCESS;
333 : }
334 :
335 54674 : void AssembleCeedGeometryData(Ceed ceed, CeedElemRestriction mesh_restr,
336 : CeedBasis mesh_basis, CeedVector mesh_nodes,
337 : CeedElemRestriction attr_restr, CeedBasis attr_basis,
338 : CeedVector elem_attr, CeedVector geom_data,
339 : CeedElemRestriction geom_data_restr)
340 : {
341 : CeedInt dim, space_dim, num_qpts;
342 54674 : PalaceCeedCall(ceed, CeedBasisGetDimension(mesh_basis, &dim));
343 54674 : PalaceCeedCall(ceed, CeedBasisGetNumComponents(mesh_basis, &space_dim));
344 54674 : PalaceCeedCall(ceed, CeedBasisGetNumQuadraturePoints(mesh_basis, &num_qpts));
345 :
346 : // Create the QFunction that computes the quadrature data.
347 : CeedQFunction build_qf;
348 54674 : switch (10 * space_dim + dim)
349 : {
350 : case 22:
351 14963 : PalaceCeedCall(ceed, CeedQFunctionCreateInterior(
352 : ceed, 1, f_build_geom_factor_22,
353 : PalaceQFunctionRelativePath(f_build_geom_factor_22_loc),
354 : &build_qf));
355 : break;
356 : case 33:
357 13504 : PalaceCeedCall(ceed, CeedQFunctionCreateInterior(
358 : ceed, 1, f_build_geom_factor_33,
359 : PalaceQFunctionRelativePath(f_build_geom_factor_33_loc),
360 : &build_qf));
361 : break;
362 : case 21:
363 11232 : PalaceCeedCall(ceed, CeedQFunctionCreateInterior(
364 : ceed, 1, f_build_geom_factor_21,
365 : PalaceQFunctionRelativePath(f_build_geom_factor_21_loc),
366 : &build_qf));
367 : break;
368 : case 32:
369 14975 : PalaceCeedCall(ceed, CeedQFunctionCreateInterior(
370 : ceed, 1, f_build_geom_factor_32,
371 : PalaceQFunctionRelativePath(f_build_geom_factor_32_loc),
372 : &build_qf));
373 : break;
374 0 : default:
375 0 : MFEM_ABORT("Invalid value of (dim, space_dim) = ("
376 : << dim << ", " << space_dim << ") for geometry factor quadrature data!");
377 : build_qf = nullptr; // Silence compiler warning
378 : }
379 :
380 : // Inputs/outputs.
381 54674 : PalaceCeedCall(ceed, CeedQFunctionAddInput(build_qf, "attr", 1, CEED_EVAL_INTERP));
382 54674 : PalaceCeedCall(ceed, CeedQFunctionAddInput(build_qf, "q_w", 1, CEED_EVAL_WEIGHT));
383 54674 : PalaceCeedCall(
384 : ceed, CeedQFunctionAddInput(build_qf, "grad_x", space_dim * dim, CEED_EVAL_GRAD));
385 : {
386 : CeedInt geom_data_size;
387 54674 : PalaceCeedCall(ceed,
388 : CeedElemRestrictionGetNumComponents(geom_data_restr, &geom_data_size));
389 54674 : MFEM_VERIFY(geom_data_size == 2 + space_dim * dim,
390 : "Insufficient storage for geometry quadrature data!");
391 54674 : PalaceCeedCall(ceed, CeedQFunctionAddOutput(build_qf, "geom_data", geom_data_size,
392 : CEED_EVAL_NONE));
393 : }
394 :
395 : // Create the operator that builds the quadrature data.
396 : CeedOperator build_op;
397 54674 : PalaceCeedCall(ceed, CeedOperatorCreate(ceed, build_qf, nullptr, nullptr, &build_op));
398 54674 : PalaceCeedCall(ceed, CeedQFunctionDestroy(&build_qf));
399 :
400 54674 : PalaceCeedCall(ceed,
401 : CeedOperatorSetField(build_op, "attr", attr_restr, attr_basis, elem_attr));
402 54674 : PalaceCeedCall(ceed, CeedOperatorSetField(build_op, "q_w", CEED_ELEMRESTRICTION_NONE,
403 : mesh_basis, CEED_VECTOR_NONE));
404 54674 : PalaceCeedCall(ceed, CeedOperatorSetField(build_op, "grad_x", mesh_restr, mesh_basis,
405 : CEED_VECTOR_ACTIVE));
406 54674 : PalaceCeedCall(ceed, CeedOperatorSetField(build_op, "geom_data", geom_data_restr,
407 : CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
408 :
409 54674 : PalaceCeedCall(ceed, CeedOperatorCheckReady(build_op));
410 :
411 : // Compute the quadrature data for the operator.
412 54674 : PalaceCeedCall(
413 : ceed, CeedOperatorApply(build_op, mesh_nodes, geom_data, CEED_REQUEST_IMMEDIATE));
414 54674 : PalaceCeedCall(ceed, CeedOperatorDestroy(&build_op));
415 54674 : }
416 :
417 15407 : void AssembleCeedOperator(const CeedQFunctionInfo &info, void *ctx, std::size_t ctx_size,
418 : Ceed ceed, CeedElemRestriction trial_restr,
419 : CeedElemRestriction test_restr, CeedBasis trial_basis,
420 : CeedBasis test_basis, CeedVector geom_data,
421 : CeedElemRestriction geom_data_restr, CeedOperator *op)
422 : {
423 : // If we are going to be assembling the quadrature data, construct the storage vector for
424 : // it (to be owned by the operator).
425 15407 : CeedVector q_data = nullptr;
426 15407 : CeedElemRestriction q_data_restr = nullptr;
427 : std::vector<CeedInt> qf_active_sizes;
428 15407 : if (info.assemble_q_data)
429 : {
430 0 : qf_active_sizes = QuadratureDataSetup(info.trial_ops, ceed, trial_restr, trial_basis,
431 : &q_data, &q_data_restr);
432 : }
433 :
434 : // Create the QFunction that defines the action of the operator (or its setup).
435 : CeedQFunction apply_qf;
436 15407 : PalaceCeedCall(ceed, CeedQFunctionCreateInterior(ceed, 1, info.apply_qf,
437 : info.apply_qf_path.c_str(), &apply_qf));
438 :
439 : CeedQFunctionContext apply_ctx;
440 15407 : PalaceCeedCall(ceed, CeedQFunctionContextCreate(ceed, &apply_ctx));
441 15407 : PalaceCeedCall(ceed, CeedQFunctionContextSetData(apply_ctx, CEED_MEM_HOST,
442 : CEED_COPY_VALUES, ctx_size, ctx));
443 15407 : PalaceCeedCall(ceed, CeedQFunctionSetContext(apply_qf, apply_ctx));
444 15407 : PalaceCeedCall(ceed, CeedQFunctionContextDestroy(&apply_ctx));
445 :
446 : // Inputs/outputs.
447 : {
448 : CeedInt geom_data_size;
449 15407 : PalaceCeedCall(ceed,
450 : CeedElemRestrictionGetNumComponents(geom_data_restr, &geom_data_size));
451 15407 : PalaceCeedCall(
452 : ceed, CeedQFunctionAddInput(apply_qf, "geom_data", geom_data_size, CEED_EVAL_NONE));
453 : }
454 15407 : if (info.trial_ops & EvalMode::Weight)
455 : {
456 1157 : PalaceCeedCall(ceed, CeedQFunctionAddInput(apply_qf, "q_w", 1, CEED_EVAL_WEIGHT));
457 : }
458 15407 : MFEM_VERIFY(!(info.test_ops & EvalMode::Weight),
459 : "CeedOperator should not have quadrature weight output!");
460 15407 : if (!info.assemble_q_data)
461 : {
462 15407 : AddQFunctionActiveInputs(info.trial_ops, ceed, trial_basis, apply_qf);
463 30814 : AddQFunctionActiveOutputs(info.test_ops, ceed, test_basis, apply_qf);
464 : }
465 : else
466 : {
467 : CeedInt q_data_size;
468 0 : PalaceCeedCall(ceed, CeedElemRestrictionGetNumComponents(q_data_restr, &q_data_size));
469 0 : PalaceCeedCall(ceed,
470 : CeedQFunctionAddOutput(apply_qf, "q_data", q_data_size, CEED_EVAL_NONE));
471 : }
472 :
473 : // Create the operator.
474 15407 : PalaceCeedCall(ceed, CeedOperatorCreate(ceed, apply_qf, nullptr, nullptr, op));
475 15407 : PalaceCeedCall(ceed, CeedQFunctionDestroy(&apply_qf));
476 :
477 15407 : PalaceCeedCall(ceed, CeedOperatorSetField(*op, "geom_data", geom_data_restr,
478 : CEED_BASIS_NONE, geom_data));
479 15407 : if (info.trial_ops & EvalMode::Weight)
480 : {
481 1157 : PalaceCeedCall(ceed, CeedOperatorSetField(*op, "q_w", CEED_ELEMRESTRICTION_NONE,
482 : trial_basis, CEED_VECTOR_NONE));
483 : }
484 15407 : if (!info.assemble_q_data)
485 : {
486 15407 : AddOperatorActiveInputFields(info.trial_ops, ceed, trial_restr, trial_basis, *op);
487 30814 : AddOperatorActiveOutputFields(info.test_ops, ceed, test_restr, test_basis, *op);
488 : }
489 : else
490 : {
491 0 : PalaceCeedCall(ceed, CeedOperatorSetField(*op, "q_data", q_data_restr, CEED_BASIS_NONE,
492 : CEED_VECTOR_ACTIVE));
493 : }
494 :
495 15407 : PalaceCeedCall(ceed, CeedOperatorCheckReady(*op));
496 :
497 : // Assemble the quadrature data and create the actual operator.
498 15407 : if (info.assemble_q_data)
499 : {
500 0 : QuadratureDataAssembly(qf_active_sizes, info, ceed, trial_restr, test_restr,
501 : trial_basis, test_basis, q_data, q_data_restr, op);
502 :
503 : // Cleanup (these are now owned by the operator).
504 0 : PalaceCeedCall(ceed, CeedElemRestrictionDestroy(&q_data_restr));
505 0 : PalaceCeedCall(ceed, CeedVectorDestroy(&q_data));
506 : }
507 15407 : }
508 :
509 1254 : void AssembleCeedInterpolator(Ceed ceed, CeedElemRestriction trial_restr,
510 : CeedElemRestriction test_restr, CeedBasis interp_basis,
511 : CeedOperator *op, CeedOperator *op_t)
512 : {
513 : // Create the QFunction that defines the action of the operator (only an identity as
514 : // element dof multiplicity is handled outside of libCEED).
515 : CeedQFunction apply_qf, apply_qf_t;
516 1254 : PalaceCeedCall(ceed, CeedQFunctionCreateIdentity(ceed, 1, CEED_EVAL_INTERP,
517 : CEED_EVAL_NONE, &apply_qf));
518 1254 : PalaceCeedCall(ceed, CeedQFunctionCreateIdentity(ceed, 1, CEED_EVAL_NONE,
519 : CEED_EVAL_INTERP, &apply_qf_t));
520 :
521 : // Create the operator.
522 1254 : PalaceCeedCall(ceed, CeedOperatorCreate(ceed, apply_qf, nullptr, nullptr, op));
523 1254 : PalaceCeedCall(ceed, CeedQFunctionDestroy(&apply_qf));
524 :
525 1254 : PalaceCeedCall(ceed, CeedOperatorSetField(*op, "input", trial_restr, interp_basis,
526 : CEED_VECTOR_ACTIVE));
527 1254 : PalaceCeedCall(ceed, CeedOperatorSetField(*op, "output", test_restr, CEED_BASIS_NONE,
528 : CEED_VECTOR_ACTIVE));
529 :
530 1254 : PalaceCeedCall(ceed, CeedOperatorCheckReady(*op));
531 :
532 : // Create the transpose operator.
533 1254 : PalaceCeedCall(ceed, CeedOperatorCreate(ceed, apply_qf_t, nullptr, nullptr, op_t));
534 1254 : PalaceCeedCall(ceed, CeedQFunctionDestroy(&apply_qf_t));
535 :
536 1254 : PalaceCeedCall(ceed, CeedOperatorSetField(*op_t, "input", test_restr, CEED_BASIS_NONE,
537 : CEED_VECTOR_ACTIVE));
538 1254 : PalaceCeedCall(ceed, CeedOperatorSetField(*op_t, "output", trial_restr, interp_basis,
539 : CEED_VECTOR_ACTIVE));
540 :
541 1254 : PalaceCeedCall(ceed, CeedOperatorCheckReady(*op_t));
542 1254 : }
543 :
544 0 : void AssembleCeedElementErrorIntegrator(
545 : const CeedQFunctionInfo &info, void *ctx, std::size_t ctx_size, Ceed ceed,
546 : CeedVector input1, CeedVector input2, CeedElemRestriction input1_restr,
547 : CeedElemRestriction input2_restr, CeedBasis input1_basis, CeedBasis input2_basis,
548 : CeedElemRestriction mesh_elem_restr, CeedVector geom_data,
549 : CeedElemRestriction geom_data_restr, CeedOperator *op)
550 : {
551 0 : MFEM_VERIFY(!info.assemble_q_data,
552 : "Quadrature interpolator does not support quadrature data assembly!");
553 :
554 : // Create basis for summing contributions from all quadrature points on the element.
555 : CeedInt num_qpts;
556 0 : PalaceCeedCall(ceed, CeedBasisGetNumQuadraturePoints(input1_basis, &num_qpts));
557 : CeedBasis mesh_elem_basis;
558 : {
559 : // Note: ceed::GetCeedTopology(CEED_TOPOLOGY_LINE) == 1.
560 0 : mfem::Vector Bt(num_qpts), Gt(num_qpts), qX(num_qpts), qW(num_qpts);
561 0 : Bt = 1.0;
562 0 : Gt = 0.0;
563 0 : qX = 0.0;
564 0 : qW = 0.0;
565 0 : PalaceCeedCall(ceed, CeedBasisCreateH1(ceed, CEED_TOPOLOGY_LINE, 1, 1, num_qpts,
566 : Bt.GetData(), Gt.GetData(), qX.GetData(),
567 : qW.GetData(), &mesh_elem_basis));
568 : }
569 :
570 : // Create the QFunction that defines the action of the operator.
571 : CeedQFunction apply_qf;
572 0 : PalaceCeedCall(ceed, CeedQFunctionCreateInterior(ceed, 1, info.apply_qf,
573 : info.apply_qf_path.c_str(), &apply_qf));
574 :
575 : CeedQFunctionContext apply_ctx;
576 0 : PalaceCeedCall(ceed, CeedQFunctionContextCreate(ceed, &apply_ctx));
577 0 : PalaceCeedCall(ceed, CeedQFunctionContextSetData(apply_ctx, CEED_MEM_HOST,
578 : CEED_COPY_VALUES, ctx_size, ctx));
579 0 : PalaceCeedCall(ceed, CeedQFunctionSetContext(apply_qf, apply_ctx));
580 0 : PalaceCeedCall(ceed, CeedQFunctionContextDestroy(&apply_ctx));
581 :
582 : // Inputs/outputs. "Test" operations are the operations for the second input vector.
583 : {
584 : CeedInt geom_data_size;
585 0 : PalaceCeedCall(ceed,
586 : CeedElemRestrictionGetNumComponents(geom_data_restr, &geom_data_size));
587 0 : PalaceCeedCall(
588 : ceed, CeedQFunctionAddInput(apply_qf, "geom_data", geom_data_size, CEED_EVAL_NONE));
589 : }
590 0 : if (info.trial_ops & EvalMode::Weight)
591 : {
592 0 : PalaceCeedCall(ceed, CeedQFunctionAddInput(apply_qf, "q_w", 1, CEED_EVAL_WEIGHT));
593 : }
594 0 : AddQFunctionActiveInputs(info.trial_ops, ceed, input1_basis, apply_qf, "u_1");
595 0 : AddQFunctionActiveInputs(info.test_ops, ceed, input2_basis, apply_qf, "u_2");
596 0 : PalaceCeedCall(ceed, CeedQFunctionAddOutput(apply_qf, "v", 1, CEED_EVAL_INTERP));
597 :
598 : // Create the operator.
599 0 : PalaceCeedCall(ceed, CeedOperatorCreate(ceed, apply_qf, nullptr, nullptr, op));
600 0 : PalaceCeedCall(ceed, CeedQFunctionDestroy(&apply_qf));
601 :
602 0 : PalaceCeedCall(ceed, CeedOperatorSetField(*op, "geom_data", geom_data_restr,
603 : CEED_BASIS_NONE, geom_data));
604 0 : if (info.trial_ops & EvalMode::Weight)
605 : {
606 0 : PalaceCeedCall(ceed, CeedOperatorSetField(*op, "q_w", CEED_ELEMRESTRICTION_NONE,
607 : input1_basis, CEED_VECTOR_NONE));
608 : }
609 0 : AddOperatorActiveInputFields(info.trial_ops, ceed, input1_restr, input1_basis, *op, "u_1",
610 : input1);
611 0 : AddOperatorActiveInputFields(info.test_ops, ceed, input2_restr, input2_basis, *op, "u_2",
612 : input2);
613 0 : PalaceCeedCall(ceed, CeedOperatorSetField(*op, "v", mesh_elem_restr, mesh_elem_basis,
614 : CEED_VECTOR_ACTIVE));
615 :
616 0 : PalaceCeedCall(ceed, CeedOperatorCheckReady(*op));
617 :
618 : // Cleanup (this is now owned by the operator).
619 0 : PalaceCeedCall(ceed, CeedBasisDestroy(&mesh_elem_basis));
620 0 : }
621 :
622 : } // namespace palace::ceed
|
