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