[petsc-users] realCoords for DOFs
Yann Jobic
yann.jobic at univ-amu.fr
Mon Dec 12 17:05:29 CST 2022
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-maint at 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-maint at 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
-------------- next part --------------
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;
}
More information about the petsc-users
mailing list