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