DMPlex: fix bad declaration order in point location routine

DMPlex: complex and C++ fixes to recent point location routines

DMPlex: reduce max iterations for point location.

As I mentioned in a previous commit, there are valid trilinear hexahedra
for which Newton's method will inevitable converge to roots outside the
ref

DMPlex: reduce max iterations for point location.

As I mentioned in a previous commit, there are valid trilinear hexahedra
for which Newton's method will inevitable converge to roots outside the
reference element, even when there is also a root in the reference element.
So, given that we are accepting imperfection from point location (a method
that is not guaranteed to work for all valid element maps), I chose a
maximum number of iterations that seems to be sufficient for even some
hard-case deformed elements.

show more ...

DMPlex: indexing and work array size corrections for point location routines

DMPlex: add debugging info in point location routines

They were useful for demonstrating quadratic convergence, but we don't want them otherwise

DMPlex: fix work array declarations in reference/coordinate transformations

Merge branch 'master' into tisaac/dmplex-feature-real-to-ref

I need updated PetscFE interfaces in what I'm going to do

* master: (603 commits)
parmetis: fix typo in testsuite
Added PETSC_HAVE_A

Merge branch 'master' into tisaac/dmplex-feature-real-to-ref

I need updated PetscFE interfaces in what I'm going to do

* master: (603 commits)
parmetis: fix typo in testsuite
Added PETSC_HAVE_ATTRIBUTEALIGNED flag to makefile
PCSetUp_SVD() used incorrect matrix for matrix size in parallel, also need to skip LAPACK call if matrix size is zero
strumpack: update to latest version 1.0.4
Remove valgrind warning about MPI_Allreduce() with uninitialized data
sed for test output that has atol convergence with quad precision
updated ex1_nest tests
fixed incorrect key value
fix: bad merge a69d34a
test: move runex78 to TESTEXAMPLES_C_NOTSINGLE
build: fix warnings
sys: update test to produce same result - irrespective of debug/opt build - and presence of ~/.petscrc
sys: src/sys/examples/tests/ex30.c -> src/sys/examples/tests/ex32.c to avoid conflict in master
Sed output, since quad precision reduces below atol on these tests
remove automatically generated matrix names from test output files since can be different
Fixes for MatMPIXXXSetPreallocation(), they were not reentrant, since they kept the information used for matrix-vector products
GAMG: fix PCView_GAMG() format - and update output/ex56_0.out
superlu_dist: 64-bit-indices fix https://github.com/xiaoyeli/superlu_dist/issues/4
Add SNESSetDivergenceTolerance / -snes_divergence_tolerance.
fixed unit read in Plex
...

show more ...

DMPlex: implement point inversion for multilinear maps.

This is implemented by a fixed number of Newton iterations. In 2D, I
know that this is a contraction map for a non-degenerate element, and
wi

DMPlex: implement point inversion for multilinear maps.

This is implemented by a fixed number of Newton iterations. In 2D, I
know that this is a contraction map for a non-degenerate element, and
will always converge to the right solution. The same cannot be said in
3D. There are non-degenerate trilinearly mapped hexahedra for which any
Newton iteration, damped or undamped, can converge to a spurious
solution outside of the element, even if there is a solution inside the
element. One could go about a costly procedure of trying to determine
if the target coordinate is inside of the region defined by the mapped
boundary of the reference cell, but this is an expensive procedure. For
now, I am relying on the elements not being so deformed that Newton
doesn't converge. I spent much longer than I should have trying to come
up with a better solution to this problem. I should write the whole
thing up.

WIP: implementation of all point location routines

DMPlex: added complementary DMPlexReferenceToCoordinates() methods

For affine, this just wraps CoordinatesRefToReal.
For tensor, it uses the tensor structure.
In general, the basis space of the PetscDS is used to map points.

WIP: DMPlex: test point inversion.

It seems to work, but the test is flawed for simplicial elements

DMPlex: point location is working for implicit discretizations

Now to do the work for non-implicit (e.g. isoparametric) meshes

show more ...

Merge branch 'stefano_zampini/feature-convert-matnest-matis' into stefano_zampini/feature-pcbddc-saddlepoint

I have to merge it here, since there's an API addition for MatIS in stefano_zampini/featu

Merge branch 'stefano_zampini/feature-convert-matnest-matis' into stefano_zampini/feature-pcbddc-saddlepoint

I have to merge it here, since there's an API addition for MatIS in stefano_zampini/feature-pcbddc-saddlepoint
that prevents me to smoothly merge stefano_zampini/feature-convert-matnest-matis into next.

This will also facilitate merging stefano_zampini/feature-pcbddc-saddlepoint into master

Conflicts:
src/docs/website/documentation/changes/dev.html

show more ...

merged new quadrature behavior

Merge branch 'tisaac/petscfe-matching-quadrature-order'

* tisaac/petscfe-matching-quadrature-order:
Share quadrature instead of sharing quadrature order
SNES tutorials ex12: update regression te

Merge branch 'tisaac/petscfe-matching-quadrature-order'

* tisaac/petscfe-matching-quadrature-order:
Share quadrature instead of sharing quadrature order
SNES tutorials ex12: update regression tests for new quadrature
DMPlex tests ex3: finish updating tests for more quadrature
DMPlex tests ex3: update regression test for discontinuous elements
DMPlex: silence compiler warning for uninitialized variable.
DMPlex tests ex3: update interpolation/convergence tests for higher-order quadrature
DMPlex: use temporary detJ everywhere
fixing Toby's fix
DMPlexComputeGeometryFEM(): clarify what is expected if fe != NULL and correct computations accordingly
PetscFECreateDefault(): increase default quadrature order, make configurable

show more ...

DMPlex: propose DMPlexCoordinatesToReference() interface

DMPlex: silence compiler warning for uninitialized variable.

DMPlex: use temporary detJ everywhere

Merge branch 'tisaac/dmplex-fix-nonaffine-geom' into tisaac/petscfe-matching-quadrature-order

* tisaac/dmplex-fix-nonaffine-geom:
fixing Toby's fix
DMPlexComputeGeometryFEM(): clarify what is ex

Merge branch 'tisaac/dmplex-fix-nonaffine-geom' into tisaac/petscfe-matching-quadrature-order

* tisaac/dmplex-fix-nonaffine-geom:
fixing Toby's fix
DMPlexComputeGeometryFEM(): clarify what is expected if fe != NULL and correct computations accordingly
added checks for AIJ in MIS and GAMG
Bib: Added refs
loop upper bound was one too high
force runex22_4 to produce same output in __float128 precision where it wants to do one more timestep
remove duplicate .bib item, fix typos in .bib items
Fortran wrapper for DMGetType (copy-pasted from src/mat/interface/ftn-custom/zmatregf.c).
Variants of example also need to handle atol convergence with quad precision in test
Add rtol and atol to output of TSView ASCII
Fortran wrapper for MatNullSpaceGetVecs
MATIS: fix compiling error from nightly builds
MATIS: set block size if it is greater than zero
MATIS: added MatTranspose_IS
Quadrature: Fix creation init
Plex: Fix output optimization - We obviously have no test for this logic

show more ...

fixing Toby's fix

Merge branch 'master' into barry/fix-successful-numerical-factorization-after-failure/master

The maint patch for the given problem cannot merge cleanly into next due to changes
in the same code in m

Merge branch 'master' into barry/fix-successful-numerical-factorization-after-failure/master

The maint patch for the given problem cannot merge cleanly into next due to changes
in the same code in master (the introduction of MatFactorGetError()) therefor need
to merge in master and use the new construct before merging to next

show more ...

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 ...

Merge branch 'master' into psanan/pc-telescope-nearnullspace

fixed merge with master

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 ...

DMPlex: added DMPlexComputePointGeometry_Internal()

Seems trivial, but it makes it easier not to have to handle points
differently.

Merge branch 'master' into stefano_zampini/feature-pcbddc-saddlepoint

Conflicts:
src/docs/website/documentation/changes/dev.html
src/ksp/ksp/examples/tutorials/ex59.c
src/ksp/pc/impls/bddc/bddc.c

Merge branch 'master' into stefano_zampini/feature-pcbddc-saddlepoint

Conflicts:
src/docs/website/documentation/changes/dev.html
src/ksp/ksp/examples/tutorials/ex59.c
src/ksp/pc/impls/bddc/bddc.c
src/ksp/pc/impls/bddc/bddcfetidp.c
src/ksp/pc/impls/bddc/bddcprivate.c
src/ksp/pc/impls/bddc/bddcschurs.c
src/mat/impls/is/matis.c
src/mat/impls/is/matis.h

show more ...

Merge branch 'knepley/feature-plex-particle-example'

* knepley/feature-plex-particle-example: (103 commits)
PetscTypes: * Added newly types into fortran definition
DMSwarmMigrate: printf fix for

Merge branch 'knepley/feature-plex-particle-example'

* knepley/feature-plex-particle-example: (103 commits)
PetscTypes: * Added newly types into fortran definition
DMSwarmMigrate: printf fix for __float128
fix c++ comments for c89
fix scanf warnings using PetscInt_FMT
Sys+DM: Fixing 64-bit format strings
Swarm: Fixed 64-bit integer problems
Swarm: Fixed PetscFree() usage
Swarm: Correct format in scanf
Plex: Small cleanup
Plex: Fixes to draw for complex
Plex: Fixed type for GetNeighbors
Swarm: Cleanup of Dave's code - Guy has a kid and think he doesn't have to play by the rules anymore
SNES ex63: Stokes example with Particles
Swarm: Fixed prototypes
Swarm: Added missing PetscFunctionBegin
Swarm: Fixed DMGetLocalVecTor_Swarm()
Swarm: Added DMSwarmCreate/DestroyLocalVectorFromField()
Swarm: Added missing prototype
Plex: Since we are using an SF, ranks that are not found default to 0, not -1 - You need to check the index for the negative value
Plex: Added a new point location type to remove entries which are not found - This was the default, but its not what DMSwarm wanted and its non-intuitive I agree
...

Conflicts:
src/dm/impls/plex/plexcreate.c
src/dm/interface/dm.c

show more ...

Builder+Doc: Fix interpolation option and improve error message