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