1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 #undef __FUNCT__ 4 #define __FUNCT__ "DMPlexCreateClosureIndex" 5 /*@ 6 DMPlexCreateClosureIndex - Calculate an index for the given PetscSection for the closure operation on the DM 7 8 Not collective 9 10 Input Parameters: 11 + dm - The DM 12 - section - The section describing the layout in v, or NULL to use the default section 13 14 Note: 15 This should greatly improve the performance of the closure operations, at the cost of additional memory. 16 17 Level: intermediate 18 19 .seealso DMPlexVecGetClosure(), DMPlexVecRestoreClosure(), DMPlexVecSetClosure(), DMPlexMatSetClosure() 20 @*/ 21 PetscErrorCode DMPlexCreateClosureIndex(DM dm, PetscSection section) 22 { 23 PetscSection closureSection; 24 IS closureIS; 25 PetscInt *clPoints; 26 PetscInt pStart, pEnd, sStart, sEnd, point, clSize; 27 PetscErrorCode ierr; 28 29 PetscFunctionBegin; 30 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 31 if (!section) {ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr);} 32 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 2); 33 ierr = PetscSectionGetChart(section, &sStart, &sEnd);CHKERRQ(ierr); 34 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 35 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) section), &closureSection);CHKERRQ(ierr); 36 ierr = PetscSectionSetChart(closureSection, pStart, pEnd);CHKERRQ(ierr); 37 for (point = pStart; point < pEnd; ++point) { 38 PetscInt *points = NULL, numPoints, p, dof, cldof = 0; 39 40 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 41 for (p = 0; p < numPoints*2; p += 2) { 42 if ((points[p] >= sStart) && (points[p] < sEnd)) { 43 ierr = PetscSectionGetDof(section, points[p], &dof);CHKERRQ(ierr); 44 if (dof) cldof += 2; 45 } 46 } 47 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 48 ierr = PetscSectionSetDof(closureSection, point, cldof);CHKERRQ(ierr); 49 } 50 ierr = PetscSectionSetUp(closureSection);CHKERRQ(ierr); 51 ierr = PetscSectionGetStorageSize(closureSection, &clSize);CHKERRQ(ierr); 52 ierr = PetscMalloc1(clSize, &clPoints);CHKERRQ(ierr); 53 for (point = pStart; point < pEnd; ++point) { 54 PetscInt *points = NULL, numPoints, p, q, dof, cldof, cloff; 55 56 ierr = PetscSectionGetDof(closureSection, point, &cldof);CHKERRQ(ierr); 57 ierr = PetscSectionGetOffset(closureSection, point, &cloff);CHKERRQ(ierr); 58 ierr = DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 59 for (p = 0, q = 0; p < numPoints*2; p += 2) { 60 if ((points[p] >= sStart) && (points[p] < sEnd)) { 61 ierr = PetscSectionGetDof(section, points[p], &dof);CHKERRQ(ierr); 62 if (dof) { 63 clPoints[cloff+q*2] = points[p]; 64 clPoints[cloff+q*2+1] = points[p+1]; 65 ++q; 66 } 67 } 68 } 69 ierr = DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 70 if (q*2 != cldof) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_PLIB, "Invalid size for closure %d should be %d", q*2, cldof); 71 } 72 ierr = ISCreateGeneral(PETSC_COMM_SELF, clSize, clPoints, PETSC_OWN_POINTER, &closureIS);CHKERRQ(ierr); 73 ierr = PetscSectionSetClosureIndex(section, (PetscObject) dm, closureSection, closureIS); 74 PetscFunctionReturn(0); 75 } 76