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