MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
integrator.hpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12#ifndef MFEM_LIBCEED_INTEG
13#define MFEM_LIBCEED_INTEG
14
16#include "../../fespace.hpp"
17#include "../../gridfunc.hpp"
18#include "operator.hpp"
19#include "coefficient.hpp"
20#include "restriction.hpp"
21#include "util.hpp"
22#include "ceed.hpp"
23
24namespace mfem
25{
26
27namespace ceed
28{
29
30/** The different evaluation modes available for PA and MF CeedIntegrator. */
32
33#ifdef MFEM_USE_CEED
34/** This structure is a template interface for the Assemble methods of
35 PAIntegrator and MFIntegrator. See ceed/mass.cpp for an example. */
37{
38 /** The path to the qFunction header. */
39 const char *header;
40 /** The name of the qFunction to build a partially assembled CeedOperator
41 with a constant Coefficient. */
42 const char *build_func_const;
43 /** The qFunction to build a partially assembled CeedOperator with a constant
44 Coefficient. */
45 CeedQFunctionUser build_qf_const;
46 /** The name of the qFunction to build a partially assembled CeedOperator
47 with a variable Coefficient. */
48 const char *build_func_quad;
49 /** The qFunction to build a partially assembled CeedOperator with a variable
50 Coefficient. */
51 CeedQFunctionUser build_qf_quad;
52 /** The name of the qFunction to apply a partially assembled CeedOperator. */
53 const char *apply_func;
54 /** The qFunction to apply a partially assembled CeedOperator. */
55 CeedQFunctionUser apply_qf;
56 /** The name of the qFunction to apply a matrix-free CeedOperator with a
57 constant Coefficient. */
59 /** The qFunction to apply a matrix-free CeedOperator with a constant
60 Coefficient. */
61 CeedQFunctionUser apply_qf_mf_const;
62 /** The name of the qFunction to apply a matrix-free CeedOperator with a
63 variable Coefficient. */
64 const char *apply_func_mf_quad;
65 /** The qFunction to apply a matrix-free CeedOperator with a variable
66 Coefficient. */
67 CeedQFunctionUser apply_qf_mf_quad;
68 /** The EvalMode on the trial basis functions. */
70 /** The EvalMode on the test basis functions. */
72 /** The size of the data at each quadrature point. */
74};
75#endif
76
77/** This class represent a partially assembled operator using libCEED. */
79{
80#ifdef MFEM_USE_CEED
81protected:
83 CeedElemRestriction trial_restr, test_restr, mesh_restr, restr_i;
84 CeedQFunction build_qfunc, apply_qfunc;
85 CeedVector node_coords, qdata;
87 CeedQFunctionContext build_ctx;
88 CeedOperator build_oper;
89
90public:
92 : Operator(),
93 trial_basis(nullptr), test_basis(nullptr), mesh_basis(nullptr),
94 trial_restr(nullptr), test_restr(nullptr), mesh_restr(nullptr),
95 restr_i(nullptr),
96 build_qfunc(nullptr), apply_qfunc(nullptr), node_coords(nullptr),
97 qdata(nullptr), coeff(nullptr), build_ctx(nullptr), build_oper(nullptr)
98 { }
99
100 /** @brief This method assembles the `PAIntegrator` with the given
101 `CeedOperatorInfo` @a info, an `mfem::FiniteElementSpace` @a fes, an
102 `mfem::IntegrationRule` @a ir, and `mfem::Coefficient` or
103 `mfem::VectorCoefficient` @a Q.
104 The `CeedOperatorInfo` type is expected to inherit from `OperatorInfo`,
105 and contain a `Context` type relevant to the qFunctions.
106
107 @param[in] info is the structure describing the CeedOperator to assemble.
108 @param[in] fes is the finite element space.
109 @param[in] ir is the integration rule for the operator.
110 @param[in] Q is the coefficient from the `Integrator`. */
111 template <typename CeedOperatorInfo, typename CoeffType>
112 void Assemble(CeedOperatorInfo &info,
113 const mfem::FiniteElementSpace &fes,
114 const mfem::IntegrationRule &ir,
115 CoeffType *Q)
116 {
117 Assemble(info, fes, ir, fes.GetNE(), nullptr, Q);
118 }
119
120 /** @brief This method assembles the `PAIntegrator` with the given
121 `CeedOperatorInfo` @a info, an `mfem::FiniteElementSpace` @a fes, an
122 `mfem::IntegrationRule` @a ir, and `mfem::Coefficient` or
123 `mfem::VectorCoefficient` @a Q for the elements given by the indices
124 @a indices.
125 The `CeedOperatorInfo` type is expected to inherit from `OperatorInfo`,
126 and contain a `Context` type relevant to the qFunctions.
127
128 @param[in] info is the structure describing the CeedOperator to assemble.
129 @param[in] fes is the finite element space.
130 @param[in] ir is the integration rule for the operator.
131 @param[in] nelem The number of elements.
132 @param[in] indices The indices of the elements of same type in the
133 `FiniteElementSpace`. If `indices == nullptr`, assumes
134 that the `FiniteElementSpace` is not mixed.
135 @param[in] Q is the coefficient from the `Integrator`. */
136 template <typename CeedOperatorInfo, typename CoeffType>
137 void Assemble(CeedOperatorInfo &info,
138 const mfem::FiniteElementSpace &fes,
139 const mfem::IntegrationRule &ir,
140 int nelem,
141 const int* indices,
142 CoeffType *Q)
143 {
144 Assemble(info, fes, fes, ir, nelem, indices, Q);
145 }
146
147 /** This method assembles the PAIntegrator for mixed forms.
148
149 @param[in] info the `CeedOperatorInfo` describing the `CeedOperator`,
150 the `CeedOperatorInfo` type is expected to inherit from
151 `OperatorInfo` and contain a `Context` type relevant to
152 the qFunctions.
153 @param[in] trial_fes the trial `FiniteElementSpace` for the form,
154 @param[in] test_fes the test `FiniteElementSpace` for the form,
155 @param[in] ir the `IntegrationRule` for the numerical integration,
156 @param[in] Q `Coefficient` or `VectorCoefficient`. */
157 template <typename CeedOperatorInfo, typename CoeffType>
158 void Assemble(CeedOperatorInfo &info,
159 const mfem::FiniteElementSpace &trial_fes,
160 const mfem::FiniteElementSpace &test_fes,
161 const mfem::IntegrationRule &ir,
162 CoeffType *Q)
163 {
164 Assemble(info, trial_fes, test_fes, ir, trial_fes.GetNE(), nullptr, Q);
165 }
166
167 /** This method assembles the PAIntegrator for mixed forms on mixed meshes.
168
169 @param[in] info the `CeedOperatorInfo` describing the `CeedOperator`,
170 the `CeedOperatorInfo` type is expected to inherit from
171 `OperatorInfo` and contain a `Context` type relevant to
172 the qFunctions.
173 @param[in] trial_fes the trial `FiniteElementSpace` for the form,
174 @param[in] test_fes the test `FiniteElementSpace` for the form,
175 @param[in] ir the `IntegrationRule` for the numerical integration,
176 @param[in] nelem The number of elements,
177 @param[in] indices The indices of the elements of same type in the
178 `FiniteElementSpace`. If `indices == nullptr`, assumes
179 that the `FiniteElementSpace` is not mixed,
180 @param[in] Q `Coefficient` or `VectorCoefficient`. */
181 template <typename CeedOperatorInfo, typename CoeffType>
182 void Assemble(CeedOperatorInfo &info,
183 const mfem::FiniteElementSpace &trial_fes,
184 const mfem::FiniteElementSpace &test_fes,
185 const mfem::IntegrationRule &ir,
186 int nelem,
187 const int* indices,
188 CoeffType *Q)
189 {
190 Ceed ceed(internal::ceed);
191 mfem::Mesh &mesh = *trial_fes.GetMesh();
192 MFEM_VERIFY(!(!indices && mesh.GetNumGeometries(mesh.Dimension()) > 1),
193 "Use ceed::MixedIntegrator on mixed meshes.");
194 InitCoefficient(Q, mesh, ir, nelem, indices, coeff, info.ctx);
195 bool const_coeff = coeff->IsConstant();
196 std::string build_func = const_coeff ? info.build_func_const
197 : info.build_func_quad;
198 CeedQFunctionUser build_qf = const_coeff ? info.build_qf_const
199 : info.build_qf_quad;
200 PAOperator op {info.qdatasize, info.header,
201 build_func, build_qf,
202 info.apply_func, info.apply_qf,
203 info.trial_op,
204 info.test_op
205 };
206 CeedInt dim = mesh.SpaceDimension();
207 CeedInt trial_vdim = trial_fes.GetVDim();
208 CeedInt test_vdim = test_fes.GetVDim();
209
210 mesh.EnsureNodes();
211 if ( &trial_fes == &test_fes )
212 {
213 InitBasisAndRestriction(trial_fes, ir, nelem, indices,
214 ceed, &trial_basis, &trial_restr);
217 }
218 else
219 {
220 InitBasisAndRestriction(trial_fes, ir, nelem, indices,
221 ceed, &trial_basis, &trial_restr);
222 InitBasisAndRestriction(test_fes, ir, nelem, indices,
223 ceed, &test_basis, &test_restr);
224 }
225
226 const mfem::FiniteElementSpace *mesh_fes = mesh.GetNodalFESpace();
227 MFEM_VERIFY(mesh_fes, "the Mesh has no nodal FE space");
228 InitBasisAndRestriction(*mesh_fes, ir, nelem, indices,
229 ceed, &mesh_basis, &mesh_restr);
230
231 CeedInt trial_nqpts, test_nqpts;
232 CeedBasisGetNumQuadraturePoints(trial_basis, &trial_nqpts);
233 CeedBasisGetNumQuadraturePoints(test_basis, &test_nqpts);
234 MFEM_VERIFY(trial_nqpts == test_nqpts,
235 "Trial and test basis must have the same number of quadrature"
236 " points.");
237 CeedInt nqpts = trial_nqpts;
238
239 const int qdatasize = op.qdatasize;
240 InitStridedRestriction(*mesh_fes, nelem, nqpts, qdatasize,
241 CEED_STRIDES_BACKEND,
242 &restr_i);
243
245
246 CeedVectorCreate(ceed, nelem * nqpts * qdatasize, &qdata);
247
248 // Context data to be passed to the Q-function.
249 info.ctx.dim = mesh.Dimension();
250 info.ctx.space_dim = mesh.SpaceDimension();
251 info.ctx.vdim = trial_fes.GetVDim();
252
253 std::string qf_file = GetCeedPath() + op.header;
254 std::string qf = qf_file + op.build_func;
255 CeedQFunctionCreateInterior(ceed, 1, op.build_qf, qf.c_str(),
256 &build_qfunc);
257
258 // Create the Q-function that builds the operator (i.e. computes its
259 // quadrature data) and set its context data.
260 if (VariableCoefficient *var_coeff =
261 dynamic_cast<VariableCoefficient*>(coeff))
262 {
263 CeedQFunctionAddInput(build_qfunc, "coeff", coeff->ncomp,
264 var_coeff->emode);
265 }
266 CeedQFunctionAddInput(build_qfunc, "dx", dim * dim, CEED_EVAL_GRAD);
267 CeedQFunctionAddInput(build_qfunc, "weights", 1, CEED_EVAL_WEIGHT);
268 CeedQFunctionAddOutput(build_qfunc, "qdata", qdatasize, CEED_EVAL_NONE);
269
270 CeedQFunctionContextCreate(ceed, &build_ctx);
271 CeedQFunctionContextSetData(build_ctx, CEED_MEM_HOST,
272 CEED_COPY_VALUES,
273 sizeof(info.ctx),
274 &info.ctx);
275 CeedQFunctionSetContext(build_qfunc, build_ctx);
276
277 // Create the operator that builds the quadrature data for the operator.
278 CeedOperatorCreate(ceed, build_qfunc, NULL, NULL, &build_oper);
279 if (GridCoefficient *gridCoeff = dynamic_cast<GridCoefficient*>(coeff))
280 {
281 InitBasisAndRestriction(*gridCoeff->gf.FESpace(), ir,
282 nelem, indices, ceed,
283 &gridCoeff->basis,
284 &gridCoeff->restr);
285 CeedOperatorSetField(build_oper, "coeff", gridCoeff->restr,
286 gridCoeff->basis, gridCoeff->coeffVector);
287 }
288 else if (QuadCoefficient *quadCoeff =
289 dynamic_cast<QuadCoefficient*>(coeff))
290 {
291 const int ncomp = quadCoeff->ncomp;
292 CeedInt strides[3] = {ncomp, 1, ncomp*nqpts};
294 nelem, nqpts, ncomp, strides,
295 &quadCoeff->restr);
296 CeedOperatorSetField(build_oper, "coeff", quadCoeff->restr,
297 CEED_BASIS_NONE, quadCoeff->coeffVector);
298 }
299 CeedOperatorSetField(build_oper, "dx", mesh_restr,
300 mesh_basis, CEED_VECTOR_ACTIVE);
301 CeedOperatorSetField(build_oper, "weights", CEED_ELEMRESTRICTION_NONE,
302 mesh_basis, CEED_VECTOR_NONE);
303 CeedOperatorSetField(build_oper, "qdata", restr_i,
304 CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
305
306 // Compute the quadrature data for the operator.
307 CeedOperatorApply(build_oper, node_coords, qdata, CEED_REQUEST_IMMEDIATE);
308
309 // Create the Q-function that defines the action of the operator.
310 qf = qf_file + op.apply_func;
311 CeedQFunctionCreateInterior(ceed, 1, op.apply_qf, qf.c_str(),
312 &apply_qfunc);
313 // input
314 switch (op.trial_op)
315 {
316 case EvalMode::None:
317 CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim, CEED_EVAL_NONE);
318 break;
319 case EvalMode::Interp:
320 CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim, CEED_EVAL_INTERP);
321 break;
322 case EvalMode::Grad:
323 CeedQFunctionAddInput(apply_qfunc, "gu", trial_vdim*dim, CEED_EVAL_GRAD);
324 break;
326 CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim, CEED_EVAL_INTERP);
327 CeedQFunctionAddInput(apply_qfunc, "gu", trial_vdim*dim, CEED_EVAL_GRAD);
328 break;
329 }
330 // qdata
331 CeedQFunctionAddInput(apply_qfunc, "qdata", qdatasize, CEED_EVAL_NONE);
332 // output
333 switch (op.test_op)
334 {
335 case EvalMode::None:
336 CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim, CEED_EVAL_NONE);
337 break;
338 case EvalMode::Interp:
339 CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim, CEED_EVAL_INTERP);
340 break;
341 case EvalMode::Grad:
342 CeedQFunctionAddOutput(apply_qfunc, "gv", test_vdim*dim, CEED_EVAL_GRAD);
343 break;
345 CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim, CEED_EVAL_INTERP);
346 CeedQFunctionAddOutput(apply_qfunc, "gv", test_vdim*dim, CEED_EVAL_GRAD);
347 break;
348 }
349 CeedQFunctionSetContext(apply_qfunc, build_ctx);
350
351 // Create the operator.
352 CeedOperatorCreate(ceed, apply_qfunc, NULL, NULL, &oper);
353 // input
354 switch (op.trial_op)
355 {
356 case EvalMode::None:
357 CeedOperatorSetField(oper, "u", trial_restr,
358 CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
359 break;
360 case EvalMode::Interp:
361 CeedOperatorSetField(oper, "u", trial_restr, trial_basis, CEED_VECTOR_ACTIVE);
362 break;
363 case EvalMode::Grad:
364 CeedOperatorSetField(oper, "gu", trial_restr, trial_basis, CEED_VECTOR_ACTIVE);
365 break;
367 CeedOperatorSetField(oper, "u", trial_restr, trial_basis, CEED_VECTOR_ACTIVE);
368 CeedOperatorSetField(oper, "gu", trial_restr, trial_basis, CEED_VECTOR_ACTIVE);
369 break;
370 }
371 // qdata
372 CeedOperatorSetField(oper, "qdata", restr_i, CEED_BASIS_NONE,
373 qdata);
374 // output
375 switch (op.test_op)
376 {
377 case EvalMode::None:
378 CeedOperatorSetField(oper, "v", test_restr,
379 CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
380 break;
381 case EvalMode::Interp:
382 CeedOperatorSetField(oper, "v", test_restr, test_basis, CEED_VECTOR_ACTIVE);
383 break;
384 case EvalMode::Grad:
385 CeedOperatorSetField(oper, "gv", test_restr, test_basis, CEED_VECTOR_ACTIVE);
386 break;
388 CeedOperatorSetField(oper, "v", test_restr, test_basis, CEED_VECTOR_ACTIVE);
389 CeedOperatorSetField(oper, "gv", test_restr, test_basis, CEED_VECTOR_ACTIVE);
390 break;
391 }
392
393 CeedVectorCreate(ceed, trial_vdim*trial_fes.GetNDofs(), &u);
394 CeedVectorCreate(ceed, test_vdim*test_fes.GetNDofs(), &v);
395 }
396
398 {
399 CeedQFunctionDestroy(&build_qfunc);
400 CeedQFunctionDestroy(&apply_qfunc);
401 CeedQFunctionContextDestroy(&build_ctx);
402 CeedVectorDestroy(&node_coords);
403 CeedVectorDestroy(&qdata);
404 delete coeff;
405 CeedOperatorDestroy(&build_oper);
406 }
407
408private:
409 /** This structure contains the data to assemble a partially assembled
410 operator with libCEED. */
411 struct PAOperator
412 {
413 /** The number of quadrature data at each quadrature point. */
414 int qdatasize;
415 /** The path to the header containing the functions for libCEED. */
416 std::string header;
417 /** The name of the Qfunction to build the quadrature data. */
418 std::string build_func;
419 /** The Qfunction to build the quadrature data. */
420 CeedQFunctionUser build_qf;
421 /** The name of the Qfunction to apply the operator. */
422 std::string apply_func;
423 /** The Qfunction to apply the operator. */
424 CeedQFunctionUser apply_qf;
425 /** The evaluation mode to apply to the trial function (CEED_EVAL_INTERP,
426 CEED_EVAL_GRAD, etc.) */
427 EvalMode trial_op;
428 /** The evaluation mode to apply to the test function ( CEED_EVAL_INTERP,
429 CEED_EVAL_GRAD, etc.)*/
430 EvalMode test_op;
431 };
432#endif
433};
434
435/** This class represent a matrix-free operator using libCEED. */
437{
438#ifdef MFEM_USE_CEED
439protected:
441 CeedElemRestriction trial_restr, test_restr, mesh_restr, restr_i;
442 CeedQFunction apply_qfunc;
443 CeedVector node_coords, qdata;
445 CeedQFunctionContext build_ctx;
446
447public:
449 : Operator(),
450 trial_basis(nullptr), test_basis(nullptr), mesh_basis(nullptr),
451 trial_restr(nullptr), test_restr(nullptr), mesh_restr(nullptr),
452 restr_i(nullptr),
453 apply_qfunc(nullptr), node_coords(nullptr),
454 qdata(nullptr), coeff(nullptr), build_ctx(nullptr) { }
455
456 /** @brief This method assembles the `MFIntegrator` with the given
457 `CeedOperatorInfo` @a info, an `mfem::FiniteElementSpace` @a fes, an
458 `mfem::IntegrationRule` @a ir, and `mfem::Coefficient` or
459 `mfem::VectorCoefficient` @a Q.
460 The `CeedOperatorInfo` type is expected to inherit from `OperatorInfo`,
461 and contain a `Context` type relevant to the qFunctions.
462
463 @param[in] info is the structure describing the CeedOperator to assemble.
464 @param[in] fes is the finite element space.
465 @param[in] ir is the integration rule for the operator.
466 @param[in] Q is the coefficient from the `Integrator`. */
467 template <typename CeedOperatorInfo, typename CoeffType>
468 void Assemble(CeedOperatorInfo &info,
469 const mfem::FiniteElementSpace &fes,
470 const mfem::IntegrationRule &ir,
471 CoeffType *Q)
472 {
473 Assemble(info, fes, ir, fes.GetNE(), nullptr, Q);
474 }
475
476 /** @brief This method assembles the `MFIntegrator` with the given
477 `CeedOperatorInfo` @a info, an `mfem::FiniteElementSpace` @a fes, an
478 `mfem::IntegrationRule` @a ir, and `mfem::Coefficient` or
479 `mfem::VectorCoefficient` @a Q for the elements given by the indices
480 @a indices.
481 The `CeedOperatorInfo` type is expected to inherit from `OperatorInfo`,
482 and contain a `Context` type relevant to the qFunctions.
483
484 @param[in] info is the structure describing the CeedOperator to assemble.
485 @param[in] fes is the finite element space.
486 @param[in] ir is the integration rule for the operator.
487 @param[in] nelem The number of elements.
488 @param[in] indices The indices of the elements of same type in the
489 `FiniteElementSpace`. If `indices == nullptr`, assumes
490 that the `FiniteElementSpace` is not mixed.
491 @param[in] Q is the coefficient from the `Integrator`. */
492 template <typename CeedOperatorInfo, typename CoeffType>
493 void Assemble(CeedOperatorInfo &info,
494 const mfem::FiniteElementSpace &fes,
495 const mfem::IntegrationRule &ir,
496 int nelem,
497 const int* indices,
498 CoeffType *Q)
499 {
500 Assemble(info, fes, fes, ir, nelem, indices, Q);
501 }
502
503 /** This method assembles the MFIntegrator for mixed forms.
504
505 @param[in] info the `CeedOperatorInfo` describing the `CeedOperator`,
506 the `CeedOperatorInfo` type is expected to inherit from
507 `OperatorInfo` and contain a `Context` type relevant to
508 the qFunctions.
509 @param[in] trial_fes the trial `FiniteElementSpace` for the form,
510 @param[in] test_fes the test `FiniteElementSpace` for the form,
511 @param[in] ir the `IntegrationRule` for the numerical integration,
512 @param[in] Q `Coefficient` or `VectorCoefficient`. */
513 template <typename CeedOperatorInfo, typename CoeffType>
514 void Assemble(CeedOperatorInfo &info,
515 const mfem::FiniteElementSpace &trial_fes,
516 const mfem::FiniteElementSpace &test_fes,
517 const mfem::IntegrationRule &ir,
518 CoeffType *Q)
519 {
520 Assemble(info, trial_fes, test_fes, ir, trial_fes.GetNE(), nullptr, Q);
521 }
522
523 /** This method assembles the MFIntegrator for mixed forms.
524
525 @param[in] info the `CeedOperatorInfo` describing the `CeedOperator`,
526 the `CeedOperatorInfo` type is expected to inherit from
527 `OperatorInfo` and contain a `Context` type relevant to
528 the qFunctions.
529 @param[in] trial_fes the trial `FiniteElementSpace` for the form,
530 @param[in] test_fes the test `FiniteElementSpace` for the form,
531 @param[in] ir the `IntegrationRule` for the numerical integration,
532 @param[in] nelem The number of elements,
533 @param[in] indices The indices of the elements of same type in the
534 `FiniteElementSpace`. If `indices == nullptr`, assumes
535 that the `FiniteElementSpace` is not mixed,
536 @param[in] Q `Coefficient` or `VectorCoefficient`. */
537 template <typename CeedOperatorInfo, typename CoeffType>
538 void Assemble(CeedOperatorInfo &info,
539 const mfem::FiniteElementSpace &trial_fes,
540 const mfem::FiniteElementSpace &test_fes,
541 const mfem::IntegrationRule &ir,
542 int nelem,
543 const int* indices,
544 CoeffType *Q)
545 {
546 Ceed ceed(internal::ceed);
547 Mesh &mesh = *trial_fes.GetMesh();
548 MFEM_VERIFY(!(!indices && mesh.GetNumGeometries(mesh.Dimension()) > 1),
549 "Use ceed::MixedIntegrator on mixed meshes.");
550 InitCoefficient(Q, mesh, ir, nelem, indices, coeff, info.ctx);
551 bool const_coeff = coeff->IsConstant();
552 std::string apply_func = const_coeff ? info.apply_func_mf_const
553 : info.apply_func_mf_quad;
554 CeedQFunctionUser apply_qf = const_coeff ? info.apply_qf_mf_const
555 : info.apply_qf_mf_quad;
556 MFOperator op {info.header,
557 apply_func, apply_qf,
558 info.trial_op,
559 info.test_op
560 };
561
562 CeedInt dim = mesh.SpaceDimension();
563 CeedInt trial_vdim = trial_fes.GetVDim();
564 CeedInt test_vdim = test_fes.GetVDim();
565
566 mesh.EnsureNodes();
567 if ( &trial_fes == &test_fes )
568 {
569 InitBasisAndRestriction(trial_fes, ir, nelem, indices, ceed,
573 }
574 else
575 {
576 InitBasisAndRestriction(trial_fes, ir, nelem, indices, ceed,
578 InitBasisAndRestriction(test_fes, ir, nelem, indices, ceed,
580 }
581
582 const mfem::FiniteElementSpace *mesh_fes = mesh.GetNodalFESpace();
583 MFEM_VERIFY(mesh_fes, "the Mesh has no nodal FE space");
584 InitBasisAndRestriction(*mesh_fes, ir, nelem, indices, ceed, &mesh_basis,
585 &mesh_restr);
586
587 CeedInt trial_nqpts, test_nqpts;
588 CeedBasisGetNumQuadraturePoints(trial_basis, &trial_nqpts);
589 CeedBasisGetNumQuadraturePoints(trial_basis, &test_nqpts);
590 MFEM_VERIFY(trial_nqpts == test_nqpts,
591 "Trial and test basis must have the same number of quadrature"
592 " points.");
593 CeedInt nqpts = trial_nqpts;
594
596
597 // Context data to be passed to the Q-function.
598 info.ctx.dim = mesh.Dimension();
599 info.ctx.space_dim = mesh.SpaceDimension();
600 info.ctx.vdim = trial_fes.GetVDim();
601
602 std::string qf_file = GetCeedPath() + op.header;
603 std::string qf = qf_file + op.apply_func;
604 CeedQFunctionCreateInterior(ceed, 1, op.apply_qf, qf.c_str(),
605 &apply_qfunc);
606
607 // Create the Q-function that builds the operator (i.e. computes its
608 // quadrature data) and set its context data.
609 if (VariableCoefficient *var_coeff =
610 dynamic_cast<VariableCoefficient*>(coeff))
611 {
612 CeedQFunctionAddInput(apply_qfunc, "coeff", coeff->ncomp,
613 var_coeff->emode);
614 }
615 // input
616 switch (op.trial_op)
617 {
618 case EvalMode::None:
619 CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim,
620 CEED_EVAL_NONE);
621 break;
622 case EvalMode::Interp:
623 CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim,
624 CEED_EVAL_INTERP);
625 break;
626 case EvalMode::Grad:
627 CeedQFunctionAddInput(apply_qfunc, "gu", trial_vdim*dim,
628 CEED_EVAL_GRAD);
629 break;
631 CeedQFunctionAddInput(apply_qfunc, "u", trial_vdim,
632 CEED_EVAL_INTERP);
633 CeedQFunctionAddInput(apply_qfunc, "gu", trial_vdim*dim,
634 CEED_EVAL_GRAD);
635 break;
636 }
637 CeedQFunctionAddInput(apply_qfunc, "dx", dim * dim, CEED_EVAL_GRAD);
638 CeedQFunctionAddInput(apply_qfunc, "weights", 1, CEED_EVAL_WEIGHT);
639 // output
640 switch (op.test_op)
641 {
642 case EvalMode::None:
643 CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim,
644 CEED_EVAL_NONE);
645 break;
646 case EvalMode::Interp:
647 CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim,
648 CEED_EVAL_INTERP);
649 break;
650 case EvalMode::Grad:
651 CeedQFunctionAddOutput(apply_qfunc, "gv", test_vdim*dim,
652 CEED_EVAL_GRAD);
653 break;
655 CeedQFunctionAddOutput(apply_qfunc, "v", test_vdim,
656 CEED_EVAL_INTERP);
657 CeedQFunctionAddOutput(apply_qfunc, "gv", test_vdim*dim,
658 CEED_EVAL_GRAD);
659 break;
660 }
661
662 CeedQFunctionContextCreate(ceed, &build_ctx);
663 CeedQFunctionContextSetData(build_ctx, CEED_MEM_HOST,
664 CEED_COPY_VALUES,
665 sizeof(info.ctx),
666 &info.ctx);
667 CeedQFunctionSetContext(apply_qfunc, build_ctx);
668
669 // Create the operator.
670 CeedOperatorCreate(ceed, apply_qfunc, NULL, NULL, &oper);
671 // coefficient
672 if (GridCoefficient *gridCoeff = dynamic_cast<GridCoefficient*>(coeff))
673 {
674 InitBasisAndRestriction(*gridCoeff->gf.FESpace(), ir, nelem, indices,
675 ceed, &gridCoeff->basis, &gridCoeff->restr);
676 CeedOperatorSetField(oper, "coeff", gridCoeff->restr,
677 gridCoeff->basis, gridCoeff->coeffVector);
678 }
679 else if (QuadCoefficient *quadCoeff =
680 dynamic_cast<QuadCoefficient*>(coeff))
681 {
682 const int ncomp = quadCoeff->ncomp;
683 CeedInt strides[3] = {ncomp, 1, ncomp*nqpts};
685 nelem, nqpts, ncomp, strides,
686 &quadCoeff->restr);
687 CeedOperatorSetField(oper, "coeff", quadCoeff->restr,
688 CEED_BASIS_NONE, quadCoeff->coeffVector);
689 }
690 // input
691 switch (op.trial_op)
692 {
693 case EvalMode::None:
694 CeedOperatorSetField(oper, "u", trial_restr,
695 CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
696 break;
697 case EvalMode::Interp:
698 CeedOperatorSetField(oper, "u", trial_restr, trial_basis,
699 CEED_VECTOR_ACTIVE);
700 break;
701 case EvalMode::Grad:
702 CeedOperatorSetField(oper, "gu", trial_restr, trial_basis,
703 CEED_VECTOR_ACTIVE);
704 break;
706 CeedOperatorSetField(oper, "u", trial_restr, trial_basis,
707 CEED_VECTOR_ACTIVE);
708 CeedOperatorSetField(oper, "gu", trial_restr, trial_basis,
709 CEED_VECTOR_ACTIVE);
710 break;
711 }
712 CeedOperatorSetField(oper, "dx", mesh_restr,
714 CeedOperatorSetField(oper, "weights", CEED_ELEMRESTRICTION_NONE,
715 mesh_basis, CEED_VECTOR_NONE);
716 // output
717 switch (op.test_op)
718 {
719 case EvalMode::None:
720 CeedOperatorSetField(oper, "v", test_restr,
721 CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
722 break;
723 case EvalMode::Interp:
724 CeedOperatorSetField(oper, "v", test_restr, test_basis,
725 CEED_VECTOR_ACTIVE);
726 break;
727 case EvalMode::Grad:
728 CeedOperatorSetField(oper, "gv", test_restr, test_basis,
729 CEED_VECTOR_ACTIVE);
730 break;
732 CeedOperatorSetField(oper, "v", test_restr, test_basis,
733 CEED_VECTOR_ACTIVE);
734 CeedOperatorSetField(oper, "gv", test_restr, test_basis,
735 CEED_VECTOR_ACTIVE);
736 break;
737 }
738
739 CeedVectorCreate(ceed, trial_vdim*trial_fes.GetNDofs(), &u);
740 CeedVectorCreate(ceed, test_vdim*test_fes.GetNDofs(), &v);
741 }
742
744 {
745 CeedQFunctionDestroy(&apply_qfunc);
746 CeedQFunctionContextDestroy(&build_ctx);
747 CeedVectorDestroy(&node_coords);
748 CeedVectorDestroy(&qdata);
749 delete coeff;
750 }
751
752private:
753 /** This structure contains the data to assemble a matrix-free operator with
754 libCEED. */
755 struct MFOperator
756 {
757 /** The path to the header containing the functions for libCEED. */
758 std::string header;
759 /** The name of the Qfunction to apply the operator. */
760 std::string apply_func;
761 /** The Qfunction to apply the operator. */
762 CeedQFunctionUser apply_qf;
763 /** The evaluation mode to apply to the trial function (CEED_EVAL_INTERP,
764 CEED_EVAL_GRAD, etc.) */
765 EvalMode trial_op;
766 /** The evaluation mode to apply to the test function ( CEED_EVAL_INTERP,
767 CEED_EVAL_GRAD, etc.) */
768 EvalMode test_op;
769 };
770#endif
771};
772
773} // namespace ceed
774
775} // namespace mfem
776
777#endif // MFEM_LIBCEED_INTEG
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
Definition fespace.hpp:710
int GetNE() const
Returns number of elements in the mesh.
Definition fespace.hpp:740
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:559
int GetVDim() const
Returns vector dimension.
Definition fespace.hpp:706
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
Mesh data type.
Definition mesh.hpp:56
void EnsureNodes()
Make sure that the mesh has valid nodes, i.e. its geometry is described by a vector finite element gr...
Definition mesh.cpp:6159
const FiniteElementSpace * GetNodalFESpace() const
Definition mesh.cpp:6206
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1163
void GetNodes(Vector &node_coord) const
Definition mesh.cpp:8973
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Definition mesh.cpp:6961
void Assemble(CeedOperatorInfo &info, const mfem::FiniteElementSpace &fes, const mfem::IntegrationRule &ir, int nelem, const int *indices, CoeffType *Q)
This method assembles the MFIntegrator with the given CeedOperatorInfo info, an mfem::FiniteElementSp...
CeedElemRestriction trial_restr
void Assemble(CeedOperatorInfo &info, const mfem::FiniteElementSpace &trial_fes, const mfem::FiniteElementSpace &test_fes, const mfem::IntegrationRule &ir, CoeffType *Q)
CeedElemRestriction mesh_restr
CeedQFunctionContext build_ctx
void Assemble(CeedOperatorInfo &info, const mfem::FiniteElementSpace &trial_fes, const mfem::FiniteElementSpace &test_fes, const mfem::IntegrationRule &ir, int nelem, const int *indices, CoeffType *Q)
void Assemble(CeedOperatorInfo &info, const mfem::FiniteElementSpace &fes, const mfem::IntegrationRule &ir, CoeffType *Q)
This method assembles the MFIntegrator with the given CeedOperatorInfo info, an mfem::FiniteElementSp...
CeedElemRestriction restr_i
CeedElemRestriction test_restr
CeedOperator oper
Definition operator.hpp:29
CeedQFunction build_qfunc
CeedElemRestriction restr_i
CeedElemRestriction test_restr
void Assemble(CeedOperatorInfo &info, const mfem::FiniteElementSpace &trial_fes, const mfem::FiniteElementSpace &test_fes, const mfem::IntegrationRule &ir, CoeffType *Q)
CeedQFunctionContext build_ctx
CeedElemRestriction mesh_restr
CeedQFunction apply_qfunc
CeedElemRestriction trial_restr
void Assemble(CeedOperatorInfo &info, const mfem::FiniteElementSpace &fes, const mfem::IntegrationRule &ir, int nelem, const int *indices, CoeffType *Q)
This method assembles the PAIntegrator with the given CeedOperatorInfo info, an mfem::FiniteElementSp...
void Assemble(CeedOperatorInfo &info, const mfem::FiniteElementSpace &fes, const mfem::IntegrationRule &ir, CoeffType *Q)
This method assembles the PAIntegrator with the given CeedOperatorInfo info, an mfem::FiniteElementSp...
void Assemble(CeedOperatorInfo &info, const mfem::FiniteElementSpace &trial_fes, const mfem::FiniteElementSpace &test_fes, const mfem::IntegrationRule &ir, int nelem, const int *indices, CoeffType *Q)
int dim
Definition ex24.cpp:53
const std::string & GetCeedPath()
Return the path to the libCEED q-function headers.
Definition util.cpp:255
void InitVector(const mfem::Vector &v, CeedVector &cv)
Initialize a CeedVector from an mfem::Vector.
Definition util.cpp:75
void InitStridedRestriction(const mfem::FiniteElementSpace &fes, CeedInt nelem, CeedInt nqpts, CeedInt qdatasize, const CeedInt *strides, CeedElemRestriction *restr)
Initialize a strided CeedElemRestriction.
void InitCoefficient(mfem::Coefficient *Q, mfem::Mesh &mesh, const mfem::IntegrationRule &ir, Coefficient *&coeff_ptr, Context &ctx)
Initializes an mfem::ceed::Coefficient coeff_ptr from an mfem::Coefficient Q, an mfem::Mesh mesh,...
void InitBasisAndRestriction(const FiniteElementSpace &fes, const IntegrationRule &irm, Ceed ceed, CeedBasis *basis, CeedElemRestriction *restr)
Initialize a CeedBasis and a CeedElemRestriction based on an mfem::FiniteElementSpace fes,...
Definition util.cpp:93
virtual bool IsConstant() const
CeedQFunctionUser apply_qf_mf_const
const char * apply_func_mf_quad
const char * apply_func_mf_const
CeedQFunctionUser build_qf_const
CeedQFunctionUser apply_qf
CeedQFunctionUser build_qf_quad
CeedQFunctionUser apply_qf_mf_quad
const char * build_func_const
const char * build_func_quad