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