1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 #undef __FUNCT__ 4 #define __FUNCT__ "DMCoarsen_Plex" 5 PetscErrorCode DMCoarsen_Plex(DM dm, MPI_Comm comm, DM *dmCoarsened) 6 { 7 const PetscReal coarseRatio = PetscSqr(0.5); 8 DM udm, coordDM; 9 Mat A; 10 Vec metricVec, coordinates, mb, mx; 11 PetscSection coordSection; 12 PetscScalar *metric; 13 PetscScalar *eqns; 14 PetscInt dim, cStart, cEnd, c, vStart, vEnd, numVertices, v, size; 15 char bdLabelName[PETSC_MAX_PATH_LEN]; 16 PetscErrorCode ierr; 17 18 PetscFunctionBegin; 19 if (!dm->coarseMesh) { 20 /* Create metric */ 21 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr); 22 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 23 ierr = DMGetCoordinateDM(dm, &coordDM);CHKERRQ(ierr); 24 ierr = DMGetDefaultSection(coordDM, &coordSection);CHKERRQ(ierr); 25 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 26 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 27 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 28 numVertices = vEnd - vStart; 29 ierr = VecCreateSeq(PETSC_COMM_SELF, numVertices*PetscSqr(dim), &metricVec);CHKERRQ(ierr); 30 ierr = VecGetArray(metricVec, &metric);CHKERRQ(ierr); 31 size = (dim*(dim+1))/2; 32 ierr = PetscMalloc1(PetscSqr(size), &eqns);CHKERRQ(ierr); 33 ierr = MatCreateSeqDense(PETSC_COMM_SELF, size, size, eqns, &A);CHKERRQ(ierr); 34 ierr = MatCreateVecs(A, &mx, &mb);CHKERRQ(ierr); 35 ierr = VecSet(mb, 1.0);CHKERRQ(ierr); 36 for (c = cStart; c < cEnd; ++c) { 37 const PetscScalar *sol; 38 PetscScalar *cellCoords = NULL; 39 PetscReal e[3], vol; 40 const PetscInt *cone; 41 PetscInt coneSize, cl, i, j, d, r; 42 43 ierr = DMPlexVecGetClosure(dm, coordSection, coordinates, c, NULL, &cellCoords);CHKERRQ(ierr); 44 /* Only works for simplices */ 45 for (i = 0, r = 0; i < dim+1; ++i) { 46 for (j = 0; j < i; ++j, ++r) { 47 for (d = 0; d < dim; ++d) e[d] = PetscRealPart(cellCoords[i*dim+d] - cellCoords[j*dim+d]); 48 /* FORTRAN ORDERING */ 49 if (dim == 2) { 50 eqns[0*size+r] = PetscSqr(e[0]); 51 eqns[1*size+r] = 2.0*e[0]*e[1]; 52 eqns[2*size+r] = PetscSqr(e[1]); 53 } else { 54 eqns[0*size+r] = PetscSqr(e[0]); 55 eqns[1*size+r] = 2.0*e[0]*e[1]; 56 eqns[2*size+r] = 2.0*e[0]*e[2]; 57 eqns[3*size+r] = PetscSqr(e[1]); 58 eqns[4*size+r] = 2.0*e[1]*e[2]; 59 eqns[5*size+r] = PetscSqr(e[2]); 60 } 61 } 62 } 63 ierr = MatSetUnfactored(A);CHKERRQ(ierr); 64 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, NULL, &cellCoords);CHKERRQ(ierr); 65 ierr = MatLUFactor(A, NULL, NULL, NULL);CHKERRQ(ierr); 66 ierr = MatSolve(A, mb, mx);CHKERRQ(ierr); 67 ierr = VecGetArrayRead(mx, &sol);CHKERRQ(ierr); 68 ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL);CHKERRQ(ierr); 69 ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr); 70 ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr); 71 for (cl = 0; cl < coneSize; ++cl) { 72 const PetscInt v = cone[cl] - vStart; 73 74 if (dim == 2) { 75 metric[v*4+0] += vol*coarseRatio*sol[0]; 76 metric[v*4+1] += vol*coarseRatio*sol[1]; 77 metric[v*4+2] += vol*coarseRatio*sol[1]; 78 metric[v*4+3] += vol*coarseRatio*sol[2]; 79 } else { 80 metric[v*9+0] += vol*coarseRatio*sol[0]; 81 metric[v*9+1] += vol*coarseRatio*sol[1]; 82 metric[v*9+3] += vol*coarseRatio*sol[1]; 83 metric[v*9+2] += vol*coarseRatio*sol[2]; 84 metric[v*9+6] += vol*coarseRatio*sol[2]; 85 metric[v*9+4] += vol*coarseRatio*sol[3]; 86 metric[v*9+5] += vol*coarseRatio*sol[4]; 87 metric[v*9+7] += vol*coarseRatio*sol[4]; 88 metric[v*9+8] += vol*coarseRatio*sol[5]; 89 } 90 } 91 ierr = VecRestoreArrayRead(mx, &sol);CHKERRQ(ierr); 92 } 93 for (v = 0; v < numVertices; ++v) { 94 const PetscInt *support; 95 PetscInt supportSize, s; 96 PetscReal vol, totVol = 0.0; 97 98 ierr = DMPlexGetSupport(udm, v+vStart, &support);CHKERRQ(ierr); 99 ierr = DMPlexGetSupportSize(udm, v+vStart, &supportSize);CHKERRQ(ierr); 100 for (s = 0; s < supportSize; ++s) {ierr = DMPlexComputeCellGeometryFVM(dm, support[s], &vol, NULL, NULL);CHKERRQ(ierr); totVol += vol;} 101 for (s = 0; s < PetscSqr(dim); ++s) metric[v*PetscSqr(dim)+s] /= totVol; 102 } 103 ierr = VecRestoreArray(metricVec, &metric);CHKERRQ(ierr); 104 ierr = VecDestroy(&mx);CHKERRQ(ierr); 105 ierr = VecDestroy(&mb);CHKERRQ(ierr); 106 ierr = MatDestroy(&A);CHKERRQ(ierr); 107 ierr = DMDestroy(&udm);CHKERRQ(ierr); 108 ierr = PetscFree(eqns);CHKERRQ(ierr); 109 110 bdLabelName[0] = '\0'; 111 ierr = PetscOptionsGetString(NULL, dm->hdr.prefix, "-plex_coarsen_bd_label", bdLabelName, PETSC_MAX_PATH_LEN-1, NULL);CHKERRQ(ierr); 112 ierr = DMPlexRemesh_Internal(dm, metricVec, bdLabelName, &dm->coarseMesh);CHKERRQ(ierr); 113 ierr = VecDestroy(&metricVec);CHKERRQ(ierr); 114 } 115 ierr = PetscObjectReference((PetscObject) dm->coarseMesh);CHKERRQ(ierr); 116 *dmCoarsened = dm->coarseMesh; 117 PetscFunctionReturn(0); 118 } 119 120 #undef __FUNCT__ 121 #define __FUNCT__ "DMCoarsenHierarchy_Plex" 122 PetscErrorCode DMCoarsenHierarchy_Plex(DM dm, PetscInt nlevels, DM dmCoarsened[]) 123 { 124 DM rdm = dm; 125 PetscInt c; 126 PetscErrorCode ierr; 127 128 PetscFunctionBegin; 129 for (c = nlevels-1; c >= 0; --c) { 130 ierr = DMCoarsen(rdm, PetscObjectComm((PetscObject) dm), &dmCoarsened[c]);CHKERRQ(ierr); 131 ierr = DMCopyBoundary(rdm, dmCoarsened[c]);CHKERRQ(ierr); 132 ierr = DMSetCoarseDM(rdm, dmCoarsened[c]);CHKERRQ(ierr); 133 rdm = dmCoarsened[c]; 134 } 135 PetscFunctionReturn(0); 136 } 137