1 static char help[] = "Example program demonstrating projection between particle and finite element spaces\n\n"; 2 3 #include <petscdmplex.h> 4 #include <petscds.h> 5 #include <petscdmswarm.h> 6 #include <petscksp.h> 7 8 int main(int argc, char **argv) 9 { 10 DM dm, sw; 11 PetscFE fe; 12 Vec u_f; 13 DMPolytopeType ct; 14 PetscInt dim, Nc = 1, faces[3]; 15 PetscInt Np = 10, field = 0, zero = 0, bs, cStart; 16 PetscReal energy_0 = 0, energy_1 = 0; 17 PetscReal lo[3], hi[3], h[3]; 18 PetscBool removePoints = PETSC_TRUE; 19 PetscReal *wq, *coords; 20 PetscDataType dtype; 21 22 PetscFunctionBeginUser; 23 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 24 /* Create a mesh */ 25 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm)); 26 PetscCall(DMSetType(dm, DMPLEX)); 27 PetscCall(DMSetFromOptions(dm)); 28 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view")); 29 30 PetscCall(DMGetDimension(dm, &dim)); 31 bs = dim; 32 PetscCall(PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", faces, &bs, NULL)); 33 PetscCall(PetscOptionsGetInt(NULL, NULL, "-np", &Np, NULL)); 34 PetscCall(DMGetBoundingBox(dm, lo, hi)); 35 for (PetscInt i = 0; i < dim; ++i) { 36 h[i] = (hi[i] - lo[i]) / faces[i]; 37 PetscCall(PetscPrintf(PETSC_COMM_SELF, " lo = %g hi = %g n = %" PetscInt_FMT " h = %g\n", (double)lo[i], (double)hi[i], faces[i], (double)h[i])); 38 } 39 // Create FE space 40 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, NULL)); 41 PetscCall(DMPlexGetCellType(dm, cStart, &ct)); 42 PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, Nc, ct, NULL, PETSC_DECIDE, &fe)); 43 PetscCall(PetscFESetFromOptions(fe)); 44 PetscCall(PetscObjectSetName((PetscObject)fe, "fe")); 45 PetscCall(DMSetField(dm, field, NULL, (PetscObject)fe)); 46 PetscCall(DMCreateDS(dm)); 47 PetscCall(PetscFEDestroy(&fe)); 48 PetscCall(DMCreateGlobalVector(dm, &u_f)); 49 // Create particle swarm 50 PetscCall(DMCreate(PETSC_COMM_SELF, &sw)); 51 PetscCall(DMSetType(sw, DMSWARM)); 52 PetscCall(DMSetDimension(sw, dim)); 53 PetscCall(DMSwarmSetType(sw, DMSWARM_PIC)); 54 PetscCall(DMSwarmSetCellDM(sw, dm)); 55 PetscCall(DMSwarmRegisterPetscDatatypeField(sw, "w_q", Nc, PETSC_SCALAR)); 56 PetscCall(DMSwarmFinalizeFieldRegister(sw)); 57 PetscCall(DMSwarmSetLocalSizes(sw, Np, zero)); 58 PetscCall(DMSetFromOptions(sw)); 59 PetscCall(DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&wq)); 60 PetscCall(DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords)); 61 for (PetscInt p = 0; p < Np; ++p) { 62 coords[p * 2 + 0] = -PetscCosReal((PetscReal)(p + 1) / (PetscReal)(Np + 1) * PETSC_PI); 63 coords[p * 2 + 1] = PetscSinReal((PetscReal)(p + 1) / (PetscReal)(Np + 1) * PETSC_PI); 64 wq[p] = 1.0; 65 energy_0 += wq[p] * (PetscSqr(coords[p * 2 + 0]) + PetscSqr(coords[p * 2 + 1])); 66 } 67 PetscCall(DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords)); 68 PetscCall(DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void **)&wq)); 69 PetscCall(DMSwarmMigrate(sw, removePoints)); 70 PetscCall(PetscObjectSetName((PetscObject)sw, "Particle Grid")); 71 PetscCall(DMViewFromOptions(sw, NULL, "-swarm_view")); 72 73 // Project between particles and continuum field 74 const char *fieldnames[1] = {"w_q"}; 75 Vec fields[1] = {u_f}; 76 PetscCall(DMSwarmProjectFields(sw, NULL, 1, fieldnames, fields, SCATTER_FORWARD)); 77 PetscCall(DMSwarmProjectFields(sw, NULL, 1, fieldnames, fields, SCATTER_REVERSE)); 78 79 // Compute energy 80 PetscCall(DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&wq)); 81 PetscCall(DMSwarmGetField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords)); 82 for (PetscInt p = 0; p < Np; ++p) energy_1 += wq[p] * (PetscSqr(coords[p * 2 + 0]) + PetscSqr(coords[p * 2 + 1])); 83 PetscCall(DMSwarmRestoreField(sw, "DMSwarmPIC_coor", &bs, &dtype, (void **)&coords)); 84 PetscCall(DMSwarmRestoreField(sw, "w_q", &bs, &dtype, (void **)&wq)); 85 PetscCall(PetscPrintf(PETSC_COMM_SELF, "Energy = %20.12e error = %20.12e\n", (double)energy_0, (double)((energy_1 - energy_0) / energy_0))); 86 87 // Cleanup 88 PetscCall(VecDestroy(&u_f)); 89 PetscCall(DMDestroy(&sw)); 90 PetscCall(DMDestroy(&dm)); 91 PetscCall(PetscFinalize()); 92 return 0; 93 } 94 95 /*TEST 96 97 build: 98 requires: !complex 99 100 test: 101 suffix: 0 102 requires: double triangle 103 args: -dm_plex_simplex 0 -dm_plex_box_faces 4,2 -dm_plex_box_lower -2.0,0.0 -dm_plex_box_upper 2.0,2.0 \ 104 -np 50 -petscspace_degree 2 \ 105 -ptof_ksp_type cg -ptof_pc_type ilu -ptof_ksp_rtol 1.e-14 \ 106 -ftop_ksp_type lsqr -ftop_pc_type none -ftop_ksp_rtol 1.e-14 \ 107 -dm_view -swarm_view 108 filter: grep -v DM_ | grep -v atomic 109 110 test: 111 suffix: bjacobi 112 requires: double triangle 113 args: -dm_plex_simplex 0 -dm_plex_box_faces 4,2 -dm_plex_box_lower -2.0,0.0 -dm_plex_box_upper 2.0,2.0 \ 114 -np 50 -petscspace_degree 2 -dm_plex_hash_location \ 115 -ptof_ksp_type cg -ptof_pc_type ilu -ptof_ksp_rtol 1.e-14 \ 116 -ftop_ksp_type lsqr -ftop_pc_type bjacobi -ftop_sub_pc_type lu -ftop_sub_pc_factor_shift_type nonzero \ 117 -dm_view -swarm_view -ftop_ksp_rtol 1.e-14 118 filter: grep -v DM_ | grep -v atomic 119 120 TEST*/ 121