[petsc-dev] Periodic meshes with <3 elements per edge?

[Jed Brown]

Jed Brown via petsc-dev <petsc-dev at mcs.anl.gov> writes:

> Matthew Knepley <knepley at gmail.com> writes:
>
>> On Tue, Aug 13, 2019 at 7:35 PM Stefano Zampini <stefano.zampini at gmail.com>
>> wrote:
>>
>>>
>>>
>>> On Aug 14, 2019, at 1:19 AM, Jed Brown via petsc-dev <
>>> petsc-dev at mcs.anl.gov> wrote:
>>>
>>> [Cc: petsc-dev]
>>>
>>> Also, why is our current mode of localized coordinates preferred over
>>> the coordinate DM being non-periodic?  Is the intent to preserve that
>>> for every point in a DM, the point is also valid in the coordinate DM?
>>>
>>> Yes.
>>
>>> Can there be "gaps" in a chart?
>>>
>>> Yes.
>>
>>> I've been digging around in the implementation because there is no
>>> documentation of localized coordinates, but it feels more complicated
>>> than I'd have hoped.
>>>
>>>
>>> A while ago, “localization” of coordinates was supporting only very simple
>>> cases, where periodic points were identified though the ‘maxCell’ parameter
>>> (used to compute  the proper cell coordinates). I think this is the reason
>>> why you need at least 3 cells to support periodicity, since the BoxMesh
>>> constructor uses the maxCell trick.
>>>
>>
>> This was my original conception since it was the only fully automatic way
>> to apply periodicity on an unstructured mesh that
>> was read from a file or generator.
>
> I take it this is in reference to this logic.
>
> PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
> {
>   PetscInt d;
>
>   PetscFunctionBegin;
>   if (!dm->maxCell) {
>     for (d = 0; d < dim; ++d) out[d] = in[d];
>   } else {
>     for (d = 0; d < dim; ++d) {
>       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
>         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
>       } else {
>         out[d] = in[d];
>       }
>     }
>   }
>   PetscFunctionReturn(0);
> }
>
>
> This implies that the mesh be aligned to coordinates (at least in the
> periodic dimension) and that dm->L[d] be set (perhaps by the user?).
> And if the mesh generator specifies periodicity, as with Gmsh below,
> this logic isn't needed.

If we know dm->L[d], we can also know the min and max values in that
dimension.  Then, any time we come to an inverted element, we can
localize such that the Jacobian is positive and coordinates lie in the
bounding box [min,max].

I guess you're concerned about doubly-periodic meshes, such as this ("a"
is the implicit alias for "A"), where the periodic element Dcab has
positive Jacobian since it is a rotation of ABCD.

 a -- b -- a
 |    |    |
 C -- D -- c
 |    |    |
 A -- B -- a

Note that BacD and CDba are correctly flagged as negative Jacobian in
this case, but (without some biasing choice) we can't determine whether
to fix an element given an BACD to become BacD or baCD.  Let's use the
current algorithm for doubly-periodic.

But when there is only one dimension of periodicity, what do you think
of using the bounding box?

 E -- F -- e
 |    |    |
 C -- D -- c
 |    |    |
 A -- B -- a

In this case, we have elements like BacD and know that we can't move B
and D because they would go outside the bounding box.  We need to move
two vertices and the only choices that stay in the bounding box are
x[A]+(L,0) and x[C]+(L,0).

We can even have

 E -- e
 |    |
 C -- c
 |    |
 A -- a

which gives us zero Jacobian (AacC), and we know we have to move two
vertices.  The only positive orientation is AacC because aAcC is tangled
and aACc has negative Jacobian.

Where does this sort of technique fail?