Actual source code: ex1.c

  1: static char help[] = "Test PCBDDCGraphCreateLocalSubdomainAdjacency() against the DMPlex facet graph.\n\n"
  2:                      "Use -nfields to select the number of fields, with arrays -dofs_type and -p describing their layouts.\n"
  3:                      "Use -bc_dir_dofs and -bc_neu_dofs to mark boundary dofs for each field, and -bc_label to select the boundary label.\n";

  5: #include <petscdmplex.h>
  6: #include <petscdt.h>
  7: #include <petsc/private/pcbddcprivateimpl.h>
  8: #include <petsc/private/pcbddcstructsimpl.h>

 10: typedef enum {
 11:   DOFS_H1,
 12:   DOFS_HCURL,
 13:   DOFS_HDIV
 14: } DofType;

 16: static const char *const DofTypes[] = {"H1", "Hcurl", "Hdiv", "DofType", "DOFS_", NULL};

 18: typedef struct {
 19:   PetscInt  nfields;
 20:   DofType  *dofsTypes;
 21:   PetscInt *p;
 22:   PetscInt *bcDirDofs;
 23:   PetscInt *bcNeuDofs;
 24:   PetscBool bcDirSet;
 25:   PetscBool bcNeuSet;
 26:   char      bcLabel[PETSC_MAX_PATH_LEN];
 27:   PetscBool view_plex_graph;
 28:   PetscBool view_bddc_graph;
 29: } TestOptions;

 31: static PetscErrorCode ComputeNumDof(DM dm, DofType dofsType, PetscInt p, PetscInt numDof[], PetscInt *numComp)
 32: {
 33:   DMPolytopeType cellType, ct;
 34:   PetscInt       dim, cStart, cEnd, formDegree;
 35:   PetscBool      simplex;

 37:   PetscFunctionBeginUser;
 38:   PetscCall(DMGetDimension(dm, &dim));
 39:   PetscCheck(dim == 2 || dim == 3, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only dimensions 2 and 3 are supported, not dimension %" PetscInt_FMT, dim);
 40:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
 41:   PetscCheck(cEnd > cStart, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The mesh has no cells");
 42:   PetscCall(DMPlexGetCellType(dm, cStart, &cellType));
 43:   simplex = (PetscBool)(DMPolytopeTypeGetNumVertices(cellType) == dim + 1);
 44:   PetscCheck(cellType == DMPolytopeTypeSimpleShape(dim, simplex), PETSC_COMM_SELF, PETSC_ERR_SUP, "Only simplex and tensor-product cells are supported");
 45:   for (PetscInt c = cStart + 1; c < cEnd; c++) {
 46:     PetscCall(DMPlexGetCellType(dm, c, &ct));
 47:     PetscCheck(ct == cellType, PETSC_COMM_SELF, PETSC_ERR_SUP, "Meshes with mixed cell types are not supported");
 48:   }

 50:   switch (dofsType) {
 51:   case DOFS_H1:
 52:     formDegree = 0;
 53:     *numComp   = 1;
 54:     break;
 55:   case DOFS_HCURL:
 56:     formDegree = 1;
 57:     *numComp   = dim;
 58:     break;
 59:   case DOFS_HDIV:
 60:     formDegree = dim - 1;
 61:     *numComp   = dim;
 62:     break;
 63:   default:
 64:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown dof type");
 65:   }

 67:   for (PetscInt d = 0; d <= dim; d++) {
 68:     PetscInt value, forms;

 70:     if (d < formDegree) {
 71:       numDof[d] = 0;
 72:       continue;
 73:     }
 74:     PetscCall(PetscDTBinomialInt(d, formDegree, &forms));
 75:     if (simplex) {
 76:       if (d > p + formDegree - 1) {
 77:         numDof[d] = 0;
 78:         continue;
 79:       }
 80:       PetscCall(PetscDTBinomialInt(p + formDegree - 1, d, &value));
 81:       PetscCall(PetscIntMultError(value, forms, &numDof[d]));
 82:     } else {
 83:       value = forms;
 84:       for (PetscInt i = 0; i < formDegree; i++) PetscCall(PetscIntMultError(value, p, &value));
 85:       for (PetscInt i = formDegree; i < d; i++) PetscCall(PetscIntMultError(value, p - 1, &value));
 86:       numDof[d] = value;
 87:     }
 88:   }
 89:   PetscFunctionReturn(PETSC_SUCCESS);
 90: }

 92: static PetscErrorCode ProcessOptions(TestOptions *options)
 93: {
 94:   PetscInt  n;
 95:   PetscBool set;

 97:   PetscFunctionBeginUser;
 98:   options->nfields = 1;
 99:   PetscOptionsBegin(PETSC_COMM_SELF, "", "Section options", "DMPLEX");
100:   PetscCall(PetscOptionsBoundedInt("-nfields", "Number of fields", "ex1.c", 1, &options->nfields, NULL, 1));
101:   PetscCall(PetscMalloc4(options->nfields, &options->dofsTypes, options->nfields, &options->p, options->nfields, &options->bcDirDofs, options->nfields, &options->bcNeuDofs));
102:   for (PetscInt f = 0; f < options->nfields; f++) {
103:     options->dofsTypes[f] = DOFS_H1;
104:     options->p[f]         = 1;
105:     options->bcDirDofs[f] = -1;
106:     options->bcNeuDofs[f] = -1;
107:   }
108:   n = options->nfields;
109:   PetscCall(PetscOptionsEnumArray("-dofs_type", "Finite element de Rham spaces", "ex1.c", DofTypes, (PetscEnum *)options->dofsTypes, &n, &set));
110:   PetscCheck(!set || n == options->nfields, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "-dofs_type requires exactly %" PetscInt_FMT " values, got %" PetscInt_FMT, options->nfields, n);
111:   n = options->nfields;
112:   PetscCall(PetscOptionsIntArray("-p", "Polynomial orders", "ex1.c", options->p, &n, &set));
113:   PetscCheck(!set || n == options->nfields, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "-p requires exactly %" PetscInt_FMT " values, got %" PetscInt_FMT, options->nfields, n);
114:   n = options->nfields;
115:   PetscCall(PetscOptionsIntArray("-bc_dir_dofs", "Dirichlet boundary label values for each field (-1 for the entire boundary)", "ex1.c", options->bcDirDofs, &n, &options->bcDirSet));
116:   PetscCheck(!options->bcDirSet || n == options->nfields, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "-bc_dir_dofs requires exactly %" PetscInt_FMT " values, got %" PetscInt_FMT, options->nfields, n);
117:   n = options->nfields;
118:   PetscCall(PetscOptionsIntArray("-bc_neu_dofs", "Neumann boundary label values for each field (-1 for the entire boundary)", "ex1.c", options->bcNeuDofs, &n, &options->bcNeuSet));
119:   PetscCheck(!options->bcNeuSet || n == options->nfields, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "-bc_neu_dofs requires exactly %" PetscInt_FMT " values, got %" PetscInt_FMT, options->nfields, n);
120:   PetscCall(PetscStrncpy(options->bcLabel, "Face Sets", sizeof(options->bcLabel)));
121:   PetscCall(PetscOptionsString("-bc_label", "Boundary label for nonnegative Dirichlet and Neumann values", "ex1.c", options->bcLabel, options->bcLabel, sizeof(options->bcLabel), NULL));
122:   options->view_plex_graph = PETSC_FALSE;
123:   PetscCall(PetscOptionsBool("-view_plex_graph", "View DMPLEX connectivity graph", "ex1.c", options->view_plex_graph, &options->view_plex_graph, NULL));
124:   options->view_bddc_graph = PETSC_FALSE;
125:   PetscCall(PetscOptionsBool("-view_bddc_graph", "View BDDC connectivity graph", "ex1.c", options->view_bddc_graph, &options->view_bddc_graph, NULL));
126:   PetscOptionsEnd();

128:   for (PetscInt f = 0; f < options->nfields; f++) {
129:     PetscCheck(options->p[f] >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Polynomial order for field %" PetscInt_FMT " must be positive, got %" PetscInt_FMT, f, options->p[f]);
130:     PetscCheck(!options->bcDirSet || options->bcDirDofs[f] >= -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Dirichlet boundary value for field %" PetscInt_FMT " must be at least -1, got %" PetscInt_FMT, f, options->bcDirDofs[f]);
131:     PetscCheck(!options->bcNeuSet || options->bcNeuDofs[f] >= -1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Neumann boundary value for field %" PetscInt_FMT " must be at least -1, got %" PetscInt_FMT, f, options->bcNeuDofs[f]);
132:   }
133:   PetscFunctionReturn(PETSC_SUCCESS);
134: }

136: static PetscErrorCode CreateSection(DM dm, const TestOptions *options, PetscSection *section)
137: {
138:   PetscInt *numComp, *numDof;
139:   PetscInt  dim;

141:   PetscFunctionBeginUser;
142:   PetscCall(DMGetDimension(dm, &dim));
143:   PetscCall(PetscMalloc2(options->nfields, &numComp, options->nfields * (dim + 1), &numDof));
144:   for (PetscInt f = 0; f < options->nfields; f++) PetscCall(ComputeNumDof(dm, options->dofsTypes[f], options->p[f], &numDof[f * (dim + 1)], &numComp[f]));
145:   PetscCall(DMSetNumFields(dm, options->nfields));
146:   PetscCall(DMPlexCreateSection(dm, NULL, numComp, numDof, 0, NULL, NULL, NULL, NULL, section));
147:   PetscCall(PetscFree2(numComp, numDof));
148:   PetscFunctionReturn(PETSC_SUCCESS);
149: }

151: static PetscErrorCode CreateElementMapping(DM dm, PetscSection section, PetscInt nfields, PetscInt *numCells, PetscInt *numDofs, ISLocalToGlobalMapping *l2g, PetscInt **localSubs, IS *fieldIS[])
152: {
153:   PetscInt  *indices, *l2gIndices, *subs, *globalField, *fieldCounts, *fieldCursor;
154:   PetscInt **fieldIndices;
155:   PetscInt   cStart, cEnd, pStart, pEnd, c, n, nlocal = 0, off = 0;

157:   PetscFunctionBeginUser;
158:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
159:   for (c = cStart; c < cEnd; c++) {
160:     PetscCall(DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, &n, &indices, NULL, NULL));
161:     nlocal += n;
162:     PetscCall(DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, &n, &indices, NULL, NULL));
163:   }
164:   PetscCall(PetscSectionGetStorageSize(section, numDofs));
165:   PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
166:   PetscCall(PetscMalloc1(nlocal, &l2gIndices));
167:   PetscCall(PetscMalloc1(nlocal, &subs));
168:   PetscCall(PetscMalloc1(*numDofs, &globalField));
169:   for (PetscInt i = 0; i < *numDofs; i++) globalField[i] = -1;
170:   for (PetscInt p = pStart; p < pEnd; p++) {
171:     for (PetscInt f = 0; f < nfields; f++) {
172:       PetscInt fdof, foff;

174:       PetscCall(PetscSectionGetFieldDof(section, p, f, &fdof));
175:       PetscCall(PetscSectionGetFieldOffset(section, p, f, &foff));
176:       for (PetscInt i = 0; i < fdof; i++) {
177:         PetscCheck(foff + i >= 0 && foff + i < *numDofs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid field offset %" PetscInt_FMT " not in [0,%" PetscInt_FMT ")", foff + i, *numDofs);
178:         globalField[foff + i] = f;
179:       }
180:     }
181:   }
182:   for (c = cStart; c < cEnd; c++) {
183:     PetscCall(DMPlexGetClosureIndices(dm, section, section, c, PETSC_TRUE, &n, &indices, NULL, NULL));
184:     for (PetscInt i = 0; i < n; i++) {
185:       PetscCheck(indices[i] >= 0 && indices[i] < *numDofs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid section index %" PetscInt_FMT " not in [0,%" PetscInt_FMT ")", indices[i], *numDofs);
186:       PetscCheck(globalField[indices[i]] >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Section index %" PetscInt_FMT " does not belong to a field", indices[i]);
187:       l2gIndices[off] = indices[i];
188:       subs[off++]     = c - cStart;
189:     }
190:     PetscCall(DMPlexRestoreClosureIndices(dm, section, section, c, PETSC_TRUE, &n, &indices, NULL, NULL));
191:   }
192:   PetscCheck(off == nlocal, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Mapped %" PetscInt_FMT " dofs, expected %" PetscInt_FMT, off, nlocal);
193:   PetscCall(PetscCalloc2(nfields, &fieldCounts, nfields, &fieldCursor));
194:   PetscCall(PetscMalloc1(nfields, &fieldIndices));
195:   PetscCall(PetscMalloc1(nfields, fieldIS));
196:   for (PetscInt i = 0; i < nlocal; i++) fieldCounts[globalField[l2gIndices[i]]]++;
197:   for (PetscInt f = 0; f < nfields; f++) PetscCall(PetscMalloc1(fieldCounts[f], &fieldIndices[f]));
198:   for (PetscInt i = 0; i < nlocal; i++) {
199:     const PetscInt f = globalField[l2gIndices[i]];

201:     fieldIndices[f][fieldCursor[f]++] = i;
202:   }
203:   for (PetscInt f = 0; f < nfields; f++) PetscCall(ISCreateGeneral(PETSC_COMM_SELF, fieldCounts[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fieldIS)[f]));
204:   PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, 1, nlocal, l2gIndices, PETSC_OWN_POINTER, l2g));
205:   PetscCall(PetscFree(fieldIndices));
206:   PetscCall(PetscFree2(fieldCounts, fieldCursor));
207:   PetscCall(PetscFree(globalField));
208:   *numCells  = cEnd - cStart;
209:   *localSubs = subs;
210:   PetscFunctionReturn(PETSC_SUCCESS);
211: }

213: static PetscErrorCode CreateBoundaryIS(DM dm, PetscSection section, PetscInt nfields, const PetscInt boundaryValues[], const char boundaryLabelName[], ISLocalToGlobalMapping l2g, const char optionName[], IS *boundaryIS)
214: {
215:   DMLabel         boundaryLabel;
216:   const PetscInt *l2gIndices;
217:   PetscInt       *indices;
218:   PetscInt        pStart, pEnd, numDofs, nlocal, n = 0;
219:   PetscBool      *marked;

221:   PetscFunctionBeginUser;
222:   PetscCall(DMGetLabel(dm, boundaryLabelName, &boundaryLabel));
223:   PetscCheck(boundaryLabel, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "%s requires the mesh label '%s'", optionName, boundaryLabelName);
224:   PetscCall(DMPlexLabelComplete(dm, boundaryLabel));
225:   for (PetscInt f = 0; f < nfields; f++) {
226:     if (boundaryValues[f] >= 0) {
227:       PetscInt labelSize;

229:       PetscCall(DMLabelGetStratumSize(boundaryLabel, boundaryValues[f], &labelSize));
230:       PetscCheck(labelSize, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "%s value %" PetscInt_FMT " for field %" PetscInt_FMT " is not present in '%s'", optionName, boundaryValues[f], f, boundaryLabelName);
231:     }
232:   }

234:   PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
235:   PetscCall(PetscSectionGetStorageSize(section, &numDofs));
236:   PetscCall(PetscCalloc1(numDofs, &marked));
237:   for (PetscInt p = pStart; p < pEnd; p++) {
238:     for (PetscInt f = 0; f < nfields; f++) {
239:       PetscInt  fdof, foff;
240:       PetscBool hasPoint;

242:       if (boundaryValues[f] == -1) PetscCall(DMLabelHasPoint(boundaryLabel, p, &hasPoint));
243:       else PetscCall(DMLabelStratumHasPoint(boundaryLabel, boundaryValues[f], p, &hasPoint));
244:       if (!hasPoint) continue;
245:       PetscCall(PetscSectionGetFieldDof(section, p, f, &fdof));
246:       PetscCall(PetscSectionGetFieldOffset(section, p, f, &foff));
247:       for (PetscInt i = 0; i < fdof; i++) marked[foff + i] = PETSC_TRUE;
248:     }
249:   }

251:   PetscCall(ISLocalToGlobalMappingGetSize(l2g, &nlocal));
252:   PetscCall(ISLocalToGlobalMappingGetIndices(l2g, &l2gIndices));
253:   for (PetscInt i = 0; i < nlocal; i++)
254:     if (marked[l2gIndices[i]]) n++;
255:   PetscCall(PetscMalloc1(n, &indices));
256:   for (PetscInt i = 0, j = 0; i < nlocal; i++)
257:     if (marked[l2gIndices[i]]) indices[j++] = i;
258:   PetscCall(ISLocalToGlobalMappingRestoreIndices(l2g, &l2gIndices));
259:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, indices, PETSC_OWN_POINTER, boundaryIS));
260:   PetscCall(PetscFree(marked));
261:   PetscFunctionReturn(PETSC_SUCCESS);
262: }

264: static PetscErrorCode CreateReferenceGraph(DM dm, PetscInt *numCells, PetscInt **xadj, PetscInt **adjncy)
265: {
266:   PetscFunctionBeginUser;
267:   PetscCall(DMPlexCreatePartitionerGraph(dm, 0, numCells, xadj, adjncy, NULL));
268:   for (PetscInt c = 0; c < *numCells; c++) PetscCall(PetscSortInt((*xadj)[c + 1] - (*xadj)[c], *adjncy + (*xadj)[c]));
269:   PetscFunctionReturn(PETSC_SUCCESS);
270: }

272: static PetscErrorCode CompareGraphs(PetscInt numCells, const PetscInt plexXadj[], const PetscInt plexAdjncy[], const PetscInt bddcXadj[], const PetscInt bddcAdjncy[])
273: {
274:   PetscFunctionBeginUser;
275:   for (PetscInt c = 0; c < numCells; c++) {
276:     const PetscInt plexDegree = plexXadj[c + 1] - plexXadj[c];
277:     const PetscInt bddcDegree = bddcXadj[c + 1] - bddcXadj[c];

279:     PetscCheck(plexDegree == bddcDegree, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %" PetscInt_FMT " has DMPlex degree %" PetscInt_FMT " and PCBDDCGraph degree %" PetscInt_FMT, c, plexDegree, bddcDegree);
280:     for (PetscInt i = 0; i < plexDegree; i++)
281:       PetscCheck(plexAdjncy[plexXadj[c] + i] == bddcAdjncy[bddcXadj[c] + i], PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cell %" PetscInt_FMT " adjacency %" PetscInt_FMT " is %" PetscInt_FMT " in DMPlex and %" PetscInt_FMT " in PCBDDCGraph", c, i, plexAdjncy[plexXadj[c] + i], bddcAdjncy[bddcXadj[c] + i]);
282:   }
283:   PetscFunctionReturn(PETSC_SUCCESS);
284: }

286: int main(int argc, char **argv)
287: {
288:   PCBDDCGraph            graph;
289:   DM                     dm;
290:   TestOptions            options = {0};
291:   IS                    *fieldIS;
292:   IS                     dirichletIS = NULL, neumannIS = NULL;
293:   ISLocalToGlobalMapping l2g;
294:   PetscSection           section;
295:   PetscInt              *plexXadj, *plexAdjncy, *bddcXadj, *bddcAdjncy, *localSubs = NULL;
296:   PetscInt               dim, depth, numCells = 0, graphCells, bddcCells, numDofs;

298:   PetscFunctionBeginUser;
299:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
300:   PetscCall(DMCreate(PETSC_COMM_SELF, &dm));
301:   PetscCall(DMSetType(dm, DMPLEX));
302:   PetscCall(DMSetFromOptions(dm));
303:   PetscCall(DMViewFromOptions(dm, NULL, "-mesh_view"));
304:   PetscCall(DMGetDimension(dm, &dim));
305:   PetscCall(DMPlexGetDepth(dm, &depth));
306:   PetscCheck(depth == dim, PETSC_COMM_SELF, PETSC_ERR_SUP, "The mesh must be fully interpolated so that DMPlex facet connectivity is explicit");

308:   PetscCall(ProcessOptions(&options));
309:   PetscCall(CreateSection(dm, &options, &section));
310:   PetscCall(CreateElementMapping(dm, section, options.nfields, &numCells, &numDofs, &l2g, &localSubs, &fieldIS));
311:   if (options.bcDirSet) PetscCall(CreateBoundaryIS(dm, section, options.nfields, options.bcDirDofs, options.bcLabel, l2g, "-bc_dir_dofs", &dirichletIS));
312:   if (options.bcNeuSet) PetscCall(CreateBoundaryIS(dm, section, options.nfields, options.bcNeuDofs, options.bcLabel, l2g, "-bc_neu_dofs", &neumannIS));
313:   PetscCall(CreateReferenceGraph(dm, &graphCells, &plexXadj, &plexAdjncy));
314:   PetscCheck(graphCells == numCells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "DMPlex graph has %" PetscInt_FMT " cells, element map has %" PetscInt_FMT, graphCells, numCells);

316:   PetscCall(PCBDDCGraphCreate(&graph));
317:   PetscCall(PCBDDCGraphInit(graph, l2g, numDofs, PETSC_INT_MAX));
318:   graph->n_local_subs = numCells;
319:   graph->local_subs   = localSubs;
320:   PetscCall(PCBDDCGraphSetUp(graph, 1, neumannIS, dirichletIS, options.nfields, fieldIS, NULL));
321:   PetscCall(PCBDDCGraphComputeConnectedComponents(graph));
322:   PetscCall(PCBDDCGraphCreateLocalSubdomainAdjacency(graph, &bddcCells, &bddcXadj, &bddcAdjncy));
323:   PetscCheck(bddcCells == numCells, PETSC_COMM_SELF, PETSC_ERR_PLIB, "PCBDDC graph has %" PetscInt_FMT " cells, element map has %" PetscInt_FMT, bddcCells, numCells);
324:   if (options.view_plex_graph) PetscCall(PetscIntCSRView(numCells, plexXadj, plexAdjncy, NULL));
325:   if (options.view_bddc_graph) PetscCall(PetscIntCSRView(bddcCells, bddcXadj, bddcAdjncy, NULL));
326:   PetscCall(CompareGraphs(numCells, plexXadj, plexAdjncy, bddcXadj, bddcAdjncy));
327:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "Graphs match\n"));

329:   PetscCall(PetscFree(bddcXadj));
330:   PetscCall(PetscFree(bddcAdjncy));
331:   PetscCall(PetscFree(plexXadj));
332:   PetscCall(PetscFree(plexAdjncy));
333:   PetscCall(PCBDDCGraphDestroy(&graph));
334:   for (PetscInt f = 0; f < options.nfields; f++) PetscCall(ISDestroy(&fieldIS[f]));
335:   PetscCall(PetscFree(fieldIS));
336:   PetscCall(ISDestroy(&dirichletIS));
337:   PetscCall(ISDestroy(&neumannIS));
338:   PetscCall(ISLocalToGlobalMappingDestroy(&l2g));
339:   PetscCall(PetscSectionDestroy(&section));
340:   PetscCall(PetscFree4(options.dofsTypes, options.p, options.bcDirDofs, options.bcNeuDofs));
341:   PetscCall(DMDestroy(&dm));
342:   PetscCall(PetscFinalize());
343:   return 0;
344: }

346: /*TEST

348:   testset:
349:     nsize: 1
350:     output_file: output/ex1.out
351:     args: -dm_plex_interpolate 1 -dm_plex_csr_alg graph

353:     test:
354:       requires: triangle
355:       suffix: h1_2d
356:       args: -dm_plex_dim 2 -dm_plex_simplex {{0 1}} -dm_plex_box_faces 3,3 -dofs_type H1 -p {{1 2 3}}

358:     test:
359:       requires: triangle
360:       suffix: hcurl_2d
361:       args: -dm_plex_dim 2 -dm_plex_simplex {{0 1}} -dm_plex_box_faces 3,3 -dofs_type Hcurl -p {{1 2 3}}

363:     test:
364:       requires: triangle
365:       suffix: hdiv_2d
366:       args: -dm_plex_dim 2 -dm_plex_simplex {{0 1}} -dm_plex_box_faces 3,3 -dofs_type Hdiv -p {{1 2 3}}

368:     test:
369:       requires: ctetgen
370:       suffix: h1_3d
371:       args: -dm_plex_dim 3 -dm_plex_simplex {{0 1}} -dm_plex_box_faces 3,3,3 -dofs_type H1 -p {{1 2 3}}

373:     test:
374:       requires: ctetgen
375:       suffix: hcurl_3d
376:       args: -dm_plex_dim 3 -dm_plex_simplex {{0 1}} -dm_plex_box_faces 3,3,3 -dofs_type Hcurl -p {{1 2 3}}

378:     test:
379:       requires: ctetgen
380:       suffix: hdiv_3d
381:       args: -dm_plex_dim 3 -dm_plex_simplex {{0 1}} -dm_plex_box_faces 3,3,3 -dofs_type Hdiv -p {{1 2 3}}

383:     test:
384:       requires: triangle
385:       suffix: multifield_2d
386:       args: -dm_plex_dim 2 -dm_plex_simplex {{0 1}} -dm_plex_box_faces 3,3 -nfields 3 -dofs_type H1,Hcurl,Hdiv -p 1,2,3

388:     test:
389:       requires: ctetgen
390:       suffix: multifield_3d
391:       args: -dm_plex_dim 3 -dm_plex_simplex {{0 1}} -dm_plex_box_faces 3,3,3 -nfields 3 -dofs_type H1,Hcurl,Hdiv -p 1,2,3

393:     test:
394:       requires: triangle
395:       suffix: dirichlet_all_2d
396:       args: -dm_plex_dim 2 -dm_plex_simplex {{0 1}} -dm_plex_box_faces 3,3 -dm_plex_box_label -dofs_type H1 -p {{1 2 3}} -bc_dir_dofs -1

398:     test:
399:       requires: ctetgen
400:       suffix: dirichlet_all_3d
401:       args: -dm_plex_dim 3 -dm_plex_simplex {{0 1}} -dm_plex_box_faces 3,3,3 -dm_plex_box_label -dofs_type H1 -p {{1 2 3}} -bc_dir_dofs -1

403:     test:
404:       requires: triangle
405:       suffix: boundary_fields_2d
406:       args: -dm_plex_dim 2 -dm_plex_simplex {{0 1}} -dm_plex_box_faces 3,3 -dm_plex_box_label -nfields 3 -dofs_type H1,Hcurl,Hdiv -p 1,2,3 -bc_dir_dofs -1,1,2 -bc_neu_dofs 3,-1,4

408:     test:
409:       requires: ctetgen
410:       suffix: boundary_fields_3d
411:       args: -dm_plex_dim 3 -dm_plex_simplex {{0 1}} -dm_plex_box_faces 3,3,3 -dm_plex_box_label -nfields 3 -dofs_type H1,Hcurl,Hdiv -p 1,2,3 -bc_dir_dofs -1,1,2 -bc_neu_dofs 4,-1,6

413:     test:
414:       requires: triangle
415:       suffix: boundary_label
416:       args: -dm_plex_dim 2 -dm_plex_simplex 0 -dm_plex_box_faces 3,3 -dm_plex_boundary_label boundary -bc_label boundary -bc_dir_dofs 1

418: TEST*/