Actual source code: ex100.c
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
5: #include <petscdmplex.h>
6: #include <petscviewerhdf5.h>
7: #define EX "ex100.c"
9: static PetscErrorCode project_function(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx)
10: {
11: PetscReal x_tot = 0;
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: }
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;
27: PetscFunctionBeginUser;
28: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
29: comm = PETSC_COMM_WORLD;
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();
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"));
43: { // Setup fe to load in the initial condition data
44: PetscFE fe;
45: PetscInt dim;
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: }
54: PetscErrorCode (*funcs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *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));
63: PetscCall(DMDestroy(&dm_create));
64: }
66: {
67: PetscSection coord_section;
68: PetscInt cStart, cEnd;
69: Vec coords;
70: PetscInt cdim;
71: PetscReal ref_cell_bounding_box_size[3];
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};
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: }
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));
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: }
108: PetscCall(PetscFinalize());
109: return 0;
110: }
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: TEST*/