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
main(int argc,char ** argv)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, PetscObjectComm((PetscObject)dm), 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