1 static char help[] = "Tests for submesh with both CG and DG coordinates\n\n"; 2 3 #include <petscdmplex.h> 4 #include <petsc/private/dmimpl.h> 5 6 int main(int argc, char **argv) 7 { 8 PetscInt dim = 1, d, cStart, cEnd, c, q, degree = 2, coordSize, offset; 9 PetscReal R = 1.0; 10 DM dm, coordDM, subdm; 11 PetscSection coordSec; 12 Vec coordVec; 13 PetscScalar *coords; 14 DMLabel filter; 15 const PetscInt filterValue = 1; 16 MPI_Comm comm; 17 PetscMPIInt size; 18 19 PetscFunctionBeginUser; 20 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 21 comm = PETSC_COMM_WORLD; 22 PetscCallMPI(MPI_Comm_size(comm, &size)); 23 if (size != 1) { 24 PetscCall(PetscPrintf(comm, "This example is specifically designed for comm size == 1.\n")); 25 PetscCall(PetscFinalize()); 26 return 0; 27 } 28 /* Make a unit circle with 16 elements. */ 29 PetscCall(DMPlexCreateSphereMesh(comm, dim, PETSC_TRUE, R, &dm)); 30 PetscUseTypeMethod(dm, createcellcoordinatedm, &coordDM); 31 PetscCall(DMSetCellCoordinateDM(dm, coordDM)); 32 PetscCall(DMDestroy(&coordDM)); 33 PetscCall(DMGetCellCoordinateSection(dm, &coordSec)); 34 PetscCall(PetscSectionSetNumFields(coordSec, 1)); 35 PetscCall(PetscSectionSetFieldComponents(coordSec, 0, dim)); 36 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 37 PetscCall(PetscSectionSetChart(coordSec, cStart, cEnd)); 38 for (c = cStart; c < cEnd; ++c) { 39 PetscCall(PetscSectionSetDof(coordSec, c, dim * (degree + 1))); 40 PetscCall(PetscSectionSetFieldDof(coordSec, c, 0, dim * (degree + 1))); 41 } 42 PetscCall(PetscSectionSetUp(coordSec)); 43 PetscCall(PetscSectionGetStorageSize(coordSec, &coordSize)); 44 PetscCall(VecCreate(PETSC_COMM_SELF, &coordVec)); 45 PetscCall(PetscObjectSetName((PetscObject)coordVec, "cellcoordinates")); 46 PetscCall(VecSetBlockSize(coordVec, dim)); 47 PetscCall(VecSetSizes(coordVec, coordSize, PETSC_DETERMINE)); 48 PetscCall(VecSetType(coordVec, VECSTANDARD)); 49 PetscCall(VecGetArray(coordVec, &coords)); 50 for (c = cStart; c < cEnd; ++c) { 51 PetscCall(PetscSectionGetOffset(coordSec, c, &offset)); 52 for (q = 0; q < degree + 1; ++q) { 53 /* Make some DG coordinates. Note that dim = 1.*/ 54 for (d = 0; d < dim; ++d) coords[offset + dim * q + d] = 100. + (PetscScalar)c + (1.0 / (PetscScalar)degree) * (PetscScalar)q; 55 } 56 } 57 PetscCall(VecRestoreArray(coordVec, &coords)); 58 PetscCall(DMSetCellCoordinatesLocal(dm, coordVec)); 59 PetscCall(VecDestroy(&coordVec)); 60 /* Make submesh. */ 61 PetscCall(DMLabelCreate(PETSC_COMM_SELF, "filter", &filter)); 62 PetscCall(DMLabelSetValue(filter, 15, filterValue)); /* last cell */ 63 PetscCall(DMPlexFilter(dm, filter, filterValue, PETSC_FALSE, PETSC_FALSE, NULL, &subdm)); 64 PetscCall(DMLabelDestroy(&filter)); 65 PetscCall(PetscObjectSetName((PetscObject)subdm, "Example_SubDM")); 66 PetscCall(DMViewFromOptions(subdm, NULL, "-dm_view")); 67 PetscCall(DMDestroy(&subdm)); 68 PetscCall(DMDestroy(&dm)); 69 PetscCall(PetscFinalize()); 70 return 0; 71 } 72 73 /*TEST 74 75 test: 76 suffix: 0 77 args: -dm_view ascii::ascii_info_detail 78 79 TEST*/ 80