DMPlex: propose DMPlexCoordinatesToReference() interface

Merge branch 'tisaac/dualspace-feature-symmetry'

* tisaac/dualspace-feature-symmetry:
DMPlexTree: fix maxFields bound used in loops
SNES tutorials ex12: update output for dual space changes

PetscViewerVTK: Move DM reference increment into the VTK viewer and out of the DM implementation, as suggested by Jed Brown

Merge branch 'tisaac/dualspace-feature-symmetry' into tisaac/petscfe-matching-quadrature-order

* tisaac/dualspace-feature-symmetry:
DMPlexTree: fix maxFields bound used in loops

DMPlexTree: fix maxFields bound used in loops

DMPlex tests ex3: finish updating tests for more quadrature

DMPlex tests ex3: update regression test for discontinuous elements

I was disconcerted that an interpolation result for p1 discontinuous
elements went from pass to fail when I changed the quadrature

DMPlex tests ex3: update regression test for discontinuous elements

I was disconcerted that an interpolation result for p1 discontinuous
elements went from pass to fail when I changed the quadrature, but what
actually happened was that the *restriction* result went from pass to
fail.

The restriction is computed by applying the transpose of interpolation,
and then scaling the result. The scaling of the result is designed to
be correct for a constant fine vector. So the fine vector can be
written as f = L + C, where C is piecewise constant on each coarse cell,
and L is orthogonal to any vector that is constant on each coarse cell.
Applying the restriction operator gives Rf = RL + RC, where RC is exact.

When we underintegrated to determine the element-difference, we were
using a one-point centroid rule f(). So f(Rf - exact) = f(RL) + f(RC -
exact). Because it is a centroid rule, f(RC-exact) = 0. Because it is
a centroid rule and the nodes are balance around the centroid, fR is a
constant vector on each coarse cell, so f(RL) = 0.

Therefore, the fact that the restriction was exact before was an example
of superconvergence, and I'm comfortable updating the test in a way that
takes that superconvergence away.

show more ...

DMPlex: silence compiler warning for uninitialized variable.

DMPlex tests ex3: update interpolation/convergence tests for higher-order quadrature

Some of the convergence rates for derivatives have decreased, because I
think we were seeing a super-convergence

DMPlex tests ex3: update interpolation/convergence tests for higher-order quadrature

Some of the convergence rates for derivatives have decreased, because I
think we were seeing a super-convergence effect from under integration.

I didn't do all of the tests yet: a discontinuous test switched from
pass to FAIL, and I don't trust it yet.

show more ...

DMPlex: use temporary detJ everywhere

fixing Toby's fix

DMPlexComputeGeometryFEM(): clarify what is expected if fe != NULL and correct computations accordingly

It makes the most sense to ask for the image of each quadrature point, and the the J, invJ and

DMPlexComputeGeometryFEM(): clarify what is expected if fe != NULL and correct computations accordingly

It makes the most sense to ask for the image of each quadrature point, and the the J, invJ and detJ of the transformation at those same points.

show more ...

Plex ex6: Updated output for fix

The bottom orientation was incorrect. For order 3, the orientation was listed as:
v1 v2
+---------------+
| |
| 8 10 | top-to-bottom
|

The bottom orientation was incorrect. For order 3, the orientation was listed as:
v1 v2
+---------------+
| |
| 8 10 | top-to-bottom
| | enumeration
| |
| 9 11 |
| |
+---------------+
v0 v3
However, it should have been bottom to top
v1 v2
+---------------+
| |
| 9 11 | bottom-to-top
| | enumeration
| |
| 8 10 |
| |
+---------------+
v0 v3
This resulted in an incorrect permutation. From plex.c:DMPlexCreateSpectralClosurePermutation(), the perm had values:
perm[5] = 9
perm[6] = 11
perm[9] = 8
perm[10] = 10
They should have been
perm[5] = 8
perm[6] = 10
perm[9] = 9
perm[10] = 11
The fixed variable: o = ofb+n*(k-1)+(k-2)-i; was changed to o = ofb+n*(k-1)+i;, reversing the top-to-bottom to bottom-to-top. This continues to work for higher orders as well.

show more ...

DMPlexTree: compute constraint matrix entries in blocks

For high order matrices, calling MatGetValue() and MatSetValue() was
costly, so the critical section of code instead uses one
MatDenseGet/Rest

DMPlexTree: compute constraint matrix entries in blocks

For high order matrices, calling MatGetValue() and MatSetValue() was
costly, so the critical section of code instead uses one
MatDenseGet/RestoreArray() and MatSetValues().

show more ...

Merge branch 'tisaac/dualspace-feature-symmetry'

In finite element assembly, to compute a residual or assemble a Jacobian
on an element, the dual space functionals of a reference finite element are

Merge branch 'tisaac/dualspace-feature-symmetry'

In finite element assembly, to compute a residual or assemble a Jacobian
on an element, the dual space functionals of a reference finite element are related
(i.e., pushed forward) to the functionals that make of the degrees of
freedom of the discretized system.

In doing so, we have to draw a correspondence between the boundary
points of the reference element (its edges and vertices) and the
boundary points of each cell in the mesh. Sometimes, this
correspondence is not perfect, because the orientations of the boundary
points of a cell may be different than the orientations of the boundary
points of the reference element. For example, an edge that runs
clockwise around the reference element may map onto an edge that runs
counter-clockwise around the real cell.

For some finite elements (such as the most common P1 Lagrange finite
elements) this doesn't matter. But for others---high-order elements,
H-div and H-curl conforming elements---a difference in orientation
affects the mapping between reference functionals and degrees of
freedom.

The simplest example is a high-order nodal Lagrange finite element in
2D: the degrees of freedom on an edge of the reference element are
numbered left-to-right, but if the orientation of the real edge is
reverse, then we must read/write degrees of freedom right-to-left.

Things get more complicated in 3D: a quadrilateral face between
hexahedra may have any of 8 different orientations (that make of the
dihedral symmetries), and each one corresponds to a different reordering
of the 2D grid of nodes supported on that quadrilateral.

Things also get more complicated for H-div elements. Functionals on the
interfaces between cells represent the flux from one cell to another,
and thus have direction. Reversing the orientation of an edge changes
the sign of the mapping between the reference functional and the degree
of freedom.

To assemble generic finite elements, we thus need to accommodate
transforms between reference functionals and degrees of freedom that
encompass (1) permutations, and (2) scalar multipliers.

This branch accomplishes this in two steps:

- We add PetscDualSpaceGetSymmetries() to the interface, to allow a dual
space to describe how symmetries of the mesh points affect
referent-to-real maps of functionals.

- We add PetscSectionSym, an object that encapsulates symmetries of the
degrees of freedom whose layout is described by a PetscSection.

When a section is created from a DM that has a finite element
discretization in its PetscDS, it will automatically construct the
appropriate PetscSectionSym when it creates a PetscSection to describe a
discretized function space.

If a PetsFE finite element discretization is not used, the user can
create a PetscSectionSym for their degrees of freedom for themselves.
We give an example of this for spectal elements in
src/dm/impls/plex/examples/tutoerials/ex6.c.

* tisaac/dualspace-feature-symmetry: (22 commits)
cast memzeros to void to avoid compiler complaints about const qualifiers
DMPlex: fix stray PetscScalar to PetscReal
DMPlex: fix stray PetscScalar to PetscReal
Plex tutorials ex6: use symmetries in SEM example
Plex tests ex3: added high order tests
Plex tests ex3: add matrix-free near null space test
DMPlex: fix PetscInt %D PetscPrintf character
DMPlexTree: use symmetries when computing constraint matrices
DMPlex: use symmetries in DMPlexGetIndicesPoint{Fields}_Internal() and DMPlexAnchorsModifyMat()
DMPlex: use PetscSectionGet/RestorePointSyms() in DMPlexVecGet/SetClosure()
DMPlex: added DMPlexGet/RestoreCompressedClosure()
DMPlex: added DMPlexGetPointDualSpaceFEM()
DMPlex: set PetscSectionSym in DMPlexCreateSectionInitial()
DMLabel: added PETSCSECTIONSYMLABEL
PetscSection: add PetscSectionSym
PetscDualSpace: added PetscDualSpaceGetSymmetries()
PetscDualSpace_Lagrange: setup recursively on dimension
PetscDualSpace_Lagrange: order functions now go after tuple used
PetscDualSpace_Lagrange: use lexicographic order for nodes
DMPlex: added DMPlexComputePointGeometry_Internal()
...

show more ...

cast memzeros to void to avoid compiler complaints about const qualifiers

Merge branch 'barry/fix-xml-logging'

DMPlex: fix stray PetscScalar to PetscReal

Plex tutorials ex6: use symmetries in SEM example

Other usage thus far has determined the symmetry from a dual space. In
this example we are demonstrating SEM node ordering *without* relying on
Pet

Plex tutorials ex6: use symmetries in SEM example

Other usage thus far has determined the symmetry from a dual space. In
this example we are demonstrating SEM node ordering *without* relying on
PetscDS/PetscFE/PetscDualSpace, etc. In this case, the user has to
describe the symmetry of their dofs for themselves. So ex6 now
demonstrates how to describe the symmetries for tensor products grids of
degrees of freedom.

show more ...

Plex tests ex3: added high order tests

High order (in 3D) relies on the correcy symmetry permutations, so this
is a test of the recent symmetry infrastructure.

Plex tests ex3: add matrix-free near null space test

The existing tests stressed DMPlexMatSetClosure(); the additional tests
stress DMPlexVecGet/SetClosure()

DMPlex: fix PetscInt %D PetscPrintf character

DMPlexTree: use symmetries when computing constraint matrices

DMPlex: use symmetries in DMPlexGetIndicesPoint{Fields}_Internal() and DMPlexAnchorsModifyMat()

The permutations are passed now instead of the orientations. For
DMPlexGetIndicesPointFields_Internal

DMPlex: use symmetries in DMPlexGetIndicesPoint{Fields}_Internal() and DMPlexAnchorsModifyMat()

The permutations are passed now instead of the orientations. For
DMPlexGetIndicesPointFields_Internal(), we pass permutations for every
point (assuming PetscSectionGetFieldPointSyms() was called for multiple
points, such as for a closure set) of every field, so we have to pass an
additional parameter indicating the offset of the current point in the
set of points.

show more ...