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 9 static PetscErrorCode project_function(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *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 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, 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)); 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