1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 /* Set the number of dof on each point and separate by fields */ 4 static PetscErrorCode DMPlexCreateSectionFields(DM dm, const PetscInt numComp[], PetscSection *section) 5 { 6 DMLabel depthLabel; 7 PetscInt depth, Nf, f, pStart, pEnd; 8 PetscBool *isFE; 9 PetscErrorCode ierr; 10 11 PetscFunctionBegin; 12 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 13 ierr = DMPlexGetDepthLabel(dm,&depthLabel);CHKERRQ(ierr); 14 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 15 ierr = PetscCalloc1(Nf, &isFE);CHKERRQ(ierr); 16 for (f = 0; f < Nf; ++f) { 17 PetscObject obj; 18 PetscClassId id; 19 20 ierr = DMGetField(dm, f, NULL, &obj);CHKERRQ(ierr); 21 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 22 if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;} 23 else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;} 24 } 25 26 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), section);CHKERRQ(ierr); 27 if (Nf > 0) { 28 ierr = PetscSectionSetNumFields(*section, Nf);CHKERRQ(ierr); 29 if (numComp) { 30 for (f = 0; f < Nf; ++f) { 31 ierr = PetscSectionSetFieldComponents(*section, f, numComp[f]);CHKERRQ(ierr); 32 if (isFE[f]) { 33 PetscFE fe; 34 PetscDualSpace dspace; 35 const PetscInt ***perms; 36 const PetscScalar ***flips; 37 const PetscInt *numDof; 38 39 ierr = DMGetField(dm,f,NULL,(PetscObject *) &fe);CHKERRQ(ierr); 40 ierr = PetscFEGetDualSpace(fe,&dspace);CHKERRQ(ierr); 41 ierr = PetscDualSpaceGetSymmetries(dspace,&perms,&flips);CHKERRQ(ierr); 42 ierr = PetscDualSpaceGetNumDof(dspace,&numDof);CHKERRQ(ierr); 43 if (perms || flips) { 44 DM K; 45 PetscInt h; 46 PetscSectionSym sym; 47 48 ierr = PetscDualSpaceGetDM(dspace,&K);CHKERRQ(ierr); 49 ierr = PetscSectionSymCreateLabel(PetscObjectComm((PetscObject)*section),depthLabel,&sym);CHKERRQ(ierr); 50 for (h = 0; h <= depth; h++) { 51 PetscDualSpace hspace; 52 PetscInt kStart, kEnd; 53 PetscInt kConeSize; 54 const PetscInt **perms0 = NULL; 55 const PetscScalar **flips0 = NULL; 56 57 ierr = PetscDualSpaceGetHeightSubspace(dspace,h,&hspace);CHKERRQ(ierr); 58 ierr = DMPlexGetHeightStratum(K,h,&kStart,&kEnd);CHKERRQ(ierr); 59 if (!hspace) continue; 60 ierr = PetscDualSpaceGetSymmetries(hspace,&perms,&flips);CHKERRQ(ierr); 61 if (perms) perms0 = perms[0]; 62 if (flips) flips0 = flips[0]; 63 if (!(perms0 || flips0)) continue; 64 ierr = DMPlexGetConeSize(K,kStart,&kConeSize);CHKERRQ(ierr); 65 ierr = PetscSectionSymLabelSetStratum(sym,depth - h,numDof[depth - h],-kConeSize,kConeSize,PETSC_USE_POINTER,perms0 ? &perms0[-kConeSize] : NULL,flips0 ? &flips0[-kConeSize] : NULL);CHKERRQ(ierr); 66 } 67 ierr = PetscSectionSetFieldSym(*section,f,sym);CHKERRQ(ierr); 68 ierr = PetscSectionSymDestroy(&sym);CHKERRQ(ierr); 69 } 70 } 71 } 72 } 73 } 74 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 75 ierr = PetscSectionSetChart(*section, pStart, pEnd);CHKERRQ(ierr); 76 ierr = PetscFree(isFE);CHKERRQ(ierr); 77 PetscFunctionReturn(0); 78 } 79 80 /* Set the number of dof on each point and separate by fields */ 81 static PetscErrorCode DMPlexCreateSectionDof(DM dm, DMLabel label[],const PetscInt numDof[], PetscSection section) 82 { 83 DMLabel depthLabel; 84 PetscInt *pMax; 85 PetscInt depth, cellHeight, pStart = 0, pEnd = 0; 86 PetscInt Nf, f, dim, d, dep, p; 87 PetscBool *isFE; 88 PetscErrorCode ierr; 89 90 PetscFunctionBegin; 91 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 92 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 93 ierr = DMPlexGetDepthLabel(dm,&depthLabel);CHKERRQ(ierr); 94 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 95 ierr = PetscMalloc1(Nf, &isFE);CHKERRQ(ierr); 96 for (f = 0; f < Nf; ++f) { 97 PetscObject obj; 98 PetscClassId id; 99 100 ierr = DMGetField(dm, f, NULL, &obj);CHKERRQ(ierr); 101 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 102 /* User is allowed to put a "placeholder" field in (c.f. DMCreateDS) */ 103 isFE[f] = id == PETSCFE_CLASSID ? PETSC_TRUE : PETSC_FALSE; 104 } 105 106 ierr = PetscMalloc1(depth+1, &pMax);CHKERRQ(ierr); 107 ierr = DMPlexGetHybridBounds(dm, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);CHKERRQ(ierr); 108 ierr = DMPlexGetVTKCellHeight(dm, &cellHeight);CHKERRQ(ierr); 109 for (f = 0; f < Nf; ++f) { 110 if (label && label[f]) { 111 IS pointIS; 112 const PetscInt *points; 113 PetscInt n; 114 115 ierr = DMLabelGetStratumIS(label[f], 1, &pointIS);CHKERRQ(ierr); 116 if (!pointIS) continue; 117 ierr = ISGetLocalSize(pointIS, &n);CHKERRQ(ierr); 118 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 119 for (p = 0; p < n; ++p) { 120 const PetscInt point = points[p]; 121 PetscInt dof, d; 122 123 ierr = DMLabelGetValue(depthLabel, point, &d);CHKERRQ(ierr); 124 ierr = DMPlexGetDepthStratum(dm, d, NULL, &pEnd);CHKERRQ(ierr); 125 /* If this is a hybrid point, use dof for one dimension lower */ 126 if ((point >= pMax[d]) && (point < pEnd)) --d; 127 dof = d < 0 ? 0 : numDof[f*(dim+1)+d]; 128 ierr = PetscSectionSetFieldDof(section, point, f, dof);CHKERRQ(ierr); 129 ierr = PetscSectionAddDof(section, point, dof);CHKERRQ(ierr); 130 } 131 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 132 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 133 } else { 134 for (dep = 0; dep <= depth - cellHeight; ++dep) { 135 /* Cases: dim > depth (cell-vertex mesh), dim == depth (fully interpolated), dim < depth (interpolated submesh) */ 136 d = dim <= depth ? dep : (!dep ? 0 : dim); 137 ierr = DMPlexGetDepthStratum(dm, dep, &pStart, &pEnd);CHKERRQ(ierr); 138 pMax[dep] = pMax[dep] < 0 ? pEnd : pMax[dep]; 139 for (p = pStart; p < pEnd; ++p) { 140 const PetscInt dof = numDof[f*(dim+1)+d]; 141 142 if (isFE[f] && p >= pMax[dep]) continue; 143 ierr = PetscSectionSetFieldDof(section, p, f, dof);CHKERRQ(ierr); 144 ierr = PetscSectionAddDof(section, p, dof);CHKERRQ(ierr); 145 } 146 } 147 } 148 } 149 ierr = PetscFree(pMax);CHKERRQ(ierr); 150 ierr = PetscFree(isFE);CHKERRQ(ierr); 151 PetscFunctionReturn(0); 152 } 153 154 /* Set the number of dof on each point and separate by fields 155 If bcComps is NULL or the IS is NULL, constrain every dof on the point 156 */ 157 static PetscErrorCode DMPlexCreateSectionBCDof(DM dm, PetscInt numBC, const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section) 158 { 159 PetscInt Nf; 160 PetscInt bc; 161 PetscSection aSec; 162 PetscErrorCode ierr; 163 164 PetscFunctionBegin; 165 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 166 for (bc = 0; bc < numBC; ++bc) { 167 PetscInt field = 0; 168 const PetscInt *comp; 169 const PetscInt *idx; 170 PetscInt Nc = 0, cNc = -1, n, i; 171 172 if (Nf) { 173 field = bcField[bc]; 174 ierr = PetscSectionGetFieldComponents(section, field, &Nc);CHKERRQ(ierr); 175 } 176 if (bcComps && bcComps[bc]) {ierr = ISGetLocalSize(bcComps[bc], &cNc);CHKERRQ(ierr);} 177 if (bcComps && bcComps[bc]) {ierr = ISGetIndices(bcComps[bc], &comp);CHKERRQ(ierr);} 178 ierr = ISGetLocalSize(bcPoints[bc], &n);CHKERRQ(ierr); 179 ierr = ISGetIndices(bcPoints[bc], &idx);CHKERRQ(ierr); 180 for (i = 0; i < n; ++i) { 181 const PetscInt p = idx[i]; 182 PetscInt numConst; 183 184 if (Nf) { 185 ierr = PetscSectionGetFieldDof(section, p, field, &numConst);CHKERRQ(ierr); 186 } else { 187 ierr = PetscSectionGetDof(section, p, &numConst);CHKERRQ(ierr); 188 } 189 /* If Nc <= 0, constrain every dof on the point */ 190 if (cNc > 0) { 191 /* We assume that a point may have multiple "nodes", which are collections of Nc dofs, 192 and that those dofs are numbered n*Nc+c */ 193 if (Nf) { 194 if (numConst % Nc) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %D has %D dof which is not divisible by %D field components", p, numConst, Nc); 195 numConst = (numConst/Nc) * cNc; 196 } else { 197 numConst = PetscMin(numConst, cNc); 198 } 199 } 200 if (Nf) {ierr = PetscSectionAddFieldConstraintDof(section, p, field, numConst);CHKERRQ(ierr);} 201 ierr = PetscSectionAddConstraintDof(section, p, numConst);CHKERRQ(ierr); 202 } 203 ierr = ISRestoreIndices(bcPoints[bc], &idx);CHKERRQ(ierr); 204 if (bcComps && bcComps[bc]) {ierr = ISRestoreIndices(bcComps[bc], &comp);CHKERRQ(ierr);} 205 } 206 ierr = DMPlexGetAnchors(dm, &aSec, NULL);CHKERRQ(ierr); 207 if (aSec) { 208 PetscInt aStart, aEnd, a; 209 210 ierr = PetscSectionGetChart(aSec, &aStart, &aEnd);CHKERRQ(ierr); 211 for (a = aStart; a < aEnd; a++) { 212 PetscInt dof, f; 213 214 ierr = PetscSectionGetDof(aSec, a, &dof);CHKERRQ(ierr); 215 if (dof) { 216 /* if there are point-to-point constraints, then all dofs are constrained */ 217 ierr = PetscSectionGetDof(section, a, &dof);CHKERRQ(ierr); 218 ierr = PetscSectionSetConstraintDof(section, a, dof);CHKERRQ(ierr); 219 for (f = 0; f < Nf; f++) { 220 ierr = PetscSectionGetFieldDof(section, a, f, &dof);CHKERRQ(ierr); 221 ierr = PetscSectionSetFieldConstraintDof(section, a, f, dof);CHKERRQ(ierr); 222 } 223 } 224 } 225 } 226 PetscFunctionReturn(0); 227 } 228 229 /* Set the constrained field indices on each point 230 If bcComps is NULL or the IS is NULL, constrain every dof on the point 231 */ 232 static PetscErrorCode DMPlexCreateSectionBCIndicesField(DM dm, PetscInt numBC,const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], PetscSection section) 233 { 234 PetscSection aSec; 235 PetscInt *indices; 236 PetscInt Nf, cdof, maxDof = 0, pStart, pEnd, p, bc, f, d; 237 PetscErrorCode ierr; 238 239 PetscFunctionBegin; 240 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 241 if (!Nf) PetscFunctionReturn(0); 242 /* Initialize all field indices to -1 */ 243 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr); 244 for (p = pStart; p < pEnd; ++p) {ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr); maxDof = PetscMax(maxDof, cdof);} 245 ierr = PetscMalloc1(maxDof, &indices);CHKERRQ(ierr); 246 for (d = 0; d < maxDof; ++d) indices[d] = -1; 247 for (p = pStart; p < pEnd; ++p) for (f = 0; f < Nf; ++f) {ierr = PetscSectionSetFieldConstraintIndices(section, p, f, indices);CHKERRQ(ierr);} 248 /* Handle BC constraints */ 249 for (bc = 0; bc < numBC; ++bc) { 250 const PetscInt field = bcField[bc]; 251 const PetscInt *comp, *idx; 252 PetscInt Nc, cNc = -1, n, i; 253 254 ierr = PetscSectionGetFieldComponents(section, field, &Nc);CHKERRQ(ierr); 255 if (bcComps && bcComps[bc]) {ierr = ISGetLocalSize(bcComps[bc], &cNc);CHKERRQ(ierr);} 256 if (bcComps && bcComps[bc]) {ierr = ISGetIndices(bcComps[bc], &comp);CHKERRQ(ierr);} 257 ierr = ISGetLocalSize(bcPoints[bc], &n);CHKERRQ(ierr); 258 ierr = ISGetIndices(bcPoints[bc], &idx);CHKERRQ(ierr); 259 for (i = 0; i < n; ++i) { 260 const PetscInt p = idx[i]; 261 const PetscInt *find; 262 PetscInt fdof, fcdof, c, j; 263 264 ierr = PetscSectionGetFieldDof(section, p, field, &fdof);CHKERRQ(ierr); 265 if (!fdof) continue; 266 if (cNc < 0) { 267 for (d = 0; d < fdof; ++d) indices[d] = d; 268 fcdof = fdof; 269 } else { 270 /* We assume that a point may have multiple "nodes", which are collections of Nc dofs, 271 and that those dofs are numbered n*Nc+c */ 272 ierr = PetscSectionGetFieldConstraintDof(section, p, field, &fcdof);CHKERRQ(ierr); 273 ierr = PetscSectionGetFieldConstraintIndices(section, p, field, &find);CHKERRQ(ierr); 274 /* Get indices constrained by previous bcs */ 275 for (d = 0; d < fcdof; ++d) {if (find[d] < 0) break; indices[d] = find[d];} 276 for (j = 0; j < fdof/Nc; ++j) for (c = 0; c < cNc; ++c) indices[d++] = j*Nc + comp[c]; 277 ierr = PetscSortRemoveDupsInt(&d, indices);CHKERRQ(ierr); 278 for (c = d; c < fcdof; ++c) indices[c] = -1; 279 fcdof = d; 280 } 281 ierr = PetscSectionSetFieldConstraintDof(section, p, field, fcdof);CHKERRQ(ierr); 282 ierr = PetscSectionSetFieldConstraintIndices(section, p, field, indices);CHKERRQ(ierr); 283 } 284 if (bcComps && bcComps[bc]) {ierr = ISRestoreIndices(bcComps[bc], &comp);CHKERRQ(ierr);} 285 ierr = ISRestoreIndices(bcPoints[bc], &idx);CHKERRQ(ierr); 286 } 287 /* Handle anchors */ 288 ierr = DMPlexGetAnchors(dm, &aSec, NULL);CHKERRQ(ierr); 289 if (aSec) { 290 PetscInt aStart, aEnd, a; 291 292 for (d = 0; d < maxDof; ++d) indices[d] = d; 293 ierr = PetscSectionGetChart(aSec, &aStart, &aEnd);CHKERRQ(ierr); 294 for (a = aStart; a < aEnd; a++) { 295 PetscInt dof, f; 296 297 ierr = PetscSectionGetDof(aSec, a, &dof);CHKERRQ(ierr); 298 if (dof) { 299 /* if there are point-to-point constraints, then all dofs are constrained */ 300 for (f = 0; f < Nf; f++) { 301 ierr = PetscSectionSetFieldConstraintIndices(section, a, f, indices);CHKERRQ(ierr); 302 } 303 } 304 } 305 } 306 ierr = PetscFree(indices);CHKERRQ(ierr); 307 PetscFunctionReturn(0); 308 } 309 310 /* Set the constrained indices on each point */ 311 static PetscErrorCode DMPlexCreateSectionBCIndices(DM dm, PetscSection section) 312 { 313 PetscInt *indices; 314 PetscInt Nf, maxDof, pStart, pEnd, p, f, d; 315 PetscErrorCode ierr; 316 317 PetscFunctionBegin; 318 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 319 ierr = PetscSectionGetMaxDof(section, &maxDof);CHKERRQ(ierr); 320 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr); 321 ierr = PetscMalloc1(maxDof, &indices);CHKERRQ(ierr); 322 for (d = 0; d < maxDof; ++d) indices[d] = -1; 323 for (p = pStart; p < pEnd; ++p) { 324 PetscInt cdof, d; 325 326 ierr = PetscSectionGetConstraintDof(section, p, &cdof);CHKERRQ(ierr); 327 if (cdof) { 328 if (Nf) { 329 PetscInt numConst = 0, foff = 0; 330 331 for (f = 0; f < Nf; ++f) { 332 const PetscInt *find; 333 PetscInt fcdof, fdof; 334 335 ierr = PetscSectionGetFieldDof(section, p, f, &fdof);CHKERRQ(ierr); 336 ierr = PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);CHKERRQ(ierr); 337 /* Change constraint numbering from field component to local dof number */ 338 ierr = PetscSectionGetFieldConstraintIndices(section, p, f, &find);CHKERRQ(ierr); 339 for (d = 0; d < fcdof; ++d) indices[numConst+d] = find[d] + foff; 340 numConst += fcdof; 341 foff += fdof; 342 } 343 if (cdof != numConst) {ierr = PetscSectionSetConstraintDof(section, p, numConst);CHKERRQ(ierr);} 344 } else { 345 for (d = 0; d < cdof; ++d) indices[d] = d; 346 } 347 ierr = PetscSectionSetConstraintIndices(section, p, indices);CHKERRQ(ierr); 348 } 349 } 350 ierr = PetscFree(indices);CHKERRQ(ierr); 351 PetscFunctionReturn(0); 352 } 353 354 /*@C 355 DMPlexCreateSection - Create a PetscSection based upon the dof layout specification provided. 356 357 Not Collective 358 359 Input Parameters: 360 + dm - The DMPlex object 361 . label - The label indicating the mesh support of each field, or NULL for the whole mesh 362 . numComp - An array of size numFields that holds the number of components for each field 363 . numDof - An array of size numFields*(dim+1) which holds the number of dof for each field on a mesh piece of dimension d 364 . numBC - The number of boundary conditions 365 . bcField - An array of size numBC giving the field number for each boundry condition 366 . bcComps - [Optional] An array of size numBC giving an IS holding the field components to which each boundary condition applies 367 . bcPoints - An array of size numBC giving an IS holding the Plex points to which each boundary condition applies 368 - perm - Optional permutation of the chart, or NULL 369 370 Output Parameter: 371 . section - The PetscSection object 372 373 Notes: 374 numDof[f*(dim+1)+d] gives the number of dof for field f on points of dimension d. For instance, numDof[1] is the 375 number of dof for field 0 on each edge. 376 377 The chart permutation is the same one set using PetscSectionSetPermutation() 378 379 Level: developer 380 381 TODO: How is this related to DMCreateLocalSection() 382 383 .seealso: DMPlexCreate(), PetscSectionCreate(), PetscSectionSetPermutation() 384 @*/ 385 PetscErrorCode DMPlexCreateSection(DM dm, DMLabel label[], const PetscInt numComp[],const PetscInt numDof[], PetscInt numBC, const PetscInt bcField[], const IS bcComps[], const IS bcPoints[], IS perm, PetscSection *section) 386 { 387 PetscSection aSec; 388 PetscErrorCode ierr; 389 390 PetscFunctionBegin; 391 ierr = DMPlexCreateSectionFields(dm, numComp, section);CHKERRQ(ierr); 392 ierr = DMPlexCreateSectionDof(dm, label, numDof, *section);CHKERRQ(ierr); 393 ierr = DMPlexCreateSectionBCDof(dm, numBC, bcField, bcComps, bcPoints, *section);CHKERRQ(ierr); 394 if (perm) {ierr = PetscSectionSetPermutation(*section, perm);CHKERRQ(ierr);} 395 ierr = PetscSectionSetFromOptions(*section);CHKERRQ(ierr); 396 ierr = PetscSectionSetUp(*section);CHKERRQ(ierr); 397 ierr = DMPlexGetAnchors(dm,&aSec,NULL);CHKERRQ(ierr); 398 if (numBC || aSec) { 399 ierr = DMPlexCreateSectionBCIndicesField(dm, numBC, bcField, bcComps, bcPoints, *section);CHKERRQ(ierr); 400 ierr = DMPlexCreateSectionBCIndices(dm, *section);CHKERRQ(ierr); 401 } 402 ierr = PetscSectionViewFromOptions(*section,NULL,"-section_view");CHKERRQ(ierr); 403 PetscFunctionReturn(0); 404 } 405 406 PetscErrorCode DMCreateLocalSection_Plex(DM dm) 407 { 408 PetscDS probBC; 409 PetscSection section; 410 DMLabel *labels; 411 IS *bcPoints, *bcComps; 412 PetscBool *isFE; 413 PetscInt *bcFields, *numComp, *numDof; 414 PetscInt depth, dim, numBd, numBC = 0, Nf, bd, bc = 0, f; 415 PetscInt cStart, cEnd, cEndInterior; 416 PetscErrorCode ierr; 417 418 PetscFunctionBegin; 419 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 420 /* FE and FV boundary conditions are handled slightly differently */ 421 ierr = PetscMalloc1(Nf, &isFE);CHKERRQ(ierr); 422 for (f = 0; f < Nf; ++f) { 423 PetscObject obj; 424 PetscClassId id; 425 426 ierr = DMGetField(dm, f, NULL, &obj);CHKERRQ(ierr); 427 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 428 if (id == PETSCFE_CLASSID) {isFE[f] = PETSC_TRUE;} 429 else if (id == PETSCFV_CLASSID) {isFE[f] = PETSC_FALSE;} 430 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f); 431 } 432 /* Allocate boundary point storage for FEM boundaries */ 433 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 434 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 435 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 436 ierr = DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);CHKERRQ(ierr); 437 ierr = DMGetDS(dm, &probBC);CHKERRQ(ierr); 438 ierr = PetscDSGetNumBoundary(probBC, &numBd);CHKERRQ(ierr); 439 if (!Nf && numBd) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "number of fields is zero and number of boundary conditions is nonzero (this should never happen)"); 440 for (bd = 0; bd < numBd; ++bd) { 441 PetscInt field; 442 DMBoundaryConditionType type; 443 const char *labelName; 444 DMLabel label; 445 446 ierr = PetscDSGetBoundary(probBC, bd, &type, NULL, &labelName, &field, NULL, NULL, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 447 ierr = DMGetLabel(dm,labelName,&label);CHKERRQ(ierr); 448 if (label && isFE[field] && (type & DM_BC_ESSENTIAL)) ++numBC; 449 } 450 /* Add ghost cell boundaries for FVM */ 451 for (f = 0; f < Nf; ++f) if (!isFE[f] && cEndInterior >= 0) ++numBC; 452 ierr = PetscCalloc3(numBC,&bcFields,numBC,&bcPoints,numBC,&bcComps);CHKERRQ(ierr); 453 /* Constrain ghost cells for FV */ 454 for (f = 0; f < Nf; ++f) { 455 PetscInt *newidx, c; 456 457 if (isFE[f] || cEndInterior < 0) continue; 458 ierr = PetscMalloc1(cEnd-cEndInterior,&newidx);CHKERRQ(ierr); 459 for (c = cEndInterior; c < cEnd; ++c) newidx[c-cEndInterior] = c; 460 bcFields[bc] = f; 461 ierr = ISCreateGeneral(PETSC_COMM_SELF, cEnd-cEndInterior, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);CHKERRQ(ierr); 462 } 463 /* Handle FEM Dirichlet boundaries */ 464 for (bd = 0; bd < numBd; ++bd) { 465 const char *bdLabel; 466 DMLabel label; 467 const PetscInt *comps; 468 const PetscInt *values; 469 PetscInt bd2, field, numComps, numValues; 470 DMBoundaryConditionType type; 471 PetscBool duplicate = PETSC_FALSE; 472 473 ierr = PetscDSGetBoundary(probBC, bd, &type, NULL, &bdLabel, &field, &numComps, &comps, NULL, &numValues, &values, NULL);CHKERRQ(ierr); 474 ierr = DMGetLabel(dm, bdLabel, &label);CHKERRQ(ierr); 475 if (!isFE[field] || !label) continue; 476 /* Only want to modify label once */ 477 for (bd2 = 0; bd2 < bd; ++bd2) { 478 const char *bdname; 479 ierr = PetscDSGetBoundary(probBC, bd2, NULL, NULL, &bdname, NULL, NULL, NULL, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 480 ierr = PetscStrcmp(bdname, bdLabel, &duplicate);CHKERRQ(ierr); 481 if (duplicate) break; 482 } 483 if (!duplicate && (isFE[field])) { 484 /* don't complete cells, which are just present to give orientation to the boundary */ 485 ierr = DMPlexLabelComplete(dm, label);CHKERRQ(ierr); 486 } 487 /* Filter out cells, if you actually want to constrain cells you need to do things by hand right now */ 488 if (type & DM_BC_ESSENTIAL) { 489 PetscInt *newidx; 490 PetscInt n, newn = 0, p, v; 491 492 bcFields[bc] = field; 493 if (numComps) {ierr = ISCreateGeneral(PETSC_COMM_SELF, numComps, comps, PETSC_COPY_VALUES, &bcComps[bc]);CHKERRQ(ierr);} 494 for (v = 0; v < numValues; ++v) { 495 IS tmp; 496 const PetscInt *idx; 497 498 ierr = DMGetStratumIS(dm, bdLabel, values[v], &tmp);CHKERRQ(ierr); 499 if (!tmp) continue; 500 ierr = ISGetLocalSize(tmp, &n);CHKERRQ(ierr); 501 ierr = ISGetIndices(tmp, &idx);CHKERRQ(ierr); 502 if (isFE[field]) { 503 for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) ++newn; 504 } else { 505 for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) ++newn; 506 } 507 ierr = ISRestoreIndices(tmp, &idx);CHKERRQ(ierr); 508 ierr = ISDestroy(&tmp);CHKERRQ(ierr); 509 } 510 ierr = PetscMalloc1(newn,&newidx);CHKERRQ(ierr); 511 newn = 0; 512 for (v = 0; v < numValues; ++v) { 513 IS tmp; 514 const PetscInt *idx; 515 516 ierr = DMGetStratumIS(dm, bdLabel, values[v], &tmp);CHKERRQ(ierr); 517 if (!tmp) continue; 518 ierr = ISGetLocalSize(tmp, &n);CHKERRQ(ierr); 519 ierr = ISGetIndices(tmp, &idx);CHKERRQ(ierr); 520 if (isFE[field]) { 521 for (p = 0; p < n; ++p) if ((idx[p] < cStart) || (idx[p] >= cEnd)) newidx[newn++] = idx[p]; 522 } else { 523 for (p = 0; p < n; ++p) if ((idx[p] >= cStart) || (idx[p] < cEnd)) newidx[newn++] = idx[p]; 524 } 525 ierr = ISRestoreIndices(tmp, &idx);CHKERRQ(ierr); 526 ierr = ISDestroy(&tmp);CHKERRQ(ierr); 527 } 528 ierr = ISCreateGeneral(PETSC_COMM_SELF, newn, newidx, PETSC_OWN_POINTER, &bcPoints[bc++]);CHKERRQ(ierr); 529 } 530 } 531 /* Handle discretization */ 532 ierr = PetscCalloc3(Nf,&labels,Nf,&numComp,Nf*(dim+1),&numDof);CHKERRQ(ierr); 533 for (f = 0; f < Nf; ++f) { 534 labels[f] = dm->fields[f].label; 535 if (isFE[f]) { 536 PetscFE fe = (PetscFE) dm->fields[f].disc; 537 const PetscInt *numFieldDof; 538 PetscInt fedim, d; 539 540 ierr = PetscFEGetNumComponents(fe, &numComp[f]);CHKERRQ(ierr); 541 ierr = PetscFEGetNumDof(fe, &numFieldDof);CHKERRQ(ierr); 542 ierr = PetscFEGetSpatialDimension(fe, &fedim);CHKERRQ(ierr); 543 for (d = 0; d < PetscMin(dim, fedim)+1; ++d) numDof[f*(dim+1)+d] = numFieldDof[d]; 544 } else { 545 PetscFV fv = (PetscFV) dm->fields[f].disc; 546 547 ierr = PetscFVGetNumComponents(fv, &numComp[f]);CHKERRQ(ierr); 548 numDof[f*(dim+1)+dim] = numComp[f]; 549 } 550 } 551 for (f = 0; f < Nf; ++f) { 552 PetscInt d; 553 for (d = 1; d < dim; ++d) { 554 if ((numDof[f*(dim+1)+d] > 0) && (depth < dim)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mesh must be interpolated when unknowns are specified on edges or faces."); 555 } 556 } 557 ierr = DMPlexCreateSection(dm, labels, numComp, numDof, numBC, bcFields, bcComps, bcPoints, NULL, §ion);CHKERRQ(ierr); 558 for (f = 0; f < Nf; ++f) { 559 PetscFE fe; 560 const char *name; 561 562 if (isFE[f]) { 563 ierr = DMGetField(dm, f, NULL, (PetscObject *) &fe);CHKERRQ(ierr); 564 ierr = PetscObjectGetName((PetscObject) fe, &name);CHKERRQ(ierr); 565 ierr = PetscSectionSetFieldName(section, f, name);CHKERRQ(ierr); 566 } 567 } 568 ierr = DMSetLocalSection(dm, section);CHKERRQ(ierr); 569 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 570 for (bc = 0; bc < numBC; ++bc) {ierr = ISDestroy(&bcPoints[bc]);CHKERRQ(ierr);ierr = ISDestroy(&bcComps[bc]);CHKERRQ(ierr);} 571 ierr = PetscFree3(bcFields,bcPoints,bcComps);CHKERRQ(ierr); 572 ierr = PetscFree3(labels,numComp,numDof);CHKERRQ(ierr); 573 ierr = PetscFree(isFE);CHKERRQ(ierr); 574 PetscFunctionReturn(0); 575 } 576