dm_view of high-order geometry/solution
Dear all, I am playing with dm/impls/plex/tests/ex33.c and know how to set high-order geometry. Is it possible to output the final mesh to vtu, e.g. annulus example? Thanks! Junming
On Mon, Jun 12, 2023 at 4:13 AM Duan Junming via petsc-users < [email protected]> wrote:
Dear all,
I am playing with dm/impls/plex/tests/ex33.c and know how to set high-order geometry.
Is it possible to output the final mesh to vtu, e.g. annulus example?
The problem is that VTK has no nice, standard way to talk about higher order polynomials. Thus, the best way to do this is to call DMRefine() once or twice and output that mesh. Thanks, Matt
Thanks!
Junming
-- 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/>
Dear Matt, Thank you for the reply. I have a more specific question about the spectral element example. Do you have any suggestions that how to write all the nodes in each cell to .vtu? Thanks! Junming ________________________________ From: [email protected] <[email protected]> Sent: Monday, June 12, 2023 1:34:24 PM To: Duan Junming Cc: [email protected] Subject: Re: [petsc-users] dm_view of high-order geometry/solution On Mon, Jun 12, 2023 at 4:13 AM Duan Junming via petsc-users <[email protected]<mailto:[email protected]>> wrote: Dear all, I am playing with dm/impls/plex/tests/ex33.c and know how to set high-order geometry. Is it possible to output the final mesh to vtu, e.g. annulus example? The problem is that VTK has no nice, standard way to talk about higher order polynomials. Thus, the best way to do this is to call DMRefine() once or twice and output that mesh. Thanks, Matt Thanks! Junming -- 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/>
On Mon, Jun 12, 2023 at 6:01 AM Duan Junming <[email protected]> wrote:
Dear Matt,
Thank you for the reply. I have a more specific question about the spectral element example. Do you have any suggestions that how to write all the nodes in each cell to .vtu?
It is the same procedure. VTU is not a great format for this. It wants everything at first order. Thanks, Matt
Thanks!
Junming ------------------------------ *From:* [email protected] <[email protected]> *Sent:* Monday, June 12, 2023 1:34:24 PM *To:* Duan Junming *Cc:* [email protected] *Subject:* Re: [petsc-users] dm_view of high-order geometry/solution
On Mon, Jun 12, 2023 at 4:13 AM Duan Junming via petsc-users < [email protected]> wrote:
Dear all,
I am playing with dm/impls/plex/tests/ex33.c and know how to set high-order geometry.
Is it possible to output the final mesh to vtu, e.g. annulus example?
The problem is that VTK has no nice, standard way to talk about higher order polynomials. Thus, the best way to do this is to call DMRefine() once or twice and output that mesh.
Thanks,
Matt
Thanks!
Junming
-- 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/>
-- 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/>
Matthew Knepley <[email protected]> writes:
On Mon, Jun 12, 2023 at 6:01 AM Duan Junming <[email protected]> wrote:
Dear Matt,
Thank you for the reply. I have a more specific question about the spectral element example. Do you have any suggestions that how to write all the nodes in each cell to .vtu?
It is the same procedure. VTU is not a great format for this. It wants everything at first order.
I would recommend configuring with --download-cgns and running with -dm_view cgns:output.cgns. This format has efficient parallel IO and curved elements work for moderate order in Paraview.
Dear Jed, Thank you for the suggestion. When I run tests/ex33.c with ./ex33 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -mesh_transform annulus -dm_coord_space 0 -dm_coord_petscspace_degree 3 -dm_refine 1 -dm_view cgns:test.cgns and load it using Paraview, the mesh is still with straight lines. Should I modify the code to make cgns work? Or any other examples for me to start? Thanks! Junming ________________________________ From: Jed Brown <[email protected]> Sent: Monday, June 12, 2023 2:52:48 PM To: Matthew Knepley; Duan Junming Cc: [email protected] Subject: Re: [petsc-users] dm_view of high-order geometry/solution Matthew Knepley <[email protected]> writes:
On Mon, Jun 12, 2023 at 6:01 AM Duan Junming <[email protected]> wrote:
Dear Matt,
Thank you for the reply. I have a more specific question about the spectral element example. Do you have any suggestions that how to write all the nodes in each cell to .vtu?
It is the same procedure. VTU is not a great format for this. It wants everything at first order.
I would recommend configuring with --download-cgns and running with -dm_view cgns:output.cgns. This format has efficient parallel IO and curved elements work for moderate order in Paraview.
Duan Junming <[email protected]> writes:
Dear Jed,
Thank you for the suggestion.
When I run tests/ex33.c with
./ex33 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -mesh_transform annulus -dm_coord_space 0 -dm_coord_petscspace_degree 3 -dm_refine 1 -dm_view cgns:test.cgns
and load it using Paraview,
the mesh is still with straight lines.
Ah, the viewer is keyed on the field (since the CGNS as supported by Paraview specifies coordinates and fields in the same space). That doesn't exist in your case. If you apply this patch and add `-petscspace_degre 3` to your command, you'll see that high order information is present. Paraview doesn't render as curves in all views, but it has the data. diff --git i/src/dm/impls/plex/tests/ex33.c w/src/dm/impls/plex/tests/ex33.c index 803095bc082..590facfa4f4 100644 --- i/src/dm/impls/plex/tests/ex33.c +++ w/src/dm/impls/plex/tests/ex33.c @@ -198,7 +198,6 @@ PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm) default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Unknown mesh transform %d", ctx->meshTransform); } - PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); PetscFunctionReturn(PETSC_SUCCESS); } @@ -227,6 +226,7 @@ static PetscErrorCode CreateDiscretization(DM dm, AppCtx *ctx) PetscCall(DMCreateDS(dm)); PetscCall(DMGetDS(dm, &ds)); PetscCall(PetscDSSetObjective(ds, 0, volume)); + PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); PetscFunctionReturn(PETSC_SUCCESS); } I can update the viewer to handle the degenerate case of no field (all my models have fields).
And here's an MR to do what you want without any code/arg changes. https://gitlab.com/petsc/petsc/-/merge_requests/6588 Jed Brown <[email protected]> writes:
Duan Junming <[email protected]> writes:
Dear Jed,
Thank you for the suggestion.
When I run tests/ex33.c with
./ex33 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -mesh_transform annulus -dm_coord_space 0 -dm_coord_petscspace_degree 3 -dm_refine 1 -dm_view cgns:test.cgns
and load it using Paraview,
the mesh is still with straight lines.
Ah, the viewer is keyed on the field (since the CGNS as supported by Paraview specifies coordinates and fields in the same space). That doesn't exist in your case. If you apply this patch and add `-petscspace_degre 3` to your command, you'll see that high order information is present. Paraview doesn't render as curves in all views, but it has the data.
diff --git i/src/dm/impls/plex/tests/ex33.c w/src/dm/impls/plex/tests/ex33.c index 803095bc082..590facfa4f4 100644 --- i/src/dm/impls/plex/tests/ex33.c +++ w/src/dm/impls/plex/tests/ex33.c @@ -198,7 +198,6 @@ PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm) default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Unknown mesh transform %d", ctx->meshTransform); } - PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); PetscFunctionReturn(PETSC_SUCCESS); }
@@ -227,6 +226,7 @@ static PetscErrorCode CreateDiscretization(DM dm, AppCtx *ctx) PetscCall(DMCreateDS(dm)); PetscCall(DMGetDS(dm, &ds)); PetscCall(PetscDSSetObjective(ds, 0, volume)); + PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); PetscFunctionReturn(PETSC_SUCCESS); }
I can update the viewer to handle the degenerate case of no field (all my models have fields).
Dear Jed, Thank you for your help! Now I moved the line using "DMViewFromOptions" after the function "PetscDSSetObjective", and it works for "-dm_coord_petscspace_degree 3 -petscspace_degree 3". But when I tried degree 4: ./ex33 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -mesh_transform annulus -dm_coord_space 0 -dm_coord_petscspace_degree 4 -petscspace_degree 4 -dm_refine 1 -dm_view cgns:test.cgns Paraview gives an empty render. Using degree 5: ./ex33 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -mesh_transform annulus -dm_coord_space 0 -dm_coord_petscspace_degree 5 -petscspace_degree 5 -dm_refine 1 -dm_view cgns:test.cgns it reports: [0]PETSC ERROR: --------------------- Error Message -------------------------------------------------------------- [0]PETSC ERROR: No support for this operation for this object type [0]PETSC ERROR: Cell type quadrilateral with closure size 36 [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting. [0]PETSC ERROR: Petsc Release Version 3.19.2, unknown [0]PETSC ERROR: ./ex33 on a arch-darwin-c-debug named JunmingMacBook-Pro.local by Junming Mon Jun 12 20:23:04 2023 [0]PETSC ERROR: Configure options --download-cgns --download-hdf5 --download-openmpi --download-triangle --with-fc=0 PETSC_ARCH=arch-darwin-c-debug --download-cgns [0]PETSC ERROR: #1 DMPlexCGNSGetPermutation_Internal() at /Users/Junming/Packages/petsc/src/dm/impls/plex/cgns/plexcgns2.c:533 [0]PETSC ERROR: #2 DMView_PlexCGNS() at /Users/Junming/Packages/petsc/src/dm/impls/plex/cgns/plexcgns2.c:769 [0]PETSC ERROR: #3 DMView_Plex() at /Users/Junming/Packages/petsc/src/dm/impls/plex/plex.c:1801 [0]PETSC ERROR: #4 DMView() at /Users/Junming/Packages/petsc/src/dm/interface/dm.c:996 [0]PETSC ERROR: #5 PetscObjectView() at /Users/Junming/Packages/petsc/src/sys/objects/destroy.c:78 [0]PETSC ERROR: #6 PetscObjectViewFromOptions() at /Users/Junming/Packages/petsc/src/sys/objects/destroy.c:128 [0]PETSC ERROR: #7 DMViewFromOptions() at /Users/Junming/Packages/petsc/src/dm/interface/dm.c:940 [0]PETSC ERROR: #8 CreateDiscretization() at ex33.c:232 [0]PETSC ERROR: #9 main() at ex33.c:263 [0]PETSC ERROR: PETSc Option Table entries: [0]PETSC ERROR: -dm_coord_petscspace_degree 5 (source: command line) [0]PETSC ERROR: -dm_coord_space 0 (source: command line) [0]PETSC ERROR: -dm_plex_box_faces 1,1 (source: command line) [0]PETSC ERROR: -dm_plex_simplex 0 (source: command line) [0]PETSC ERROR: -dm_refine 1 (source: command line) [0]PETSC ERROR: -dm_view cgns:test.cgns (source: command line) [0]PETSC ERROR: -mesh_transform annulus (source: command line) [0]PETSC ERROR: -petscspace_degree 5 (source: command line) [0]PETSC ERROR: ----------------End of Error Message -------send entire error message to [email protected] -------------------------------------------------------------------------- MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_SELF with errorcode 56. Does cgns work for degree >= 4? Junming ________________________________ From: Jed Brown <[email protected]> Sent: Monday, June 12, 2023 19:07 To: Duan Junming; Matthew Knepley Cc: [email protected] Subject: Re: [petsc-users] dm_view of high-order geometry/solution And here's an MR to do what you want without any code/arg changes. https://gitlab.com/petsc/petsc/-/merge_requests/6588 Jed Brown <[email protected]> writes:
Duan Junming <[email protected]> writes:
Dear Jed,
Thank you for the suggestion.
When I run tests/ex33.c with
./ex33 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -mesh_transform annulus -dm_coord_space 0 -dm_coord_petscspace_degree 3 -dm_refine 1 -dm_view cgns:test.cgns
and load it using Paraview,
the mesh is still with straight lines.
Ah, the viewer is keyed on the field (since the CGNS as supported by Paraview specifies coordinates and fields in the same space). That doesn't exist in your case. If you apply this patch and add `-petscspace_degre 3` to your command, you'll see that high order information is present. Paraview doesn't render as curves in all views, but it has the data.
diff --git i/src/dm/impls/plex/tests/ex33.c w/src/dm/impls/plex/tests/ex33.c index 803095bc082..590facfa4f4 100644 --- i/src/dm/impls/plex/tests/ex33.c +++ w/src/dm/impls/plex/tests/ex33.c @@ -198,7 +198,6 @@ PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm) default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Unknown mesh transform %d", ctx->meshTransform); } - PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); PetscFunctionReturn(PETSC_SUCCESS); }
@@ -227,6 +226,7 @@ static PetscErrorCode CreateDiscretization(DM dm, AppCtx *ctx) PetscCall(DMCreateDS(dm)); PetscCall(DMGetDS(dm, &ds)); PetscCall(PetscDSSetObjective(ds, 0, volume)); + PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); PetscFunctionReturn(PETSC_SUCCESS); }
I can update the viewer to handle the degenerate case of no field (all my models have fields).
CGNS supports fourth order and it's coded in PETSc, but Paraview hasn't implemented reading it yet. I think it would not be much work for someone (maybe you) to add it to Paraview. I have lots of applications on cubics, but not much beyond that so it hasn't risen to top priority for me. There is an accepted extension, but the CGNS implementation is still in a branch. https://cgns.github.io/ProposedExtensions/CPEX0045_HighOrder_v2.pdf https://github.com/CGNS/CGNS/tree/CPEX45_high_order There have been recent (announced in a blog post, yet still undocumented) extensions to the VTU format that would support this, but the format is so bad for parallel IO and time series that I haven't been motivated to extend the PETSc writer. Of course we would welcome contributions. Duan Junming <[email protected]> writes:
Dear Jed,
Thank you for your help!
Now I moved the line using "DMViewFromOptions" after the function "PetscDSSetObjective",
and it works for "-dm_coord_petscspace_degree 3 -petscspace_degree 3".
But when I tried degree 4:
./ex33 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -mesh_transform annulus -dm_coord_space 0 -dm_coord_petscspace_degree 4 -petscspace_degree 4 -dm_refine 1 -dm_view cgns:test.cgns
Paraview gives an empty render.
Using degree 5:
./ex33 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -mesh_transform annulus -dm_coord_space 0 -dm_coord_petscspace_degree 5 -petscspace_degree 5 -dm_refine 1 -dm_view cgns:test.cgns
it reports:
[0]PETSC ERROR: --------------------- Error Message -------------------------------------------------------------- [0]PETSC ERROR: No support for this operation for this object type [0]PETSC ERROR: Cell type quadrilateral with closure size 36 [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting. [0]PETSC ERROR: Petsc Release Version 3.19.2, unknown [0]PETSC ERROR: ./ex33 on a arch-darwin-c-debug named JunmingMacBook-Pro.local by Junming Mon Jun 12 20:23:04 2023 [0]PETSC ERROR: Configure options --download-cgns --download-hdf5 --download-openmpi --download-triangle --with-fc=0 PETSC_ARCH=arch-darwin-c-debug --download-cgns [0]PETSC ERROR: #1 DMPlexCGNSGetPermutation_Internal() at /Users/Junming/Packages/petsc/src/dm/impls/plex/cgns/plexcgns2.c:533 [0]PETSC ERROR: #2 DMView_PlexCGNS() at /Users/Junming/Packages/petsc/src/dm/impls/plex/cgns/plexcgns2.c:769 [0]PETSC ERROR: #3 DMView_Plex() at /Users/Junming/Packages/petsc/src/dm/impls/plex/plex.c:1801 [0]PETSC ERROR: #4 DMView() at /Users/Junming/Packages/petsc/src/dm/interface/dm.c:996 [0]PETSC ERROR: #5 PetscObjectView() at /Users/Junming/Packages/petsc/src/sys/objects/destroy.c:78 [0]PETSC ERROR: #6 PetscObjectViewFromOptions() at /Users/Junming/Packages/petsc/src/sys/objects/destroy.c:128 [0]PETSC ERROR: #7 DMViewFromOptions() at /Users/Junming/Packages/petsc/src/dm/interface/dm.c:940 [0]PETSC ERROR: #8 CreateDiscretization() at ex33.c:232 [0]PETSC ERROR: #9 main() at ex33.c:263 [0]PETSC ERROR: PETSc Option Table entries: [0]PETSC ERROR: -dm_coord_petscspace_degree 5 (source: command line) [0]PETSC ERROR: -dm_coord_space 0 (source: command line) [0]PETSC ERROR: -dm_plex_box_faces 1,1 (source: command line) [0]PETSC ERROR: -dm_plex_simplex 0 (source: command line) [0]PETSC ERROR: -dm_refine 1 (source: command line) [0]PETSC ERROR: -dm_view cgns:test.cgns (source: command line) [0]PETSC ERROR: -mesh_transform annulus (source: command line) [0]PETSC ERROR: -petscspace_degree 5 (source: command line) [0]PETSC ERROR: ----------------End of Error Message -------send entire error message to [email protected] -------------------------------------------------------------------------- MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_SELF with errorcode 56.
Does cgns work for degree >= 4?
Junming
________________________________ From: Jed Brown <[email protected]> Sent: Monday, June 12, 2023 19:07 To: Duan Junming; Matthew Knepley Cc: [email protected] Subject: Re: [petsc-users] dm_view of high-order geometry/solution
And here's an MR to do what you want without any code/arg changes.
https://gitlab.com/petsc/petsc/-/merge_requests/6588
Jed Brown <[email protected]> writes:
Duan Junming <[email protected]> writes:
Dear Jed,
Thank you for the suggestion.
When I run tests/ex33.c with
./ex33 -dm_plex_simplex 0 -dm_plex_box_faces 1,1 -mesh_transform annulus -dm_coord_space 0 -dm_coord_petscspace_degree 3 -dm_refine 1 -dm_view cgns:test.cgns
and load it using Paraview,
the mesh is still with straight lines.
Ah, the viewer is keyed on the field (since the CGNS as supported by Paraview specifies coordinates and fields in the same space). That doesn't exist in your case. If you apply this patch and add `-petscspace_degre 3` to your command, you'll see that high order information is present. Paraview doesn't render as curves in all views, but it has the data.
diff --git i/src/dm/impls/plex/tests/ex33.c w/src/dm/impls/plex/tests/ex33.c index 803095bc082..590facfa4f4 100644 --- i/src/dm/impls/plex/tests/ex33.c +++ w/src/dm/impls/plex/tests/ex33.c @@ -198,7 +198,6 @@ PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *ctx, DM *dm) default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Unknown mesh transform %d", ctx->meshTransform); } - PetscCall(DMViewFromOptions(*dm, NULL, "-dm_view")); PetscFunctionReturn(PETSC_SUCCESS); }
@@ -227,6 +226,7 @@ static PetscErrorCode CreateDiscretization(DM dm, AppCtx *ctx) PetscCall(DMCreateDS(dm)); PetscCall(DMGetDS(dm, &ds)); PetscCall(PetscDSSetObjective(ds, 0, volume)); + PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); PetscFunctionReturn(PETSC_SUCCESS); }
I can update the viewer to handle the degenerate case of no field (all my models have fields).
participants (3)
-
Duan Junming -
Jed Brown -
Matthew Knepley