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), &section));
 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(&section));
 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*/