1 #include "petscsf.h"
2 static char help[] = "Test CGNS writing output with isoperiodic boundaries\n\n";
3 // Also tests DMSetCoordinateDisc() for isoperiodic boundaries and projection = true
4
5 #include <petscdmplex.h>
6 #include <petscviewerhdf5.h>
7 #define EX "ex100.c"
8
project_function(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar * u,PetscCtx ctx)9 static PetscErrorCode project_function(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx)
10 {
11 PetscReal x_tot = 0;
12
13 PetscFunctionBeginUser;
14 for (PetscInt d = 0; d < dim; d++) x_tot += x[d];
15 for (PetscInt c = 0; c < Nc; c++) u[c] = sin(2 * M_PI * x_tot);
16 PetscFunctionReturn(PETSC_SUCCESS);
17 }
18
main(int argc,char ** argv)19 int main(int argc, char **argv)
20 {
21 DM dm_create, dm_read;
22 Vec V;
23 char file[PETSC_MAX_PATH_LEN];
24 MPI_Comm comm;
25 PetscInt solution_degree = 2;
26
27 PetscFunctionBeginUser;
28 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
29 comm = PETSC_COMM_WORLD;
30
31 PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");
32 PetscCall(PetscOptionsInt("-solution_degree", "The input CGNS file", EX, solution_degree, &solution_degree, NULL));
33 PetscCall(PetscOptionsString("-file", "The input CGNS file", EX, file, file, sizeof(file), NULL));
34 PetscOptionsEnd();
35
36 { // Create DM
37 PetscCall(DMCreate(comm, &dm_create));
38 PetscCall(DMSetType(dm_create, DMPLEX));
39 PetscCall(DMSetOptionsPrefix(dm_create, "create_"));
40 PetscCall(DMSetFromOptions(dm_create));
41 PetscCall(DMViewFromOptions(dm_create, NULL, "-dm_create_view"));
42
43 { // Setup fe to load in the initial condition data
44 PetscFE fe;
45 PetscInt dim;
46
47 PetscCall(DMGetDimension(dm_create, &dim));
48 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, solution_degree, PETSC_DETERMINE, &fe));
49 PetscCall(DMAddField(dm_create, NULL, (PetscObject)fe));
50 PetscCall(DMCreateDS(dm_create));
51 PetscCall(PetscFEDestroy(&fe));
52 }
53
54 PetscErrorCode (*funcs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, PetscCtx ctx) = {project_function};
55 PetscCall(DMGetGlobalVector(dm_create, &V));
56 PetscCall(DMProjectFunction(dm_create, 0, &funcs, NULL, INSERT_VALUES, V));
57 PetscViewer viewer;
58 PetscCall(PetscViewerCGNSOpen(comm, file, FILE_MODE_WRITE, &viewer));
59 PetscCall(VecView(V, viewer));
60 PetscCall(PetscViewerDestroy(&viewer));
61 PetscCall(DMRestoreGlobalVector(dm_create, &V));
62
63 PetscCall(DMDestroy(&dm_create));
64 }
65
66 {
67 PetscSection coord_section;
68 PetscInt cStart, cEnd;
69 Vec coords;
70 PetscInt cdim;
71 PetscReal ref_cell_bounding_box_size[3];
72
73 PetscCall(DMPlexCreateFromFile(comm, file, "ex100_plex", PETSC_TRUE, &dm_read));
74 PetscCall(DMViewFromOptions(dm_read, NULL, "-dm_read_view"));
75 PetscCall(DMGetCoordinateDim(dm_read, &cdim));
76 PetscCall(DMGetCoordinateSection(dm_read, &coord_section));
77 PetscCall(DMGetCoordinates(dm_read, &coords));
78 PetscCall(DMPlexGetHeightStratum(dm_read, 0, &cStart, &cEnd));
79 for (PetscInt cell = cStart; cell < cEnd; cell++) {
80 PetscInt num_closure = 0;
81 PetscScalar *cell_coords = NULL;
82 PetscReal cell_bounding_box_size[3], difference;
83 PetscReal min[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
84 PetscReal max[3] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
85
86 PetscCall(DMPlexVecGetClosure(dm_read, coord_section, coords, cell, &num_closure, &cell_coords));
87 for (PetscInt n = 0; n < num_closure / cdim; n++) {
88 for (PetscInt d = 0; d < cdim; ++d) {
89 min[d] = PetscMin(min[d], PetscRealPart(cell_coords[n * cdim + d]));
90 max[d] = PetscMax(max[d], PetscRealPart(cell_coords[n * cdim + d]));
91 }
92 }
93
94 for (PetscInt d = 0; d < cdim; d++) cell_bounding_box_size[d] = max[d] - min[d];
95 if (cell == cStart) PetscCall(PetscArraycpy(ref_cell_bounding_box_size, cell_bounding_box_size, cdim));
96
97 for (PetscInt d = 0; d < cdim; d++) {
98 difference = PetscAbsReal((ref_cell_bounding_box_size[d] - cell_bounding_box_size[d]) / ref_cell_bounding_box_size[d]);
99 if (difference > PETSC_MACHINE_EPSILON * 100) {
100 PetscPrintf(comm, "Cell %" PetscInt_FMT " doesn't match bounding box size of Cell %" PetscInt_FMT " in dimension %" PetscInt_FMT ". Relative difference = %g\n", cell, cStart, d, (double)difference);
101 }
102 }
103 PetscCall(DMPlexVecRestoreClosure(dm_read, coord_section, coords, cell, &num_closure, &cell_coords));
104 }
105 PetscCall(DMDestroy(&dm_read));
106 }
107
108 PetscCall(PetscFinalize());
109 return 0;
110 }
111
112 /*TEST
113 build:
114 requires: cgns
115 test:
116 # Coordinates of dm_create are linear, but the CGNS writer projects them to quadratic to match the solution.
117 # The verification checks that all the cells of the resulting file are the same size
118 args: -create_dm_plex_shape zbox -create_dm_plex_box_faces 2,2 -create_dm_plex_box_bd periodic,periodic -file test.cgns -dm_plex_cgns_parallel
119 output_file: output/empty.out
120 TEST*/
121