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
// *****************************************************************************
/*!
  \file      src/PDE/Reconstruction.cpp
  \copyright 2012-2015 J. Bakosi,
             2016-2018 Los Alamos National Security, LLC.,
             2019-2020 Triad National Security, LLC.
             All rights reserved. See the LICENSE file for details.
  \brief     Reconstruction for reconstructed discontinuous Galerkin methods
  \details   This file contains functions that reconstruct an "n"th order
    polynomial to an "n+1"th order polynomial using a least-squares
    reconstruction procedure.
*/
// *****************************************************************************

#include <array>
#include <vector>
#include <iostream>

#include "Vector.hpp"
#include "Around.hpp"
#include "Base/HashMapReducer.hpp"
#include "Reconstruction.hpp"
#include "MultiMat/MultiMatIndexing.hpp"

void
tk::lhsLeastSq_P0P1(
  const inciter::FaceData& fd,
  const Fields& geoElem,
  const Fields& geoFace,
  std::vector< std::array< std::array< real, 3 >, 3 > >& lhs_ls )
// *****************************************************************************
//  Compute lhs matrix for the least-squares reconstruction
//! \param[in] fd Face connectivity and boundary conditions object
//! \param[in] geoElem Element geometry array
//! \param[in] geoFace Face geometry array
//! \param[in,out] lhs_ls LHS reconstruction matrix
//! \details This function computing the lhs matrix for reconstruction, is
//!   common for primitive and conserved quantities.
// *****************************************************************************
{
  const auto& esuf = fd.Esuf();
  const auto nelem = fd.Esuel().size()/4;<--- Variable 'nelem' is assigned a value that is never used.

  // Compute internal and boundary face contributions
  for (std::size_t f=0; f<esuf.size()/2; ++f)
  {
    Assert( esuf[2*f] > -1, "Left-side element detected as -1" );

    auto el = static_cast< std::size_t >(esuf[2*f]);
    auto er = esuf[2*f+1];

    std::array< real, 3 > geoElemR;
    std::size_t eR(0);

    // A second-order (piecewise linear) solution polynomial can be obtained
    // from the first-order (piecewise constant) FV solutions by using a
    // least-squares (LS) reconstruction process. LS uses the first-order
    // solutions from the cell being processed, and the cells surrounding it.
    // The LS system is obtaining by requiring the following to hold:
    // 'Taylor expansions of solution from cell-i to the centroids of each of
    // its neighboring cells should be equal to the cell average solution on
    // that neighbor cell.'
    // This gives a system of equations for the three second-order DOFs that are
    // to be determined. In 3D tetrahedral meshes, this would give four
    // equations (one for each neighbor )for the three unknown DOFs. This
    // overdetermined system is solved in the least-squares sense using the
    // normal equations approach. The normal equations approach involves
    // pre-multiplying the overdetermined system by the transpose of the system
    // matrix to obtain a square matrix (3x3 in this case).

    // get a 3x3 system by applying the normal equation approach to the
    // least-squares overdetermined system

    if (er > -1) {
    // internal face contribution
      eR = static_cast< std::size_t >(er);
      // Put in cell-centroid coordinates
      geoElemR = {{ geoElem(eR,1,0), geoElem(eR,2,0), geoElem(eR,3,0) }};
    }
    else {
    // boundary face contribution
      // Put in face-centroid coordinates
      geoElemR = {{ geoFace(f,4,0), geoFace(f,5,0), geoFace(f,6,0) }};
    }

    std::array< real, 3 > wdeltax{{ geoElemR[0]-geoElem(el,1,0),
                                    geoElemR[1]-geoElem(el,2,0),
                                    geoElemR[2]-geoElem(el,3,0) }};

    // define a lambda for contributing to lhs matrix
    auto lhs = [&]( std::size_t e ){
    for (std::size_t idir=0; idir<3; ++idir)
      for (std::size_t jdir=0; jdir<3; ++jdir)
        lhs_ls[e][idir][jdir] += wdeltax[idir] * wdeltax[jdir];
    };

    // always add left element contribution (at a boundary face, the internal
    // element is always the left element)
    lhs(el);
    // add right element contribution for internal faces only
    if (er > -1)
      if (eR < nelem) lhs(eR);

  }
}

void
tk::intLeastSq_P0P1(
  ncomp_t ncomp,
  ncomp_t offset,
  const std::size_t rdof,
  const inciter::FaceData& fd,
  const Fields& geoElem,
  const Fields& W,
  std::vector< std::vector< std::array< real, 3 > > >& rhs_ls )
// *****************************************************************************
//  \brief Compute internal surface contributions to rhs vector of the
//    least-squares reconstruction
//! \param[in] ncomp Number of scalar components in this PDE system
//! \param[in] offset Offset this PDE system operates from
//! \param[in] rdof Maximum number of reconstructed degrees of freedom
//! \param[in] fd Face connectivity and boundary conditions object
//! \param[in] geoElem Element geometry array
//! \param[in] W Solution vector to be reconstructed at recent time step
//! \param[in,out] rhs_ls RHS reconstruction vector
//! \details This function computing the internal face contributions to the rhs
//!   vector for reconstruction, is common for primitive and conserved
//!   quantities. If `W` == `U`, compute internal face contributions for the
//!   conserved variables. If `W` == `P`, compute internal face contributions
//!   for the primitive variables.
// *****************************************************************************
{
  const auto& esuf = fd.Esuf();
  const auto nelem = fd.Esuel().size()/4;<--- Variable 'nelem' is assigned a value that is never used.

  // Compute internal face contributions
  for (auto f=fd.Nbfac(); f<esuf.size()/2; ++f)
  {
    Assert( esuf[2*f] > -1 && esuf[2*f+1] > -1, "Interior element detected "
            "as -1" );

    auto el = static_cast< std::size_t >(esuf[2*f]);
    auto er = static_cast< std::size_t >(esuf[2*f+1]);

    // get a 3x3 system by applying the normal equation approach to the
    // least-squares overdetermined system

    // 'wdeltax' is the distance vector between the centroids of this element
    // and its neighbor
    std::array< real, 3 > wdeltax{{ geoElem(er,1,0)-geoElem(el,1,0),
                                    geoElem(er,2,0)-geoElem(el,2,0),
                                    geoElem(er,3,0)-geoElem(el,3,0) }};

    for (std::size_t idir=0; idir<3; ++idir)
    {
      // rhs vector
      for (ncomp_t c=0; c<ncomp; ++c)
      {
        auto mark = c*rdof;
        rhs_ls[el][c][idir] +=
          wdeltax[idir] * (W(er,mark,offset)-W(el,mark,offset));
        if (er < nelem)
          rhs_ls[er][c][idir] +=
            wdeltax[idir] * (W(er,mark,offset)-W(el,mark,offset));
      }
    }
  }
}

void
tk::bndLeastSqConservedVar_P0P1(
  ncomp_t system,
  ncomp_t ncomp,
  ncomp_t offset,
  std::size_t rdof,
  const std::vector< bcconf_t >& bcconfig,
  const inciter::FaceData& fd,
  const Fields& geoFace,
  const Fields& geoElem,
  real t,
  const StateFn& state,
  const Fields& P,
  const Fields& U,
  std::vector< std::vector< std::array< real, 3 > > >& rhs_ls,
  std::size_t nprim )
// *****************************************************************************
//  \brief Compute boundary surface contributions to rhs vector of the
//    least-squares reconstruction of conserved quantities of the PDE system
//! \param[in] system Equation system index
//! \param[in] ncomp Number of scalar components in this PDE system
//! \param[in] offset Offset this PDE system operates from
//! \param[in] rdof Maximum number of reconstructed degrees of freedom
//! \param[in] bcconfig BC configuration vector for multiple side sets
//! \param[in] fd Face connectivity and boundary conditions object
//! \param[in] geoFace Face geometry array
//! \param[in] geoElem Element geometry array
//! \param[in] t Physical time
//! \param[in] state Function to evaluate the left and right solution state at
//!   boundaries
//! \param[in] P Primitive vector to be reconstructed at recent time step
//! \param[in] U Solution vector to be reconstructed at recent time step
//! \param[in,out] rhs_ls RHS reconstruction vector
//! \param[in] nprim This is the number of primitive quantities stored for this
//!   PDE system. This is necessary to extend the state vector to the right
//!   size, so that correct boundary conditions are obtained.
//!   A default is set to 0, so that calling code for systems that do not store
//!   primitive quantities does not need to specify this argument.
//! \details This function computing the boundary face contributions to the rhs
//!   vector for reconstruction, is specialized for conserved quantities. Also
//!   see the version for primitve quantities bndLeastSqPrimitiveVar_P0P1().
// *****************************************************************************
{
  const auto& bface = fd.Bface();
  const auto& esuf = fd.Esuf();

  for (const auto& s : bcconfig) {       // for all bc sidesets
    auto bc = bface.find( std::stoi(s) );// faces for side set
    if (bc != end(bface))
    {
      // Compute boundary face contributions
      for (const auto& f : bc->second)
      {
        Assert( esuf[2*f+1] == -1, "physical boundary element not -1" );

        std::size_t el = static_cast< std::size_t >(esuf[2*f]);

        // arrays for quadrature points
        std::array< real, 3 >
          fc{{ geoFace(f,4,0), geoFace(f,5,0), geoFace(f,6,0) }};
        std::array< real, 3 >
          fn{{ geoFace(f,1,0), geoFace(f,2,0), geoFace(f,3,0) }};

        // Compute the state variables at the left element
        std::vector< real >B(1,1.0);
        auto ul = eval_state( ncomp, offset, rdof, 1, el, U, B );
        auto uprim = eval_state( nprim, offset, rdof, 1, el, P, B );

        // consolidate primitives into state vector
        ul.insert(ul.end(), uprim.begin(), uprim.end());

        Assert( ul.size() == ncomp+nprim, "Incorrect size for "
                "appended state vector" );

        // Compute the state at the face-center using BC
        auto ustate = state( system, ncomp, ul, fc[0], fc[1], fc[2], t, fn );

        std::array< real, 3 > wdeltax{{ fc[0]-geoElem(el,1,0),
                                        fc[1]-geoElem(el,2,0),
                                        fc[2]-geoElem(el,3,0) }};

        for (std::size_t idir=0; idir<3; ++idir)
        {
          // rhs vector
          for (ncomp_t c=0; c<ncomp; ++c)
            rhs_ls[el][c][idir] +=
              wdeltax[idir] * (ustate[1][c]-ustate[0][c]);
        }
      }
    }
  }
}

void
tk::bndLeastSqPrimitiveVar_P0P1(
  ncomp_t system,
  ncomp_t nprim,
  ncomp_t offset,
  std::size_t rdof,
  const std::vector< bcconf_t >& bcconfig,
  const inciter::FaceData& fd,
  const Fields& geoFace,
  const Fields& geoElem,
  real t,
  const StateFn& state,
  const Fields& P,
  const Fields& U,
  std::vector< std::vector< std::array< real, 3 > > >& rhs_ls,
  std::size_t ncomp )
// *****************************************************************************
//  \brief Compute boundary surface contributions to rhs vector of the
//    least-squares reconstruction of primitive quantities of the PDE system
//! \param[in] system Equation system index
//! \param[in] nprim Number of primitive quantities stored for this PDE system
//! \param[in] offset Offset this PDE system operates from
//! \param[in] rdof Maximum number of reconstructed degrees of freedom
//! \param[in] bcconfig BC configuration vector for multiple side sets
//! \param[in] fd Face connectivity and boundary conditions object
//! \param[in] geoFace Face geometry array
//! \param[in] geoElem Element geometry array
//! \param[in] t Physical time
//! \param[in] state Function to evaluate the left and right solution state at
//!   boundaries
//! \param[in] P Primitive vector to be reconstructed at recent time step
//! \param[in] U Conserved vector at recent time step
//! \param[in,out] rhs_ls RHS reconstruction vector
//! \param[in] ncomp This is the number of conserved quantities stored for this
//!   system. This is necessary to extend the state vector to the right size,
//!   so that correct boundary conditions are obtained.
//! \details Since this function computes rhs contributions of boundary faces
//!   for the reconstruction of primitive quantities only, it should not be
//!   called for systems of PDEs that do not store primitive quantities. Also
//!   see the version for conserved quantities bndLeastSqConservedVar_P0P1().
// *****************************************************************************
{
  const auto& bface = fd.Bface();
  const auto& esuf = fd.Esuf();

  Assert( nprim != 0, "Primitive variables reconstruction should not be called "
          "for PDE systems not storing primitive quantities");

  for (const auto& s : bcconfig) {       // for all bc sidesets
    auto bc = bface.find( std::stoi(s) );// faces for side set
    if (bc != end(bface))
    {
      // Compute boundary face contributions
      for (const auto& f : bc->second)
      {
        Assert( esuf[2*f+1] == -1, "physical boundary element not -1" );

        std::size_t el = static_cast< std::size_t >(esuf[2*f]);

        // arrays for quadrature points
        std::array< real, 3 >
          fc{{ geoFace(f,4,0), geoFace(f,5,0), geoFace(f,6,0) }};
        std::array< real, 3 >
          fn{{ geoFace(f,1,0), geoFace(f,2,0), geoFace(f,3,0) }};

        // Compute the state variables at the left element
        std::vector< real >B(1,1.0);
        auto ul = eval_state( nprim, offset, rdof, 1, el, P, B );
        auto ucons = eval_state(ncomp, offset, rdof, 1, el, U, B );

        // consolidate conserved quantities into state vector
        tk::concat< tk::real >(std::move(ul), ucons);

        Assert( ucons.size() == ncomp+nprim, "Incorrect size for "
                "appended state vector" );

        // Compute the state at the face-center using BC
        auto ustate = state( system, ncomp, ucons, fc[0], fc[1], fc[2], t, fn );

        std::array< real, 3 > wdeltax{{ fc[0]-geoElem(el,1,0),
                                        fc[1]-geoElem(el,2,0),
                                        fc[2]-geoElem(el,3,0) }};

        for (std::size_t idir=0; idir<3; ++idir)
        {
          // rhs vector
          for (ncomp_t c=0; c<nprim; ++c)
          {
            auto cp = ustate[0].size()-nprim+c;
            rhs_ls[el][c][idir] +=
              wdeltax[idir] * (ustate[1][cp]-ustate[0][cp]);
          }
        }
      }
    }
  }
}

void
tk::solveLeastSq_P0P1(
  ncomp_t ncomp,
  ncomp_t offset,
  const std::size_t rdof,
  const std::vector< std::array< std::array< real, 3 >, 3 > >& lhs,
  const std::vector< std::vector< std::array< real, 3 > > >& rhs,
  Fields& W )
// *****************************************************************************
//  Solve the 3x3 linear system for least-squares reconstruction
//! \param[in] ncomp Number of scalar components in this PDE system
//! \param[in] offset Offset this PDE system operates from
//! \param[in] rdof Maximum number of reconstructed degrees of freedom
//! \param[in] lhs LHS reconstruction matrix
//! \param[in] rhs RHS reconstruction vector
//! \param[in,out] W Solution vector to be reconstructed at recent time step
//! \details Solves the 3x3 linear system for each element, individually. For
//!   systems that require reconstructions of primitive quantities, this should
//!   be called twice, once with the argument 'W' as U (conserved), and again
//!   with 'W' as P (primitive).
// *****************************************************************************
{
  auto nelem = lhs.size();

  for (std::size_t e=0; e<nelem; ++e)
  {
    for (ncomp_t c=0; c<ncomp; ++c)
    {
      auto mark = c*rdof;

      // solve system using Cramer's rule
      auto ux = tk::cramer( lhs[e], rhs[e][c] );

      W(e,mark+1,offset) = ux[0];
      W(e,mark+2,offset) = ux[1];
      W(e,mark+3,offset) = ux[2];
    }
  }
}

void
tk::recoLeastSqExtStencil(
  std::size_t rdof,
  std::size_t offset,
  std::size_t nielem,
  const std::map< std::size_t, std::vector< std::size_t > >& esup,
  const std::vector< std::size_t >& inpoel,
  const Fields& geoElem,
  Fields& W )
// *****************************************************************************
//  \brief Reconstruct the second-order solution using least-squares approach
//    from an extended stencil involving the node-neighbors
//! \param[in] rdof Maximum number of reconstructed degrees of freedom
//! \param[in] offset Offset this PDE system operates from
//! \param[in] nielem Number of internal elements in this mesh chunk
//! \param[in] esup Elements surrounding points
//! \param[in] inpoel Element-node connectivity
//! \param[in] geoElem Element geometry array
//! \param[in,out] W Solution vector to be reconstructed at recent time step
//! \details A second-order (piecewise linear) solution polynomial is obtained
//!   from the first-order (piecewise constant) FV solutions by using a
//!   least-squares (LS) reconstruction process. This LS reconstruction function
//!   using the nodal-neighbors of a cell, to get an overdetermined system of
//!   equations for the derivatives of the solution. This overdetermined system
//!   is solved in the least-squares sense using the normal equations approach.
// *****************************************************************************
{
  const auto ncomp = W.nprop()/rdof;

  for (std::size_t e=0; e<nielem; ++e)
  {
    // lhs matrix
    std::array< std::array< tk::real, 3 >, 3 >
      lhs_ls( {{ {{0.0, 0.0, 0.0}},
                 {{0.0, 0.0, 0.0}},
                 {{0.0, 0.0, 0.0}} }} );
    // rhs matrix
    std::vector< std::array< tk::real, 3 > >
    rhs_ls( ncomp, {{ 0.0, 0.0, 0.0 }} );

    // loop over all nodes of the element e
    for (std::size_t lp=0; lp<4; ++lp)
    {
      auto p = inpoel[4*e+lp];
      const auto& pesup = cref_find(esup, p);

      // loop over all the elements surrounding this node p
      for (auto er : pesup)
      {
        // centroid distance
        std::array< real, 3 > wdeltax{{ geoElem(er,1,0)-geoElem(e,1,0),
                                        geoElem(er,2,0)-geoElem(e,2,0),
                                        geoElem(er,3,0)-geoElem(e,3,0) }};

        // contribute to lhs matrix
        for (std::size_t idir=0; idir<3; ++idir)
          for (std::size_t jdir=0; jdir<3; ++jdir)
            lhs_ls[idir][jdir] += wdeltax[idir] * wdeltax[jdir];

        // compute rhs matrix
        for (std::size_t c=0; c<ncomp; ++c)
        {
          auto mark = c*rdof;
          for (std::size_t idir=0; idir<3; ++idir)
            rhs_ls[c][idir] +=
              wdeltax[idir] * (W(er,mark,offset)-W(e,mark,offset));
        }
      }
    }

    // solve least-square normal equation system using Cramer's rule
    for (ncomp_t c=0; c<ncomp; ++c)
    {
      auto mark = c*rdof;

      auto ux = tk::cramer( lhs_ls, rhs_ls[c] );

      // Update the P1 dofs with the reconstructioned gradients.
      // Since this reconstruction does not affect the cell-averaged solution,
      // W(e,mark+0,offset) is unchanged.
      W(e,mark+1,offset) = ux[0];
      W(e,mark+2,offset) = ux[1];
      W(e,mark+3,offset) = ux[2];
    }
  }
}

void
tk::transform_P0P1( ncomp_t ncomp,
                    ncomp_t offset,
                    std::size_t rdof,
                    std::size_t nelem,
                    const std::vector< std::size_t >& inpoel,
                    const UnsMesh::Coords& coord,
                    Fields& W )
// *****************************************************************************
//  Transform the reconstructed P1-derivatives to the Dubiner dofs
//! \param[in] ncomp Number of scalar components in this PDE system
//! \param[in] offset Index for equation systems
//! \param[in] rdof Total number of reconstructed dofs
//! \param[in] nelem Total number of elements
//! \param[in] inpoel Element-node connectivity
//! \param[in] coord Array of nodal coordinates
//! \param[in,out] W Second-order reconstructed vector which gets transformed to
//!   the Dubiner reference space
//! \details Since the DG solution (and the primitive quantities) are assumed to
//!   be stored in the Dubiner space, this transformation from Taylor
//!   coefficients to Dubiner coefficients is necessary.
// *****************************************************************************
{
  const auto& cx = coord[0];
  const auto& cy = coord[1];
  const auto& cz = coord[2];

  for (std::size_t e=0; e<nelem; ++e)
  {
    // Extract the element coordinates
    std::array< std::array< real, 3>, 4 > coordel {{
      {{ cx[ inpoel[4*e  ] ], cy[ inpoel[4*e  ] ], cz[ inpoel[4*e  ] ] }},
      {{ cx[ inpoel[4*e+1] ], cy[ inpoel[4*e+1] ], cz[ inpoel[4*e+1] ] }},
      {{ cx[ inpoel[4*e+2] ], cy[ inpoel[4*e+2] ], cz[ inpoel[4*e+2] ] }},
      {{ cx[ inpoel[4*e+3] ], cy[ inpoel[4*e+3] ], cz[ inpoel[4*e+3] ] }}
    }};

    auto jacInv =
      tk::inverseJacobian( coordel[0], coordel[1], coordel[2], coordel[3] );

    // Compute the derivatives of basis function for DG(P1)
    auto dBdx = tk::eval_dBdx_p1( rdof, jacInv );

    for (ncomp_t c=0; c<ncomp; ++c)
    {
      auto mark = c*rdof;

      // solve system using Cramer's rule
      auto ux = tk::cramer( {{ {{dBdx[0][1], dBdx[0][2], dBdx[0][3]}},
                               {{dBdx[1][1], dBdx[1][2], dBdx[1][3]}},
                               {{dBdx[2][1], dBdx[2][2], dBdx[2][3]}} }},
                            {{ W(e,mark+1,offset),
                               W(e,mark+2,offset),
                               W(e,mark+3,offset) }} );

      // replace physical derivatives with transformed dofs
      W(e,mark+1,offset) = ux[0];
      W(e,mark+2,offset) = ux[1];
      W(e,mark+3,offset) = ux[2];
    }
  }
}

void
tk::safeReco( std::size_t offset,
              std::size_t rdof,
              std::size_t nmat,
              std::size_t el,
              int er,
              const Fields& U,
              std::array< std::vector< real >, 2 >& state )
// *****************************************************************************
//  Compute safe reconstructions near material interfaces
//! \param[in] offset Index for equation systems
//! \param[in] rdof Total number of reconstructed dofs
//! \param[in] nmat Total number of material is PDE system
//! \param[in] el Element on the left-side of face
//! \param[in] er Element on the right-side of face
//! \param[in] U Solution vector at recent time-stage
//! \param[in,out] state Second-order reconstructed state, at cell-face, that
//!   is being modified for safety
//! \details When the consistent limiting is applied, there is a possibility
//!   that the material densities and energies violate TVD bounds. This function
//!   enforces the TVD bounds locally
// *****************************************************************************
{
  using inciter::densityIdx;
  using inciter::densityDofIdx;

  // define a lambda for the safe limiting
  auto safeLimit = [&]( std::size_t c, real ul, real ur )
  {
    // find min/max at the face
    auto uMin = std::min(ul, ur);
    auto uMax = std::max(ul, ur);

    // left-state limiting
    if ( (state[0][c] < ul) &&
      (state[0][c] < ur) )
    {
      state[0][c] = uMin;
    }
    else if ( (state[0][c] > ul) &&
      (state[0][c] > ur) )
    {
      state[0][c] = uMax;
    }

    // right-state limiting
    if (er > -1)
    {
      if ( (state[1][c] < ul) &&
        (state[1][c] < ur) )
      {
        state[1][c] = uMin;
      }
      else if ( (state[1][c] > ul) &&
        (state[1][c] > ur) )
      {
        state[1][c] = uMax;
      }
    }
  };

  for (std::size_t k=0; k<nmat; ++k)
  {
    real ul(0.0), ur(0.0);

    // establish left- and right-hand states
    ul = U(el, densityDofIdx(nmat, k, rdof, 0), offset);
    if (er <= -1)
    {
      ur = state[1][densityIdx(nmat, k)];
    }
    else
    {
      auto eR = static_cast< std::size_t >(er);
      ur = U(eR, densityDofIdx(nmat, k, rdof, 0), offset);
    }

    // limit reconstructed density
    safeLimit(densityIdx(nmat,k), ul, ur);

  }
}