Actual source code: ex98.c
1: static char help[] = "Test FEM layout with constraints\n\n";
3: #include <petsc.h>
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;
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();
22: PetscCall(DMPlexCreateFromFile(PETSC_COMM_WORLD, ifilename, NULL, PETSC_TRUE, &dm));
23: PetscCall(DMPlexDistributeSetDefault(dm, PETSC_FALSE));
24: PetscCall(DMSetFromOptions(dm));
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"));
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));
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: }
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: }
69: PetscCall(PetscSectionSetUp(section));
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"));
89: PetscCall(PetscSectionDestroy(§ion));
90: PetscCall(DMDestroy(&dm));
92: PetscCall(PetscFinalize());
93: return 0;
94: }
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
103: test:
104: suffix: 0
105: args:
107: TEST*/