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