1 static char help[] = "Verify isoperiodic cone corrections"; 2 3 #include <petscdmplex.h> 4 5 // Creates periodic solution on a [0,1] x D domain for D dimension 6 static PetscErrorCode project_function(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) 7 { 8 PetscReal x_tot = 0; 9 10 PetscFunctionBeginUser; 11 for (PetscInt d = 0; d < dim; d++) x_tot += x[d]; 12 for (PetscInt c = 0; c < Nc; c++) { 13 PetscScalar value = PetscSinReal(2 * M_PI * x_tot); 14 if (PetscAbsScalar(value) < 1e-7) value = 0.; 15 u[c] = value; 16 } 17 PetscFunctionReturn(PETSC_SUCCESS); 18 } 19 20 PetscErrorCode CreateFEField(DM dm) 21 { 22 PetscInt degree; 23 24 PetscFunctionBeginUser; 25 { // Get degree of the coords section 26 PetscFE fe_coords; 27 PetscSpace coord_space; 28 DM cdm; 29 30 PetscCall(DMGetCoordinateDM(dm, &cdm)); 31 PetscCall(DMGetField(cdm, 0, NULL, (PetscObject *)&fe_coords)); 32 PetscCall(PetscFEGetBasisSpace(fe_coords, &coord_space)); 33 PetscCall(PetscSpaceGetDegree(coord_space, °ree, NULL)); 34 } 35 36 PetscCall(DMClearFields(dm)); 37 PetscCall(DMSetLocalSection(dm, NULL)); // See https://gitlab.com/petsc/petsc/-/issues/1669 38 39 { // Setup fe to load in the initial condition data 40 PetscFE fe; 41 PetscInt dim; 42 43 PetscCall(DMGetDimension(dm, &dim)); 44 PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, degree, PETSC_DETERMINE, &fe)); 45 PetscCall(DMAddField(dm, NULL, (PetscObject)fe)); 46 PetscCall(DMCreateDS(dm)); 47 PetscCall(PetscFEDestroy(&fe)); 48 } 49 PetscFunctionReturn(PETSC_SUCCESS); 50 } 51 52 int main(int argc, char **argv) 53 { 54 MPI_Comm comm; 55 DM dm = NULL; 56 Vec V, V_G2L, V_local; 57 PetscReal norm; 58 PetscErrorCode (*funcs)(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar *u, void *ctx) = {project_function}; 59 60 PetscFunctionBeginUser; 61 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 62 comm = PETSC_COMM_WORLD; 63 64 PetscCall(DMCreate(comm, &dm)); 65 PetscCall(DMSetType(dm, DMPLEX)); 66 PetscCall(PetscObjectSetName((PetscObject)dm, "ex101_dm")); 67 PetscCall(DMSetFromOptions(dm)); 68 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 69 PetscCall(CreateFEField(dm)); 70 71 // Verify that projected function on global vector (then projected onto local vector) is equal to projected function onto a local vector 72 PetscCall(DMGetLocalVector(dm, &V_G2L)); 73 PetscCall(DMGetGlobalVector(dm, &V)); 74 PetscCall(DMProjectFunction(dm, 0, &funcs, NULL, INSERT_VALUES, V)); 75 PetscCall(DMGlobalToLocal(dm, V, INSERT_VALUES, V_G2L)); 76 77 PetscCall(DMGetLocalVector(dm, &V_local)); 78 PetscCall(DMProjectFunctionLocal(dm, 0, &funcs, NULL, INSERT_VALUES, V_local)); 79 80 PetscCall(VecAXPY(V_local, -1, V_G2L)); 81 PetscCall(VecNorm(V_local, NORM_2, &norm)); 82 PetscReal tol = PetscDefined(USE_REAL___FLOAT128) ? 1e-12 : 1e4 * PETSC_MACHINE_EPSILON; 83 PetscCheck(norm < tol, PETSC_COMM_SELF, PETSC_ERR_PLIB, "GlobalToLocal result does not match Local projection by norm %g", (double)norm); 84 85 PetscCall(DMRestoreGlobalVector(dm, &V)); 86 PetscCall(DMRestoreLocalVector(dm, &V_G2L)); 87 PetscCall(DMRestoreLocalVector(dm, &V_local)); 88 PetscCall(DMDestroy(&dm)); 89 90 PetscCall(PetscFinalize()); 91 return 0; 92 } 93 94 /*TEST 95 96 test: 97 nsize: 2 98 args: -dm_plex_shape zbox -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 1,2,3 -dm_coord_space -dm_coord_petscspace_degree 3 99 args: -dm_plex_box_bd periodic,periodic,periodic -dm_view ::ascii_info_detail -petscpartitioner_type simple 100 101 test: 102 requires: cgns 103 suffix: cgns 104 nsize: 3 105 args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/2x2x2_Q3_wave.cgns -dm_plex_cgns_parallel -dm_view ::ascii_info_detail -dm_plex_box_label true -dm_plex_box_label_bd periodic,periodic,periodic -petscpartitioner_type simple 106 107 TEST*/ 108