High order `FEM` visualization

Nobody has claimed this yet.

Describe the bug

High order visualization seems to be incomplete / buggy (wrong number of local dofs).

Ideally, I would like to get the solution per element instead of per local dof, but AFAIK this seems un-implemented.

To Reproduce
Consider the mwe consisting of:

a reduced .xml using an acoustic solver on a bilayer mesh, using Q2 spectral elements (hence 3 points per hexahedron edge);

xml

<?xml version="1.0" ?>

<Problem>
  <Solvers>
    <AcousticSEM
      name="acousticSolver"
      discretization="FE1"
      targetRegions="{ Fluid1, Fluid2 }"
      sourceCoordinates="{ { 1500, 20, 2900 } }"
      timeSourceFrequency="10"
      receiverCoordinates="{ { 3752.29351, 20, 2933.3333 } }"/>
  </Solvers>

  <NumericalMethods>
    <FiniteElements>
      <FiniteElementSpace
        name="FE1"
        order="2"
        formulation="SEM"/>
    </FiniteElements>
  </NumericalMethods>

  <ElementRegions>
    <CellElementRegion
      name="Fluid1"
      cellBlocks="{ 1_hexahedra }"
      materialList="{ nullModel }"/>
    <CellElementRegion
      name="Fluid2"
      cellBlocks="{ 2_hexahedra }"
      materialList="{ nullModel }"/>
  </ElementRegions>

  <Mesh>
    <VTKMesh
      name="Bilayer"
      fieldsToImport="{ mediumVelocity }"
      fieldNamesInGEOSX="{ mediumVelocity }"
      file="high_order_bug.vtu"
      useGlobalIds="0"/>
  </Mesh> 

  <Outputs>
    <VTK
      name="vtkOutput"
      writeFEMFaces="1"/>
  </Outputs>

  <Events maxTime="1e-1">
    <!-- control the timestepping here with forceDt -->
    <PeriodicEvent
      name="fluidSolverApplications"
      forceDt="1e-2"
      target="/Solvers/acousticSolver"/>
    <!-- output -->
    <PeriodicEvent
      name="outputs"
      timeFrequency="5e-2"
      targetExactTimestep="0"
      target="/Outputs/vtkOutput">
    </PeriodicEvent>
  </Events>

</Problem>

a .vtu input mesh with nx=3, ny=1, nz=2 elements.

Input the .xml to GEOS, and open the output .pvd file in paraview.

Expected behavior

I expect to see the correct number of points (px = 3 * 2 + 1) * (py = 1 * 2 + 1) * (pz = 2 * 2 + 1) = 105 points, and 6 * 8 = 48 subcells.
Also, the ordering seems wrong: applying a PointDatatoCellData filter renders incorrect results (likely a connectivity issue).

Screenshots

bug-q2

On this screenshot:

the number of points in a xy plane seems ok except for the last z layer.
the number of points in the z direction is incorrect (3 instead of 5).

Start by running the provided XML with high-order Q2 spectral elements and the linked high_order_bug.vtu mesh, then inspect the generated .pvd output in ParaView. Compare the output against the expected 105 points, 48 subcells, and correct connectivity when applying PointDataToCellData; done means the high-order visualization matches those expectations.

Written by the indexing model from the issue text.

Description

Contributor guide

First steps

Research direction

Assessment