1 static char help[] = "Test FEM layout with constraints\n\n";
2
3 #include <petsc.h>
4
main(int argc,char ** argv)5 int main(int argc, char **argv)
6 {
7 DM dm, pdm;
8 PetscSection section;
9 const PetscInt field = 0;
10 char ifilename[PETSC_MAX_PATH_LEN];
11 PetscInt sdim, s, pStart, pEnd, p, numVS, numPoints;
12 PetscInt constraints[1];
13 IS setIS, pointIS;
14 const PetscInt *setID, *pointID;
15
16 PetscFunctionBeginUser;
17 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
18 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FEM Layout Options", "ex98");
19 PetscCall(PetscOptionsString("-i", "Filename to read", "ex98", ifilename, ifilename, sizeof(ifilename), NULL));
20 PetscOptionsEnd();
21
22 PetscCall(DMPlexCreateFromFile(PETSC_COMM_WORLD, ifilename, NULL, PETSC_TRUE, &dm));
23 PetscCall(DMPlexDistributeSetDefault(dm, PETSC_FALSE));
24 PetscCall(DMSetFromOptions(dm));
25
26 PetscCall(DMPlexDistribute(dm, 0, NULL, &pdm));
27 if (pdm) {
28 PetscCall(DMDestroy(&dm));
29 dm = pdm;
30 }
31 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
32
33 /* create a section */
34 PetscCall(DMGetDimension(dm, &sdim));
35 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion));
36 PetscCall(PetscSectionSetNumFields(section, 1));
37 PetscCall(PetscSectionSetFieldName(section, field, "U"));
38 PetscCall(PetscSectionSetFieldComponents(section, field, sdim));
39 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
40 PetscCall(PetscSectionSetChart(section, pStart, pEnd));
41
42 /* initialize the section storage for a P1 field */
43 PetscCall(DMPlexGetDepthStratum(dm, 0, &pStart, &pEnd));
44 for (p = pStart; p < pEnd; ++p) {
45 PetscCall(PetscSectionSetDof(section, p, sdim));
46 PetscCall(PetscSectionSetFieldDof(section, p, 0, sdim));
47 }
48
49 /* add constraints at all vertices belonging to a vertex set */
50 /* first pass is to reserve space */
51 PetscCall(DMGetLabelSize(dm, "Vertex Sets", &numVS));
52 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "# Vertex set: %" PetscInt_FMT "\n", numVS));
53 PetscCall(DMGetLabelIdIS(dm, "Vertex Sets", &setIS));
54 PetscCall(ISGetIndices(setIS, &setID));
55 for (s = 0; s < numVS; ++s) {
56 PetscCall(DMGetStratumIS(dm, "Vertex Sets", setID[s], &pointIS));
57 PetscCall(DMGetStratumSize(dm, "Vertex Sets", setID[s], &numPoints));
58 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "set %" PetscInt_FMT " size: %" PetscInt_FMT "\n", s, numPoints));
59 PetscCall(ISGetIndices(pointIS, &pointID));
60 for (p = 0; p < numPoints; ++p) {
61 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\t point %" PetscInt_FMT "\n", pointID[p]));
62 PetscCall(PetscSectionSetConstraintDof(section, pointID[p], 1));
63 PetscCall(PetscSectionSetFieldConstraintDof(section, pointID[p], field, 1));
64 }
65 PetscCall(ISRestoreIndices(pointIS, &pointID));
66 PetscCall(ISDestroy(&pointIS));
67 }
68
69 PetscCall(PetscSectionSetUp(section));
70
71 /* add constraints at all vertices belonging to a vertex set */
72 /* second pass is to assign constraints to a specific component / dof */
73 for (s = 0; s < numVS; ++s) {
74 PetscCall(DMGetStratumIS(dm, "Vertex Sets", setID[s], &pointIS));
75 PetscCall(DMGetStratumSize(dm, "Vertex Sets", setID[s], &numPoints));
76 PetscCall(ISGetIndices(pointIS, &pointID));
77 for (p = 0; p < numPoints; ++p) {
78 constraints[0] = setID[s] % sdim;
79 PetscCall(PetscSectionSetConstraintIndices(section, pointID[p], constraints));
80 PetscCall(PetscSectionSetFieldConstraintIndices(section, pointID[p], field, constraints));
81 }
82 PetscCall(ISRestoreIndices(pointIS, &pointID));
83 PetscCall(ISDestroy(&pointIS));
84 }
85 PetscCall(ISRestoreIndices(setIS, &setID));
86 PetscCall(ISDestroy(&setIS));
87 PetscCall(PetscObjectViewFromOptions((PetscObject)section, NULL, "-dm_section_view"));
88
89 PetscCall(PetscSectionDestroy(§ion));
90 PetscCall(DMDestroy(&dm));
91
92 PetscCall(PetscFinalize());
93 return 0;
94 }
95
96 /*TEST
97 build:
98 requires: exodusii pnetcdf !complex
99 testset:
100 args: -i ${wPETSC_DIR}/share/petsc/datafiles/meshes/SquareFaceSet.exo -dm_view -dm_section_view
101 nsize: 1
102
103 test:
104 suffix: 0
105 args:
106
107 TEST*/
108