Hi all,
I'm trying to get the coords of the dofs of a DMPlex for a PetscFE
discretization, for orders greater than 1.
I'm struggling to run dm/impls/plex/tutorials/ex8.c
I've got the following error (with option -view_coord) :
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR: DMGetCoordinatesLocalSetUp() has not been called
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.18.2, unknown
[0]PETSC ERROR:
/home/jobic/projet/fe-utils/petsc/fetools/cmake-build-debug/ex_closure_petsc
on a named luke by jobic Mon Dec 12 23:34:37 2022
[0]PETSC ERROR: Configure options
--prefix=/local/lib/petsc/3.18/p3/gcc/nompi_hdf5 --with-mpi=0
--with-debugging=1 --with-blacs=1 --download-zlib,--download-p4est
--download-hdf5=1 --download-triangle=1 --with-single-library=0
--with-large-file-io=1 --with-shared-libraries=0 -CFLAGS=" -g -O0"
-CXXFLAGS=" -g -O0" -FFLAGS=" -g -O0" PETSC_ARCH=nompi_gcc_hdf5
[0]PETSC ERROR: #1 DMGetCoordinatesLocalNoncollective() at
/home/devel/src_linux/petsc-3.18.0/src/dm/interface/dmcoordinates.c:621
[0]PETSC ERROR: #2 DMPlexGetCellCoordinates() at
/home/devel/src_linux/petsc-3.18.0/src/dm/impls/plex/plexgeometry.c:1291
[0]PETSC ERROR: #3 main() at
/home/jobic/projet/fe-utils/petsc/fetools/src/ex_closure_petsc.c:86
[0]PETSC ERROR: PETSc Option Table entries:
[0]PETSC ERROR: -dm_plex_box_faces 2,2
[0]PETSC ERROR: -dm_plex_dim 2
[0]PETSC ERROR: -dm_plex_simplex 0
[0]PETSC ERROR: -petscspace_degree 1
[0]PETSC ERROR: -view_coord
[0]PETSC ERROR: ----------------End of Error Message -------send entire
error message to petsc-ma...@mcs.anl.gov----------
Maybe i've done something wrong ?
Moreover, i don't quite understand the function DMPlexGetLocalOffsets,
and how to use it with DMGetCoordinatesLocal. It seems that
DMGetCoordinatesLocal do not have the coords of the dofs, but only the
nodes defining the geometry.
I've made some small modifications of ex8.c, but i still have an error :
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Invalid argument
[0]PETSC ERROR: Wrong type of object: Parameter # 1
[0]PETSC ERROR: WARNING! There are option(s) set that were not used!
Could be the program crashed before they were used or a spelling
mistake, etc!
[0]PETSC ERROR: Option left: name:-sol value: vtk:sol.vtu
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.18.2, unknown
[0]PETSC ERROR:
/home/jobic/projet/fe-utils/petsc/fetools/cmake-build-debug/view_coords
on a named luke by jobic Mon Dec 12 23:51:05 2022
[0]PETSC ERROR: Configure options
--prefix=/local/lib/petsc/3.18/p3/gcc/nompi_hdf5 --with-mpi=0
--with-debugging=1 --with-blacs=1 --download-zlib,--download-p4est
--download-hdf5=1 --download-triangle=1 --with-single-library=0
--with-large-file-io=1 --with-shared-libraries=0 -CFLAGS=" -g -O0"
-CXXFLAGS=" -g -O0" -FFLAGS=" -g -O0" PETSC_ARCH=nompi_gcc_hdf5
[0]PETSC ERROR: #1 PetscFEGetHeightSubspace() at
/home/devel/src_linux/petsc-3.18.0/src/dm/dt/fe/interface/fe.c:1692
[0]PETSC ERROR: #2 DMPlexGetLocalOffsets() at
/home/devel/src_linux/petsc-3.18.0/src/dm/impls/plex/plexceed.c:98
[0]PETSC ERROR: #3 ViewOffsets() at
/home/jobic/projet/fe-utils/petsc/fetools/src/view_coords.c:28
[0]PETSC ERROR: #4 main() at
/home/jobic/projet/fe-utils/petsc/fetools/src/view_coords.c:99
[0]PETSC ERROR: PETSc Option Table entries:
[0]PETSC ERROR: -heat_petscspace_degree 2
[0]PETSC ERROR: -sol vtk:sol.vtu
[0]PETSC ERROR: ----------------End of Error Message -------send entire
error message to petsc-ma...@mcs.anl.gov----------
Is dm/impls/plex/tutorials/ex8.c a good example for viewing the coords
of the dofs of a DMPlex ?
Thanks,
Yann
static char help[] = "Test the function \n\n";
/* run with -heat_petscspace_degree 2 -sol vtk:sol.vtu */
/* En fait petsc, pour le heat_petscspace_degree 0 garde le maillage ori,
traingulaire 3 sommets */
/* Pour les ordres superieurs, il raffine les elements, donc a 2 fois plus
d'elements par directions */
/* et garde de l'ordre 1 en geometrie
*/
/* Il faut tester avec HDF5 voir si ce format ne permet pas des elements de
plus haut degres */
#include "petscdm.h"
#include <petscdmplex.h>
#include <petscfe.h>
#include <petscmath.h>
#include "petscksp.h"
/* simple gaussian centered at (0.5,0.5) */
static void gaussian(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt
uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar
u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[],
const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar
constants[], PetscScalar uexact[]) {
uexact[0] =
PetscExpReal(-(PetscPowReal(x[0]-0.5,2)+PetscPowReal(x[1]-0.5,2))/0.05);
}
static PetscErrorCode ViewOffsets(DM dm, Vec X)
{
PetscInt num_elem, elem_size, num_comp, num_dof;
PetscInt *elem_restr_offsets;
const PetscScalar *x = NULL;
const char *name;
PetscFunctionBegin;
PetscCall(PetscObjectGetName((PetscObject)dm, &name));
PetscCall(DMPlexGetLocalOffsets(dm, NULL, 0, 0, 0, &num_elem, &elem_size,
&num_comp, &num_dof, &elem_restr_offsets));
PetscCall(PetscPrintf(PETSC_COMM_SELF, "DM %s offsets: num_elem %"
PetscInt_FMT ", size %" PetscInt_FMT ", comp %" PetscInt_FMT ", dof %"
PetscInt_FMT "\n", name, num_elem, elem_size, num_comp, num_dof));
if (X) PetscCall(VecGetArrayRead(X, &x));
for (PetscInt c = 0; c < num_elem; c++) {
PetscCall(PetscIntView(elem_size, &elem_restr_offsets[c * elem_size],
PETSC_VIEWER_STDOUT_SELF));
if (x) {
for (PetscInt i = 0; i < elem_size; i++)
PetscCall(PetscScalarView(num_comp, &x[elem_restr_offsets[c * elem_size + i]],
PETSC_VIEWER_STDERR_SELF));
}
}
if (X) PetscCall(VecRestoreArrayRead(X, &x));
PetscCall(PetscFree(elem_restr_offsets));
PetscFunctionReturn(0);
}
int main(int argc, char **argv) {
DM dm;
Vec u;
PetscFE fe;
PetscQuadrature quad;
PetscInt dim=2,NumComp=1,size;
PetscInt cells[3]={2,2,4};
PetscBool isSimplice=PETSC_FALSE;
void (**exactFields)(PetscInt, PetscInt, PetscInt, const PetscInt[],
const PetscInt[], const PetscScalar[], const PetscScalar[], const
PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const
PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt,
const PetscScalar[], PetscScalar[]);
PetscCall(PetscInitialize(&argc,&argv,(char *)0,help));
PetscCall(DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, isSimplice, cells, NULL,
NULL, NULL, PETSC_TRUE, &dm));
PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject) dm),
dim,NumComp,isSimplice,"heat_",PETSC_DECIDE,&fe));
PetscCall(PetscFEViewFromOptions(fe,NULL,"-fe_field_view"));
PetscCall(DMSetField(dm, 0, NULL, (PetscObject) fe));
PetscCall(PetscFEGetQuadrature(fe,&quad));
PetscCall(PetscQuadratureView(quad,PETSC_VIEWER_STDOUT_WORLD));
PetscCall(PetscFEDestroy(&fe));
PetscCall(DMCreateDS(dm));
PetscMalloc(1,&exactFields);
exactFields[0]=gaussian;
PetscCall(DMCreateGlobalVector(dm,&u));
PetscCall(PetscObjectSetName((PetscObject) u, "temp"));
PetscCall(VecGetSize(u,&size));
PetscPrintf(PETSC_COMM_WORLD, "Size of the global vector : %d\n",size);
{
DM cdm;
PetscSection section;
PetscInt c,cStart,cEnd;
Vec X;
PetscInt cdim;
PetscCall(DMGetLocalSection(dm, §ion));
PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
PetscCall(DMGetCoordinateDim(dm, &cdim));
PetscCall(DMGetCoordinateDM(dm, &cdm));
PetscCall(PetscObjectSetName((PetscObject)cdm, "coords"));
for (c = cStart; c < cEnd; ++c) {
const PetscScalar *array;
PetscScalar *x = NULL;
PetscInt ndof;
PetscBool isDG;
PetscCall(DMPlexGetCellCoordinates(dm, c, &isDG, &ndof, &array, &x));
PetscCheck(ndof % cdim == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "ndof
not divisible by cdim");
PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element #%" PetscInt_FMT "
coordinates\n", c - cStart));
for (PetscInt i = 0; i < ndof; i += cdim) PetscCall(PetscScalarView(cdim,
&x[i], PETSC_VIEWER_STDOUT_SELF));
PetscCall(DMPlexRestoreCellCoordinates(dm, c, &isDG, &ndof, &array, &x));
}
PetscCall(ViewOffsets(dm, NULL));
PetscCall(DMGetCoordinatesLocal(dm, &X));
VecGetSize(X,&size);
PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Taille du vecteur des coordonnes :
%d\n",size));
PetscCall(ViewOffsets(dm, X));
}
PetscCall(DMProjectField(dm,0.0,u,exactFields,INSERT_ALL_VALUES, u));
VecViewFromOptions(u,NULL,"-sol");
PetscCall(DMDestroy(&dm));
PetscCall(VecDestroy(&u));
PetscFree(exactFields);
PetscCall(PetscFinalize());
return 0;
}
