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