On Tue, Jul 5, 2016 at 4:17 AM, Morten Nobel-Jørgensen <m...@dtu.dk> wrote:
> Hi all,
>
> I hope someone can help me with the following:
>
> I’m having some problems when exporting a distributed DMPlex – the cells
> (+cell types) seems to be duplicated.
>
> When I’m running the code on a non-distributed system it works as
> expected, but when I run it on multiple processors (2 in my case) the
> output is invalid.
>
> I have attached a simple example and the output for np=1 and np=2.
>
The problem here is VTK output with overlapped meshes. If you change to
overlap = 0, it works as expected. I never fixed the VTK output
for this, but the HDF5 output works correctly. I will put this on the list
of things to do. I am attaching your code with some cleanup from me,
including assigning values in parallel.
Thanks,
Matt
> Abbreviated the code essentially does the following:
> '
>
> PetscInt dim = 3;
> PetscInt cells[] = {1, 1, 2};
> PetscInt overlap = 1;
> PetscInitialize(&argc, &argv, NULL, help);
> DMPlexCreateHexBoxMesh(PETSC_COMM_WORLD, dim, cells, DM_BOUNDARY_NONE,
> DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, &dm);
> DMPlexDistribute(dm, overlap, NULL, &dist);
> dm = dist;
> SetupDOFs(dm);
> Vec V;
> DMCreateGlobalVector(dm, &V);
> AssignSomeValues(V);
> PetscViewer viewer;
> const char* fn = "output.vtk";
> PetscViewerVTKOpen(PETSC_COMM_WORLD,fn,FILE_MODE_WRITE,&viewer);
> VecView(V,viewer);
> PetscViewerDestroy(&viewer);
>
>
> Kind regards,
> Morten
>
--
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
static char help[] = "";
#include <petscdmplex.h>
const int dof = 1;
#undef __FUNCT__
#define __FUNCT__ "SetupDOFs"
/* Assign dofs (3 for each node) */
PetscErrorCode SetupDOFs(DM dm) {
PetscSection s;
PetscInt pStart, pEnd, cStart, cEnd, c, fStart, fEnd, f, vStart, vEnd, v, eStart, eEnd, e;
PetscErrorCode ierr;
PetscFunctionBeginUser;
ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &s);CHKERRQ(ierr);
ierr = PetscSectionSetNumFields(s, 1);CHKERRQ(ierr);
ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr);
ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr);
ierr = DMPlexGetHeightStratum(dm, 2, &eStart, &eEnd);CHKERRQ(ierr);
ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
ierr = PetscSectionSetChart(s, pStart, pEnd);CHKERRQ(ierr);
for (v = vStart; v < vEnd; ++v) {
ierr = PetscSectionSetDof(s, v, dof);CHKERRQ(ierr);
ierr = PetscSectionSetFieldDof(s, v, 0, dof);CHKERRQ(ierr);
}
ierr = PetscSectionSetUp(s);CHKERRQ(ierr);
ierr = DMSetDefaultSection(dm, s);CHKERRQ(ierr);
ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "AssignSomeValues"
PetscErrorCode AssignSomeValues(Vec V){
PetscInt vectorSize, rStart, rEnd, r;
PetscErrorCode ierr;
PetscFunctionBeginUser;
ierr = VecGetLocalSize(V, &vectorSize);CHKERRQ(ierr);
ierr = VecGetOwnershipRange(V, &rStart, &rEnd);CHKERRQ(ierr);
ierr = PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Vector size %d\n",vectorSize);CHKERRQ(ierr);
ierr = PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);CHKERRQ(ierr);
for (r = rStart; r < rEnd; ++r) {
const PetscScalar value = r;
ierr = VecSetValues(V, 1, &r, &value, INSERT_VALUES);CHKERRQ(ierr);
}
ierr = VecAssemblyBegin(V);CHKERRQ(ierr);
ierr = VecAssemblyEnd(V);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc, char **argv) {
DM dm, dist;
Mat mat;
Vec V;
PetscInt dim = 3;
PetscInt cells[3] = {1, 1, 2};
PetscInt overlap = 0;
PetscMPIInt rank;
PetscViewer viewer;
const char *fn = "output.vtk";
PetscErrorCode ierr;
ierr = PetscInitialize(&argc, &argv, NULL, help);CHKERRQ(ierr);
ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);
ierr = DMPlexCreateHexBoxMesh(PETSC_COMM_WORLD, dim, cells, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, &dm); CHKERRQ(ierr);
ierr = DMPlexDistribute(dm, overlap, NULL, &dist); CHKERRQ(ierr);
if (dist) {
ierr = DMDestroy(&dm);CHKERRQ(ierr);
dm = dist;
}
ierr = DMViewFromOptions(dm, NULL, "-dm_view"); CHKERRQ(ierr);
ierr = SetupDOFs(dm);CHKERRQ(ierr);
ierr = DMCreateGlobalVector(dm, &V);CHKERRQ(ierr);
ierr = AssignSomeValues(V);CHKERRQ(ierr);
ierr = PetscViewerVTKOpen(PETSC_COMM_WORLD,fn,FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
ierr = VecView(V,viewer);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
ierr = VecDestroy(&V);CHKERRQ(ierr);
ierr = DMDestroy(&dm);CHKERRQ(ierr);
PetscFinalize();
return 0;
}
