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