1 static char help[] = "Test FEM layout with constraints\n\n"; 2 3 #include <petsc.h> 4 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