Dear Matt,
I tried it, but it doesn't seem to work.
Attached is a very small working example illustrating the problem.
I create a DMPlexBox Mesh, periodic in the Y direction. I then scale the Y coordinates with a factor 10, and add 1.0 to it. Both
DMGetCoordinatesLocal and DMGetCellCoordinatesLocal.
Then I evaluate the coordinates with DMPlexGetCellCoordinates. Most of the Y coordinates are correct, but not all of them - for
instance, the minimum Y coordinate is 0.0, and this should be 1.0.
Am I doing something wrong?
Thanks and best regards,
Berend.
On 5/17/23 17:58, Matthew Knepley wrote:
On Wed, May 17, 2023 at 11:20 AM Berend van Wachem <[email protected]
<mailto:[email protected]>> wrote:
Dear Matt,
Is there a way to 'redo' the DMLocalizeCoordinates() ? Or to undo it?
Alternatively, can we make the calling of DMLocalizeCoordinates() in the
DMPlexCreate...() routines optional?
Otherwise, we would have to copy all arrays of coordinates from
DMGetCoordinatesLocal() and DMGetCellCoordinatesLocal() before
scaling them.
I am likely not being clear. I think all you have to do is the following:
DMGetCoordinatesLocal(dm, &xl);
VecScale(xl, scale);
DMSetCoordinatesLocal(dm, xl);
DMGetCellCoordinatesLocal(dm, &xl);
VecScale(xl, scale);
DMSetCellCoordinatesLocal(dm, xl);
Does this not work?
Thanks,
Matt
Best regards, Berend.
On 5/17/23 16:35, Matthew Knepley wrote:
> On Wed, May 17, 2023 at 10:21 AM Berend van Wachem <[email protected]
<mailto:[email protected]>
<mailto:[email protected] <mailto:[email protected]>>> wrote:
>
> Dear Matt,
>
> Thanks for getting back to me so quickly.
>
> If I scale each of the coordinates of the mesh (say, I want to cube
each
> co-ordinate), and I do this for both:
>
> DMGetCoordinatesLocal();
> DMGetCellCoordinatesLocal();
>
> How do I know I am not cubing one coordinate multiple times?
>
>
> Good question. Right now, the only connection between the two sets of
coordinates is DMLocalizeCoordinates(). Since
sometimes
> people want to do non-trivial things to
> coordinates, I prefer not to push in an API for "just" scaling, but I
could be convinced
> the other way.
>
> Thanks,
>
> Matt
>
> Thanks, Berend.
>
> On 5/17/23 16:10, Matthew Knepley wrote:
> > On Wed, May 17, 2023 at 10:02 AM Berend van Wachem
> > <[email protected] <mailto:[email protected]>
<mailto:[email protected]
<mailto:[email protected]>> <mailto:[email protected]
<mailto:[email protected]>
> <mailto:[email protected]
<mailto:[email protected]>>>> wrote:
> >
> > Dear PETSc Team,
> >
> > We are using DMPlex, and we create a mesh using
> >
> > DMPlexCreateBoxMesh (.... );
> >
> > and get a uniform mesh. The mesh is periodic.
> >
> > We typically want to "scale" the coordinates (vertices) of
the mesh,
> > and
> > to achieve this, we call
> >
> > DMGetCoordinatesLocal(dm, &coordinates);
> >
> > and scale the entries in the Vector coordinates appropriately.
> >
> > and then
> >
> > DMSetCoordinatesLocal(dm, coordinates);
> >
> >
> > After this, we localise the coordinates by calling
> >
> > DMLocalizeCoordinates(dm);
> >
> > This worked fine up to PETSc 3.18, but with versions after
this, the
> > coordinates we get from the call
> >
> > DMPlexGetCellCoordinates(dm, CellID, &isDG, &CoordSize,
> > &ArrayCoordinates, &Coordinates);
> >
> > are no longer correct if the mesh is periodic. A number of the
> > coordinates returned from calling DMPlexGetCellCoordinates
are wrong.
> >
> > I think, this is because DMLocalizeCoordinates is now
automatically
> > called within the routine DMPlexCreateBoxMesh.
> >
> > So, my question is: How should we scale the coordinates from
a periodic
> > DMPlex mesh so that they are reflected correctly when calling
both
> > DMGetCoordinatesLocal and DMPlexGetCellCoordinates, with
PETSc versions
> > >= 3.18?
> >
> >
> > I think we might have to add an API function. For now, when you
scale
> > the coordinates,
> > can you scale both copies?
> >
> > DMGetCoordinatesLocal()
> > DMGetCellCoordinatesLocal();
> >
> > and then set them back.
> >
> > Thanks,
> >
> > Matt
> >
> > Many thanks, Berend.
> >
> > --
> > What most experimenters take for granted before they begin their
> > experiments is infinitely more interesting than any results to
which
> > their experiments lead.
> > -- Norbert Wiener
> >
> > https://www.cse.buffalo.edu/~knepley/
<https://www.cse.buffalo.edu/~knepley/>
<https://www.cse.buffalo.edu/~knepley/
<https://www.cse.buffalo.edu/~knepley/>> <http://www.cse.buffalo.edu/~knepley/
<http://www.cse.buffalo.edu/~knepley/>
> <http://www.cse.buffalo.edu/~knepley/
<http://www.cse.buffalo.edu/~knepley/>>>
>
>
>
> --
> What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any
results to
> which their experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
<https://www.cse.buffalo.edu/~knepley/> <http://www.cse.buffalo.edu/~knepley/
<http://www.cse.buffalo.edu/~knepley/>>
--
What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to
which their experiments lead.
-- Norbert Wiener
https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
#include <petsc.h>
#include <petscdmplex.h>
#include <assert.h>
static char help[] = "Vertices Example";
#define MF_MAX(a, b) ((a) > (b) ? (a) : (b))
#define MF_MIN(a, b) ((a) < (b) ? (a) : (b))
int main(int argc, char **argv)
{
DM dm;
PetscErrorCode ierr;
const PetscInt dim = 3;
const PetscBool cellSimplex = PETSC_FALSE;
const PetscBool interpBefore = PETSC_FALSE;
const PetscInt iSize[3] = {5, 5, 5};
const PetscScalar Xmin[3] = {0.0, 0.0, 0.0};
const PetscScalar Xmax[3] = {1.0, 1.0, 1.0};
const DMBoundaryType MeshPeriodicity[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_NONE};
Vec coordinates;
PetscScalar *coordinatesArray;
PetscInt coordVecSize;
PetscInt CellHeight, CellStart, CellEnd;
PetscReal yMinCoord = 1000.0, yMaxCoord = -1000.0;
PetscFunctionBegin;
ierr = PetscInitialize(&argc, &argv, (char *) 0, help);
CHKERRQ(ierr);
ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, cellSimplex, iSize, Xmin, Xmax, MeshPeriodicity, interpBefore, &dm);
CHKERRQ(ierr);
ierr = DMSetFromOptions(dm);
CHKERRQ(ierr);
/* First scale/set the coordinates from DMGetCoordinatesLocal*/
ierr = DMGetCoordinatesLocal(dm, &coordinates);
CHKERRQ(ierr);
ierr = VecGetLocalSize(coordinates, &coordVecSize);
CHKERRQ(ierr);
assert(coordVecSize % 3 == 0);
ierr = VecGetArray(coordinates, &coordinatesArray);
CHKERRQ(ierr);
for (PetscInt i = 0; i < coordVecSize; i += dim)
{
PetscScalar *coord = &coordinatesArray[i];
coord[1] *= 10 + 1.0;
}
ierr = VecRestoreArray(coordinates, &coordinatesArray);
CHKERRQ(ierr);
ierr = DMSetCoordinatesLocal(dm, coordinates);
CHKERRQ(ierr);
/* Secondly, scale/set the coordinates from DMGetCellCoordinatesLocal */
ierr = DMGetCellCoordinatesLocal(dm, &coordinates);
CHKERRQ(ierr);
ierr = VecGetSize(coordinates, &coordVecSize);
CHKERRQ(ierr);
assert(coordVecSize % 3 == 0);
ierr = VecGetArray(coordinates, &coordinatesArray);
CHKERRQ(ierr);
for (PetscInt i = 0; i < coordVecSize; i += dim)
{
PetscScalar *coord = &coordinatesArray[i];
coord[1] *= 10 + 1.0;
}
ierr = VecRestoreArray(coordinates, &coordinatesArray);
CHKERRQ(ierr);
ierr = DMSetCellCoordinatesLocal(dm, coordinates);
CHKERRQ(ierr);
/* Check the vertices per Cell basis: */
ierr = DMPlexGetVTKCellHeight(dm, &CellHeight);
CHKERRQ(ierr);
ierr = DMPlexGetHeightStratum(dm, CellHeight, &CellStart, &CellEnd);
CHKERRQ(ierr);
for (PetscInt LocalCellID = CellStart; LocalCellID < CellEnd; LocalCellID++)
{
PetscScalar *Coordinates;
const PetscScalar *ArrayCoordinates;
PetscInt CoordSize;
PetscBool isDG;
ierr = DMPlexGetCellCoordinates(dm, LocalCellID, &isDG, &CoordSize, &ArrayCoordinates, &Coordinates);
CHKERRQ(ierr);
assert(CoordSize % 3 == 0);
for (PetscInt v = 0; v < CoordSize / 3; v++)
{
yMinCoord = MF_MIN(yMinCoord, Coordinates[3 * v + 1]);
yMaxCoord = MF_MAX(yMaxCoord, Coordinates[3 * v + 1]);
}
ierr = DMPlexRestoreCellCoordinates(dm, LocalCellID, &isDG, &CoordSize, &ArrayCoordinates, &Coordinates);
CHKERRQ(ierr);
}
printf("Y-min coord: %e, Y-min coord: %e\n", yMinCoord, yMaxCoord);
PetscFinalize();
}
