1 static char help[] = "Create a Plex sphere from quads and create a P1 section\n\n"; 2 3 #include <petscdmplex.h> 4 5 typedef struct { 6 PetscInt dim; /* Topological problem dimension */ 7 PetscBool simplex; /* Mesh with simplices */ 8 } AppCtx; 9 10 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 11 { 12 PetscErrorCode ierr; 13 14 PetscFunctionBeginUser; 15 options->dim = 2; 16 options->simplex = PETSC_FALSE; 17 18 ierr = PetscOptionsBegin(comm, "", "Sphere Mesh Options", "DMPLEX");CHKERRQ(ierr); 19 ierr = PetscOptionsRangeInt("-dim", "Problem dimension", "ex7.c", options->dim, &options->dim, NULL,1,3);CHKERRQ(ierr); 20 ierr = PetscOptionsBool("-simplex", "Use simplices, or tensor product cells", "ex7.c", options->simplex, &options->simplex, NULL);CHKERRQ(ierr); 21 ierr = PetscOptionsEnd(); 22 PetscFunctionReturn(0); 23 } 24 25 static PetscErrorCode ProjectToUnitSphere(DM dm) 26 { 27 Vec coordinates; 28 PetscScalar *coords; 29 PetscInt Nv, v, bs, dim, d; 30 PetscErrorCode ierr; 31 32 PetscFunctionBeginUser; 33 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 34 ierr = VecGetLocalSize(coordinates, &Nv);CHKERRQ(ierr); 35 ierr = VecGetBlockSize(coordinates, &bs);CHKERRQ(ierr); 36 ierr = DMGetCoordinateDim(dm, &dim);CHKERRQ(ierr); 37 if (dim != bs) SETERRQ2(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Coordinate bs %D does not match dim %D",bs,dim); 38 Nv /= dim; 39 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 40 for (v = 0; v < Nv; ++v) { 41 PetscReal r = 0.0; 42 43 for (d = 0; d < dim; ++d) r += PetscSqr(PetscRealPart(coords[v*dim+d])); 44 r = PetscSqrtReal(r); 45 for (d = 0; d < dim; ++d) coords[v*dim+d] /= r; 46 } 47 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 48 PetscFunctionReturn(0); 49 } 50 51 static PetscErrorCode SetupSection(DM dm) 52 { 53 PetscSection s; 54 PetscInt vStart, vEnd, v; 55 PetscErrorCode ierr; 56 57 PetscFunctionBeginUser; 58 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 59 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), &s);CHKERRQ(ierr); 60 ierr = PetscSectionSetNumFields(s, 1);CHKERRQ(ierr); 61 ierr = PetscSectionSetFieldComponents(s, 0, 1);CHKERRQ(ierr); 62 ierr = PetscSectionSetChart(s, vStart, vEnd);CHKERRQ(ierr); 63 for (v = vStart; v < vEnd; ++v) { 64 ierr = PetscSectionSetDof(s, v, 1);CHKERRQ(ierr); 65 ierr = PetscSectionSetFieldDof(s, v, 0, 1);CHKERRQ(ierr); 66 } 67 ierr = PetscSectionSetUp(s);CHKERRQ(ierr); 68 ierr = DMSetLocalSection(dm, s);CHKERRQ(ierr); 69 ierr = PetscSectionDestroy(&s);CHKERRQ(ierr); 70 PetscFunctionReturn(0); 71 } 72 73 int main(int argc, char **argv) 74 { 75 DM dm; 76 Vec u; 77 AppCtx ctx; 78 PetscErrorCode ierr; 79 80 ierr = PetscInitialize(&argc, &argv, NULL,help);if (ierr) return ierr; 81 ierr = ProcessOptions(PETSC_COMM_WORLD, &ctx);CHKERRQ(ierr); 82 ierr = DMPlexCreateSphereMesh(PETSC_COMM_WORLD, ctx.dim, ctx.simplex, &dm);CHKERRQ(ierr); 83 ierr = PetscObjectSetName((PetscObject) dm, "Sphere");CHKERRQ(ierr); 84 /* Distribute mesh over processes */ 85 { 86 DM dmDist = NULL; 87 PetscPartitioner part; 88 ierr = DMPlexGetPartitioner(dm, &part);CHKERRQ(ierr); 89 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 90 ierr = DMPlexDistribute(dm, 0, NULL, &dmDist);CHKERRQ(ierr); 91 if (dmDist) { 92 ierr = DMDestroy(&dm);CHKERRQ(ierr); 93 dm = dmDist; 94 } 95 } 96 ierr = DMSetFromOptions(dm);CHKERRQ(ierr); 97 ierr = ProjectToUnitSphere(dm);CHKERRQ(ierr); 98 ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr); 99 ierr = SetupSection(dm);CHKERRQ(ierr); 100 ierr = DMGetGlobalVector(dm, &u);CHKERRQ(ierr); 101 ierr = VecSet(u, 2);CHKERRQ(ierr); 102 ierr = VecViewFromOptions(u, NULL, "-vec_view");CHKERRQ(ierr); 103 ierr = DMRestoreGlobalVector(dm, &u);CHKERRQ(ierr); 104 ierr = DMDestroy(&dm);CHKERRQ(ierr); 105 ierr = PetscFinalize(); 106 return ierr; 107 } 108 109 /*TEST 110 111 test: 112 suffix: 2d_quad 113 requires: !__float128 114 args: -dm_view 115 116 test: 117 suffix: 2d_quad_parallel 118 requires: !__float128 119 args: -dm_view -petscpartitioner_type simple 120 nsize: 2 121 122 test: 123 suffix: 2d_tri 124 requires: !__float128 125 args: -simplex -dm_view 126 127 test: 128 suffix: 2d_tri_parallel 129 requires: !__float128 130 args: -simplex -dm_view -petscpartitioner_type simple 131 nsize: 2 132 133 test: 134 suffix: 3d_tri 135 requires: !__float128 136 args: -dim 3 -simplex -dm_view 137 138 test: 139 suffix: 3d_tri_parallel 140 requires: !__float128 141 args: -dim 3 -simplex -dm_view -petscpartitioner_type simple 142 nsize: 2 143 144 TEST*/ 145