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;
17: PetscInitialize(&argc, &argv, NULL, help);
18: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FEM Layout Options", "ex98");
19: PetscOptionsString("-i", "Filename to read", "ex98", ifilename, ifilename, sizeof(ifilename), NULL);
20: PetscOptionsEnd();
22: DMPlexCreateFromFile(PETSC_COMM_WORLD, ifilename, NULL, PETSC_TRUE, &dm);
23: DMPlexDistributeSetDefault(dm, PETSC_FALSE);
24: DMSetFromOptions(dm);
26: DMPlexDistribute(dm, 0, NULL, &pdm);
27: if (pdm) {
28: DMDestroy(&dm);
29: dm = pdm;
30: }
31: DMViewFromOptions(dm, NULL, "-dm_view");
33: /* create a section */
34: DMGetDimension(dm, &sdim);
35: PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);
36: PetscSectionSetNumFields(section, 1);
37: PetscSectionSetFieldName(section, field, "U");
38: PetscSectionSetFieldComponents(section, field, sdim);
39: DMPlexGetChart(dm, &pStart, &pEnd);
40: PetscSectionSetChart(section, pStart, pEnd);
42: /* initialize the section storage for a P1 field */
43: DMPlexGetDepthStratum(dm, 0, &pStart, &pEnd);
44: for (p = pStart; p < pEnd; ++p) {
45: PetscSectionSetDof(section, p, sdim);
46: 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: DMGetLabelSize(dm, "Vertex Sets", &numVS);
52: PetscPrintf(PETSC_COMM_WORLD, "# Vertex set: %" PetscInt_FMT "\n", numVS);
53: DMGetLabelIdIS(dm, "Vertex Sets", &setIS);
54: ISGetIndices(setIS, &setID);
55: for (s = 0; s < numVS; ++s) {
56: DMGetStratumIS(dm, "Vertex Sets", setID[s], &pointIS);
57: DMGetStratumSize(dm, "Vertex Sets", setID[s], &numPoints);
58: PetscPrintf(PETSC_COMM_WORLD, "set %" PetscInt_FMT " size: %" PetscInt_FMT "\n", s, numPoints);
59: ISGetIndices(pointIS, &pointID);
60: for (p = 0; p < numPoints; ++p) {
61: PetscPrintf(PETSC_COMM_WORLD, "\t point %" PetscInt_FMT "\n", pointID[p]);
62: PetscSectionSetConstraintDof(section, pointID[p], 1);
63: PetscSectionSetFieldConstraintDof(section, pointID[p], field, 1);
64: }
65: ISRestoreIndices(pointIS, &pointID);
66: ISDestroy(&pointIS);
67: }
69: 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: DMGetStratumIS(dm, "Vertex Sets", setID[s], &pointIS);
75: DMGetStratumSize(dm, "Vertex Sets", setID[s], &numPoints);
76: ISGetIndices(pointIS, &pointID);
77: for (p = 0; p < numPoints; ++p) {
78: constraints[0] = setID[s] % sdim;
79: PetscSectionSetConstraintIndices(section, pointID[p], constraints);
80: PetscSectionSetFieldConstraintIndices(section, pointID[p], field, constraints);
81: }
82: ISRestoreIndices(pointIS, &pointID);
83: ISDestroy(&pointIS);
84: }
85: ISRestoreIndices(setIS, &setID);
86: ISDestroy(&setIS);
87: PetscObjectViewFromOptions((PetscObject)section, NULL, "-dm_section_view");
89: PetscSectionDestroy(§ion);
90: DMDestroy(&dm);
92: 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*/