1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
// *****************************************************************************
/*!
  \file      src/PDE/DGPDE.hpp
  \copyright 2012-2015 J. Bakosi,
             2016-2018 Los Alamos National Security, LLC.,
             2019-2021 Triad National Security, LLC.
             All rights reserved. See the LICENSE file for details.
  \brief     Partial differential equation base for discontinuous Galerkin PDEs
  \details   This file defines a generic partial differential equation (PDE)
    class for PDEs that use discontinuous Galerkin spatial discretization.
    The class uses runtime polymorphism without client-side inheritance:
    inheritance is confined to the internals of the class, invisible to
    client-code. The class exclusively deals with ownership enabling client-side
    value semantics. Credit goes to Sean Parent at Adobe:
    https://github.com/sean-parent/sean-parent.github.com/wiki/
    Papers-and-Presentations.
*/
// *****************************************************************************
#ifndef DGPDE_h
#define DGPDE_h

#include <array>
#include <string>
#include <vector>
#include <memory>
#include <unordered_set>
#include <functional>

#include "Types.hpp"
#include "Fields.hpp"
#include "FaceData.hpp"
#include "UnsMesh.hpp"
#include "Inciter/InputDeck/InputDeck.hpp"
#include "FunctionPrototypes.hpp"
#include "History.hpp"

namespace inciter {

extern ctr::InputDeck g_inputdeck;

using ncomp_t = tk::ncomp_t;
using BCStateFn =
  std::vector< std::pair< std::vector< std::size_t >, tk::StateFn > >;

//! Extract BC configuration ignoring if BC not specified
//! \note A more preferable way of catching errors such as this function
//!   hides is during parsing, so that we don't even get here if BCs are
//!   not correctly specified. For now we simply ignore if BCs are not
//!   specified by allowing empty BC vectors from the user input.
struct ConfigBC {
  BCStateFn& state;    //!< BC state config: sidesets + statefn
  const std::vector< tk::StateFn >& fn;    //!< BC state functions
  std::size_t c;       //!< Counts BC types configured
  //! Constructor
  ConfigBC( BCStateFn& s,
            const std::vector< tk::StateFn >& f ) :
    state(s), fn(f), c(0) {}
  //! Function to call for each BC type
  template< typename U > void operator()( brigand::type_<U> ) {
    std::vector< std::size_t > cfg, v;
    // collect sidesets across all meshes
    for (const auto& ibc : g_inputdeck.get< tag::bc >()) {
      v.insert(v.end(), ibc.get< U >().begin(), ibc.get< U >().end());
    }
    if (v.size() > 0) cfg = v;<--- Variable 'cfg' is assigned a value that is never used.
    Assert( fn.size() > c, "StateFn missing for BC type" );
    state.push_back( { cfg, fn[c++] } );
  }
};

//! State function for invalid/un-configured boundary conditions
[[noreturn]] tk::StateFn::result_type
invalidBC( ncomp_t, const std::vector< EOS >&,
           const std::vector< tk::real >&, tk::real, tk::real, tk::real,
           tk::real, const std::array< tk::real, 3> & );

//! \brief Partial differential equation base for discontinuous Galerkin PDEs
//! \details This class uses runtime polymorphism without client-side
//!   inheritance: inheritance is confined to the internals of the this class,
//!   invisible to client-code. The class exclusively deals with ownership
//!   enabling client-side value semantics. Credit goes to Sean Parent at Adobe:
//!   https://github.com/sean-parent/sean-parent.github.com/wiki/
//!   Papers-and-Presentations. For example client code that models a DGPDE,
//!   see inciter::CompFlow.
class DGPDE {

  private:
    using ncomp_t = tk::ncomp_t;

  public:
    //! Default constructor taking no arguments for Charm++
    explicit DGPDE() = default;

    //! \brief Constructor taking an object modeling Concept.
    //! \details The object of class T comes pre-constructed.
    //! \param[in] x Instantiated object of type T given by the template
    //!   argument.
    template< typename T > explicit DGPDE( T x ) :
      self( std::make_unique< Model<T> >( std::move(x) ) ) {}

    //! \brief Constructor taking a function pointer to a constructor of an
    //!   object modeling Concept.
    //! \details Passing std::function allows late execution of the constructor,
    //!   i.e., as late as inside this class' constructor, and thus usage from
    //!   a factory. Note that there are at least two different ways of using
    //!   this constructor:
    //!   - Bind T's constructor arguments and place it in std::function<T()>
    //!   and passing no arguments as args.... This case then instantiates the
    //!   model via its constructor and stores it in here.
    //!   - Bind a single placeholder argument to T's constructor and pass it in
    //!   as host's args..., which then forwards it to model's constructor. This
    //!   allows late binding, i.e., binding the argument only here.
    //! \see See also the wrapper tk::recordModel() which does the former and
    //!   tk::recordModelLate() which does the latter, both defined in
    //!   src/Base/Factory.h.
    //! \param[in] x Function pointer to a constructor of an object modeling
    //!    Concept.
    //! \param[in] args Zero or more constructor arguments
    template< typename T, typename...Args >
    explicit DGPDE( std::function<T(Args...)> x, Args&&... args ) :
      self( std::make_unique< Model<T> >(
              std::move( x( std::forward<Args>(args)... ) ) ) ) {}

    //! Public interface to find number of primitive quantities for the diff eq
    std::size_t nprim() const
    { return self->nprim(); }

    //! Public interface to find number of materials for the diff eq
    std::size_t nmat() const
    { return self->nmat(); }

    //! Public interface to find Dofs for each equation in pde system
    void numEquationDofs(std::vector< std::size_t >& numEqDof) const
    { return self->numEquationDofs(numEqDof); }

    //! Public interface to determine elements that lie inside the IC box
    void IcBoxElems( const tk::Fields& geoElem,
      std::size_t nielem,
      std::vector< std::unordered_set< std::size_t > >& inbox ) const
    { self->IcBoxElems( geoElem, nielem, inbox ); }

    //! Public interface to setting the initial conditions for the diff eq
    void initialize(
      const tk::Fields& L,
      const std::vector< std::size_t >& inpoel,
      const tk::UnsMesh::Coords& coord,
      const std::vector< std::unordered_set< std::size_t > >& inbox,
      const std::unordered_map< std::size_t, std::set< std::size_t > >&
        elemblkid,
      tk::Fields& unk,
      tk::real t,
      const std::size_t nielem ) const
    { self->initialize( L, inpoel, coord, inbox, elemblkid, unk, t, nielem ); }

    //! Public interface to computing the left-hand side matrix for the diff eq
    void lhs( const tk::Fields& geoElem, tk::Fields& l ) const
    { self->lhs( geoElem, l ); }

    //! Public interface to updating the interface cells for the diff eq
    void updateInterfaceCells( tk::Fields& unk,
                               std::size_t nielem,
                               std::vector< std::size_t >& ndofel ) const
    { self->updateInterfaceCells( unk, nielem, ndofel ); }

    //! Public interface to updating the primitives for the diff eq
    void updatePrimitives( const tk::Fields& unk,
                           const tk::Fields& L,
                           const tk::Fields& geoElem,
                           tk::Fields& prim,
                           std::size_t nielem ) const
    { self->updatePrimitives( unk, L, geoElem, prim, nielem ); }

    //! Public interface to cleaning up trace materials for the diff eq
    void cleanTraceMaterial( tk::real t,
                             const tk::Fields& geoElem,
                             tk::Fields& unk,
                             tk::Fields& prim,
                             std::size_t nielem ) const
    { self->cleanTraceMaterial( t, geoElem, unk, prim, nielem ); }

    //! Public interface to reconstructing the second-order solution
    void reconstruct( tk::real t,
                      const tk::Fields& geoFace,
                      const tk::Fields& geoElem,
                      const inciter::FaceData& fd,
                      const std::map< std::size_t, std::vector< std::size_t > >&
                        esup,
                      const std::vector< std::size_t >& inpoel,
                      const tk::UnsMesh::Coords& coord,
                      tk::Fields& U,
                      tk::Fields& P ) const
    {
      self->reconstruct( t, geoFace, geoElem, fd, esup, inpoel, coord, U, P );
    }

    //! Public interface to limiting the second-order solution
    void limit( tk::real t,
                const tk::Fields& geoFace,
                const tk::Fields& geoElem,
                const inciter::FaceData& fd,
                const std::map< std::size_t, std::vector< std::size_t > >& esup,
                const std::vector< std::size_t >& inpoel,
                const tk::UnsMesh::Coords& coord,
                const std::vector< std::size_t >& ndofel,
                const std::vector< std::size_t >& gid,
                const std::unordered_map< std::size_t, std::size_t >& bid,
                const std::vector< std::vector<tk::real> >& uNodalExtrm,
                const std::vector< std::vector<tk::real> >& pNodalExtrm,
                const std::vector< std::vector<tk::real> >& mtInv,
                tk::Fields& U,
                tk::Fields& P,
                std::vector< std::size_t >& shockmarker ) const
    {
      self->limit( t, geoFace, geoElem, fd, esup, inpoel, coord, ndofel, gid,
                   bid, uNodalExtrm, pNodalExtrm, mtInv, U, P, shockmarker );
    }

    //! Public interface to update the conservative variable solution
    void CPL( const tk::Fields& prim,
              const tk::Fields& geoElem,
              const std::vector< std::size_t >& inpoel,
              const tk::UnsMesh::Coords& coord,
              tk::Fields& unk,
              std::size_t nielem ) const
    {
      self->CPL( prim, geoElem, inpoel, coord, unk, nielem );
    }

    //! Public interface to getting the cell-averaged deformation gradients
    std::array< std::vector< tk::real >, 9 > cellAvgDeformGrad(
      const tk::Fields& U,
      std::size_t nielem ) const
    {
      return self->cellAvgDeformGrad( U, nielem );
    }

    //! Public interface to computing the P1 right-hand side vector
    void rhs( tk::real t,
              const tk::Fields& geoFace,
              const tk::Fields& geoElem,
              const inciter::FaceData& fd,
              const std::vector< std::size_t >& inpoel,
              const std::vector< std::unordered_set< std::size_t > >& boxelems,
              const tk::UnsMesh::Coords& coord,
              const tk::Fields& U,
              const tk::Fields& P,
              const std::vector< std::size_t >& ndofel,
              const tk::real dt,
              tk::Fields& R ) const
    {
      self->rhs( t, geoFace, geoElem, fd, inpoel, boxelems, coord, U, P,
                 ndofel, dt, R );
    }

    //! Evaluate the adaptive indicator and mark the ndof for each element
    void eval_ndof( std::size_t nunk,
                    const tk::UnsMesh::Coords& coord,
                    const std::vector< std::size_t >& inpoel,
                    const inciter::FaceData& fd,
                    const tk::Fields& unk,
                    const tk::Fields& prim,
                    inciter::ctr::PrefIndicatorType indicator,
                    std::size_t ndof,
                    std::size_t ndofmax,
                    tk::real tolref,
                    std::vector< std::size_t >& ndofel ) const
    {
      self->eval_ndof( nunk, coord, inpoel, fd, unk, prim, indicator, ndof,
        ndofmax, tolref, ndofel );
    }

    //! Public interface for computing the minimum time step size
    tk::real dt( const std::array< std::vector< tk::real >, 3 >& coord,
                 const std::vector< std::size_t >& inpoel,
                 const inciter::FaceData& fd,
                 const tk::Fields& geoFace,
                 const tk::Fields& geoElem,
                 const std::vector< std::size_t >& ndofel,
                 const tk::Fields& U,
                 const tk::Fields& P,
                 const std::size_t nielem ) const
    { return self->dt( coord, inpoel, fd, geoFace, geoElem, ndofel, U,
                       P, nielem ); }

    //! Public interface to returning maps of output var functions
    std::map< std::string, tk::GetVarFn > OutVarFn() const
    { return self->OutVarFn(); }

    //! Public interface to returning analytic field output labels
    std::vector< std::string > analyticFieldNames() const
    { return self->analyticFieldNames(); }

    //! Public interface to returning time history field output labels
    std::vector< std::string > histNames() const { return self->histNames(); }

    //! Public interface to returning variable names
    std::vector< std::string > names() const { return self->names(); }

    //! Public interface to returning surface field output
    std::vector< std::vector< tk::real > >
    surfOutput( const std::map< int, std::vector< std::size_t > >& bnd,
                tk::Fields& U ) const
    { return self->surfOutput( bnd, U ); }

    //! Public interface to return point history output
    std::vector< std::vector< tk::real > >
    histOutput( const std::vector< HistData >& h,
      const std::vector< std::size_t >& inpoel,
      const tk::UnsMesh::Coords& coord,
      const tk::Fields& U,
      const tk::Fields& P ) const
    { return self->histOutput( h, inpoel, coord, U, P ); }

    //! Public interface to returning analytic solution
    tk::InitializeFn::result_type
    analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
    { return self->analyticSolution( xi, yi, zi, t ); }

    //! Public interface to returning the analytic solution for conserved vars
    tk::InitializeFn::result_type
    solution( tk::real xi, tk::real yi, tk::real zi, tk::real t ) const
    { return self->solution( xi, yi, zi, t ); }

    //! Public interface to returning the specific total energy
    tk::real
    sp_totalenergy( std::size_t e, const tk::Fields& unk ) const
    { return self->sp_totalenergy( e, unk ); }

    //! Copy assignment
    DGPDE& operator=( const DGPDE& x )
    { DGPDE tmp(x); *this = std::move(tmp); return *this; }
    //! Copy constructor
    DGPDE( const DGPDE& x ) : self( x.self->copy() ) {}
    //! Move assignment
    DGPDE& operator=( DGPDE&& ) noexcept = default;
    //! Move constructor
    DGPDE( DGPDE&& ) noexcept = default;

  private:
    //! \brief Concept is a pure virtual base class specifying the requirements
    //!   of polymorphic objects deriving from it
    struct Concept {
      Concept() = default;
      Concept( const Concept& ) = default;
      virtual ~Concept() = default;
      virtual Concept* copy() const = 0;
      virtual std::size_t nprim() const = 0;
      virtual std::size_t nmat() const = 0;
      virtual void numEquationDofs(std::vector< std::size_t >&) const = 0;
      virtual void IcBoxElems( const tk::Fields&,
        std::size_t,
        std::vector< std::unordered_set< std::size_t > >& ) const = 0;
      virtual void initialize(
        const tk::Fields&,
        const std::vector< std::size_t >&,
        const tk::UnsMesh::Coords&,
        const std::vector< std::unordered_set< std::size_t > >&,
        const std::unordered_map< std::size_t, std::set< std::size_t > >&,
        tk::Fields&,
        tk::real,
        const std::size_t nielem ) const = 0;
      virtual void lhs( const tk::Fields&, tk::Fields& ) const = 0;
      virtual void updateInterfaceCells( tk::Fields&,
                                         std::size_t,
                                         std::vector< std::size_t >& ) const = 0;
      virtual void updatePrimitives( const tk::Fields&,
                                     const tk::Fields&,
                                     const tk::Fields&,
                                     tk::Fields&,
                                     std::size_t ) const = 0;
      virtual void cleanTraceMaterial( tk::real,
                                       const tk::Fields&,
                                       tk::Fields&,
                                       tk::Fields&,
                                       std::size_t ) const = 0;
      virtual void reconstruct( tk::real,
                                const tk::Fields&,
                                const tk::Fields&,
                                const inciter::FaceData&,
                                const std::map< std::size_t,
                                  std::vector< std::size_t > >&,
                                const std::vector< std::size_t >&,
                                const tk::UnsMesh::Coords&,
                                tk::Fields&,
                                tk::Fields& ) const = 0;
      virtual void limit( tk::real,
                          const tk::Fields&,
                          const tk::Fields&,
                          const inciter::FaceData&,
                          const std::map< std::size_t,
                            std::vector< std::size_t > >&,
                          const std::vector< std::size_t >&,
                          const tk::UnsMesh::Coords&,
                          const std::vector< std::size_t >&,
                          const std::vector< std::size_t >&,
                          const std::unordered_map< std::size_t, std::size_t >&,
                          const std::vector< std::vector<tk::real> >&,
                          const std::vector< std::vector<tk::real> >&,
                          const std::vector< std::vector<tk::real> >&,
                          tk::Fields&,
                          tk::Fields&,
                          std::vector< std::size_t >& ) const = 0;
      virtual void CPL( const tk::Fields&,
                        const tk::Fields&,
                        const std::vector< std::size_t >&,
                        const tk::UnsMesh::Coords&,
                        tk::Fields&,
                        std::size_t ) const = 0;
      virtual std::array< std::vector< tk::real >, 9 > cellAvgDeformGrad(
        const tk::Fields&,
        std::size_t ) const = 0;
      virtual void rhs( tk::real,
                        const tk::Fields&,
                        const tk::Fields&,
                        const inciter::FaceData&,
                        const std::vector< std::size_t >&,
                        const std::vector< std::unordered_set< std::size_t > >&,
                        const tk::UnsMesh::Coords&,
                        const tk::Fields&,
                        const tk::Fields&,
                        const std::vector< std::size_t >&,
                        const tk::real,
                        tk::Fields& ) const = 0;
      virtual void eval_ndof( std::size_t,
                              const tk::UnsMesh::Coords&,
                              const std::vector< std::size_t >&,
                              const inciter::FaceData&,
                              const tk::Fields&,
                              const tk::Fields&,
                              inciter::ctr::PrefIndicatorType,
                              std::size_t,
                              std::size_t,
                              tk::real,
                              std::vector< std::size_t >& ) const = 0;
      virtual tk::real dt( const std::array< std::vector< tk::real >, 3 >&,
                           const std::vector< std::size_t >&,
                           const inciter::FaceData&,
                           const tk::Fields&,
                           const tk::Fields&,
                           const std::vector< std::size_t >&,
                           const tk::Fields&,
                           const tk::Fields&,
                           const std::size_t ) const = 0;
      virtual std::map< std::string, tk::GetVarFn > OutVarFn() const = 0;
      virtual std::vector< std::string > analyticFieldNames() const = 0;
      virtual std::vector< std::string > histNames() const = 0;
      virtual std::vector< std::string > names() const = 0;
      virtual std::vector< std::vector< tk::real > > surfOutput(
        const std::map< int, std::vector< std::size_t > >&,
        tk::Fields& ) const = 0;
      virtual std::vector< std::vector< tk::real > > histOutput(
        const std::vector< HistData >&,
        const std::vector< std::size_t >&,
        const tk::UnsMesh::Coords&,
        const tk::Fields&,
        const tk::Fields& ) const = 0;
      virtual tk::InitializeFn::result_type analyticSolution(
        tk::real xi, tk::real yi, tk::real zi, tk::real t ) const = 0;
      virtual tk::InitializeFn::result_type solution(
        tk::real xi, tk::real yi, tk::real zi, tk::real t ) const = 0;
      virtual tk::real sp_totalenergy(
        std::size_t, const tk::Fields& ) const = 0;
    };

    //! \brief Model models the Concept above by deriving from it and overriding
    //!   the virtual functions required by Concept
    template< typename T >
    struct Model : Concept {
      explicit Model( T x ) : data( std::move(x) ) {}
      Concept* copy() const override { return new Model( *this ); }
      std::size_t nprim() const override
      { return data.nprim(); }
      std::size_t nmat() const override
      { return data.nmat(); }
      void numEquationDofs(std::vector< std::size_t >& numEqDof) const override
      { data.numEquationDofs(numEqDof); }
      void IcBoxElems( const tk::Fields& geoElem,
        std::size_t nielem,
        std::vector< std::unordered_set< std::size_t > >& inbox )
      const override { data.IcBoxElems( geoElem, nielem, inbox ); }
      void initialize(
        const tk::Fields& L,
        const std::vector< std::size_t >& inpoel,
        const tk::UnsMesh::Coords& coord,
        const std::vector< std::unordered_set< std::size_t > >& inbox,
        const std::unordered_map< std::size_t, std::set< std::size_t > >&
          elemblkid,
        tk::Fields& unk,
        tk::real t,
        const std::size_t nielem )
      const override { data.initialize( L, inpoel, coord, inbox, elemblkid, unk,
        t, nielem ); }
      void lhs( const tk::Fields& geoElem, tk::Fields& l ) const override
      { data.lhs( geoElem, l ); }
      void updateInterfaceCells( tk::Fields& unk,
                                 std::size_t nielem,
                                 std::vector< std::size_t >& ndofel )
      const override { data.updateInterfaceCells( unk, nielem, ndofel ); }
      void updatePrimitives( const tk::Fields& unk,
                             const tk::Fields& L,
                             const tk::Fields& geoElem,
                             tk::Fields& prim,
                             std::size_t nielem )
      const override { data.updatePrimitives( unk, L, geoElem, prim, nielem ); }
      void cleanTraceMaterial( tk::real t,
                               const tk::Fields& geoElem,
                               tk::Fields& unk,
                               tk::Fields& prim,
                               std::size_t nielem )
      const override { data.cleanTraceMaterial( t, geoElem, unk, prim, nielem ); }
      void reconstruct( tk::real t,
                        const tk::Fields& geoFace,
                        const tk::Fields& geoElem,
                        const inciter::FaceData& fd,
                        const std::map< std::size_t,
                          std::vector< std::size_t > >& esup,
                        const std::vector< std::size_t >& inpoel,
                        const tk::UnsMesh::Coords& coord,
                        tk::Fields& U,
                        tk::Fields& P ) const override
      {
        data.reconstruct( t, geoFace, geoElem, fd, esup, inpoel, coord, U, P );
      }
      void limit( tk::real t,
                  const tk::Fields& geoFace,
                  const tk::Fields& geoElem,
                  const inciter::FaceData& fd,
                  const std::map< std::size_t, std::vector< std::size_t > >&
                    esup,
                  const std::vector< std::size_t >& inpoel,
                  const tk::UnsMesh::Coords& coord,
                  const std::vector< std::size_t >& ndofel,
                  const std::vector< std::size_t >& gid,
                  const std::unordered_map< std::size_t, std::size_t >& bid,
                  const std::vector< std::vector<tk::real> >& uNodalExtrm,
                  const std::vector< std::vector<tk::real> >& pNodalExtrm,
                  const std::vector< std::vector<tk::real> >& mtInv,
                  tk::Fields& U,
                  tk::Fields& P,
                  std::vector< std::size_t >& shockmarker ) const override
      {
        data.limit( t, geoFace, geoElem, fd, esup, inpoel, coord, ndofel, gid,
                    bid, uNodalExtrm, pNodalExtrm, mtInv, U, P, shockmarker );
      }
      void CPL( const tk::Fields& prim,
                const tk::Fields& geoElem,
                const std::vector< std::size_t >& inpoel,
                const tk::UnsMesh::Coords& coord,
                tk::Fields& unk,
                std::size_t nielem ) const override
      {
        data.CPL( prim, geoElem, inpoel, coord, unk, nielem );
      }
      std::array< std::vector< tk::real >, 9 > cellAvgDeformGrad(
        const tk::Fields& U,
        std::size_t nielem ) const override
      {
        return data.cellAvgDeformGrad( U, nielem );
      }
      void rhs(
        tk::real t,
        const tk::Fields& geoFace,
        const tk::Fields& geoElem,
        const inciter::FaceData& fd,
        const std::vector< std::size_t >& inpoel,
        const std::vector< std::unordered_set< std::size_t > >& boxelems,
        const tk::UnsMesh::Coords& coord,
        const tk::Fields& U,
        const tk::Fields& P,
        const std::vector< std::size_t >& ndofel,
        const tk::real dt,
        tk::Fields& R ) const override
      {
        data.rhs( t, geoFace, geoElem, fd, inpoel, boxelems, coord, U, P,
                  ndofel, dt, R );
      }
      void eval_ndof( std::size_t nunk,
                      const tk::UnsMesh::Coords& coord,
                      const std::vector< std::size_t >& inpoel,
                      const inciter::FaceData& fd,
                      const tk::Fields& unk,
                      const tk::Fields& prim,
                      inciter::ctr::PrefIndicatorType indicator,
                      std::size_t ndof,
                      std::size_t ndofmax,
                      tk::real tolref,
                      std::vector< std::size_t >& ndofel ) const override
      { data.eval_ndof( nunk, coord, inpoel, fd, unk, prim, indicator, ndof,
                        ndofmax, tolref, ndofel ); }
      tk::real dt( const std::array< std::vector< tk::real >, 3 >& coord,
                   const std::vector< std::size_t >& inpoel,
                   const inciter::FaceData& fd,
                   const tk::Fields& geoFace,
                   const tk::Fields& geoElem,
                   const std::vector< std::size_t >& ndofel,
                   const tk::Fields& U,
                   const tk::Fields& P,
                   const std::size_t nielem ) const override
      { return data.dt( coord, inpoel, fd, geoFace, geoElem, ndofel,
                        U, P, nielem ); }
      std::map< std::string, tk::GetVarFn > OutVarFn() const override
      { return data.OutVarFn(); }
      std::vector< std::string > analyticFieldNames() const override
      { return data.analyticFieldNames(); }
      std::vector< std::string > histNames() const override
      { return data.histNames(); }
      std::vector< std::string > names() const override
      { return data.names(); }
      std::vector< std::vector< tk::real > > surfOutput(
        const std::map< int, std::vector< std::size_t > >& bnd,
        tk::Fields& U ) const override
      { return data.surfOutput( bnd, U ); }
      std::vector< std::vector< tk::real > > histOutput(
        const std::vector< HistData >& h,
        const std::vector< std::size_t >& inpoel,
        const tk::UnsMesh::Coords& coord,
        const tk::Fields& U,
        const tk::Fields& P ) const override
      { return data.histOutput( h, inpoel, coord, U, P ); }
      tk::InitializeFn::result_type
      analyticSolution( tk::real xi, tk::real yi, tk::real zi, tk::real t )
       const override { return data.analyticSolution( xi, yi, zi, t ); }
      tk::InitializeFn::result_type
      solution( tk::real xi, tk::real yi, tk::real zi, tk::real t )
       const override { return data.solution( xi, yi, zi, t ); }
      tk::real sp_totalenergy( std::size_t e, const tk::Fields& unk )
       const override { return data.sp_totalenergy( e, unk ); }
      T data;
    };

    std::unique_ptr< Concept > self;    //!< Base pointer used polymorphically
};

} // inciter::

#endif // DGPDE_h