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 switch (dim) { 49 case 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 break; 54 case 3: 55 eqns[0*size+r] = PetscSqr(e[0]); 56 eqns[1*size+r] = 2.0*e[0]*e[1]; 57 eqns[2*size+r] = 2.0*e[0]*e[2]; 58 eqns[3*size+r] = PetscSqr(e[1]); 59 eqns[4*size+r] = 2.0*e[1]*e[2]; 60 eqns[5*size+r] = PetscSqr(e[2]); 61 break; 62 } 63 } 64 } 65 ierr = MatSetUnfactored(A);CHKERRQ(ierr); 66 ierr = DMPlexVecRestoreClosure(dm, coordSection, coordinates, c, NULL, &cellCoords);CHKERRQ(ierr); 67 ierr = MatLUFactor(A, NULL, NULL, NULL);CHKERRQ(ierr); 68 ierr = MatSolve(A, mb, mx);CHKERRQ(ierr); 69 ierr = VecGetArrayRead(mx, &sol);CHKERRQ(ierr); 70 ierr = DMPlexComputeCellGeometryFVM(dm, c, &vol, NULL, NULL);CHKERRQ(ierr); 71 ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr); 72 ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr); 73 for (cl = 0; cl < coneSize; ++cl) { 74 const PetscInt v = cone[cl] - vStart; 75 76 if (dim == 2) { 77 metric[v*4+0] += vol*coarseRatio*sol[0]; 78 metric[v*4+1] += vol*coarseRatio*sol[1]; 79 metric[v*4+2] += vol*coarseRatio*sol[1]; 80 metric[v*4+3] += vol*coarseRatio*sol[2]; 81 } else { 82 metric[v*9+0] += vol*coarseRatio*sol[0]; 83 metric[v*9+1] += vol*coarseRatio*sol[1]; 84 metric[v*9+3] += vol*coarseRatio*sol[1]; 85 metric[v*9+2] += vol*coarseRatio*sol[2]; 86 metric[v*9+6] += vol*coarseRatio*sol[2]; 87 metric[v*9+4] += vol*coarseRatio*sol[3]; 88 metric[v*9+5] += vol*coarseRatio*sol[4]; 89 metric[v*9+7] += vol*coarseRatio*sol[4]; 90 metric[v*9+8] += vol*coarseRatio*sol[5]; 91 } 92 } 93 ierr = VecRestoreArrayRead(mx, &sol);CHKERRQ(ierr); 94 } 95 for (v = 0; v < numVertices; ++v) { 96 const PetscInt *support; 97 PetscInt supportSize, s; 98 PetscReal vol, totVol = 0.0; 99 100 ierr = DMPlexGetSupport(udm, v+vStart, &support);CHKERRQ(ierr); 101 ierr = DMPlexGetSupportSize(udm, v+vStart, &supportSize);CHKERRQ(ierr); 102 for (s = 0; s < supportSize; ++s) {ierr = DMPlexComputeCellGeometryFVM(dm, support[s], &vol, NULL, NULL);CHKERRQ(ierr); totVol += vol;} 103 for (s = 0; s < PetscSqr(dim); ++s) metric[v*PetscSqr(dim)+s] /= totVol; 104 } 105 ierr = VecRestoreArray(metricVec, &metric);CHKERRQ(ierr); 106 ierr = VecDestroy(&mx);CHKERRQ(ierr); 107 ierr = VecDestroy(&mb);CHKERRQ(ierr); 108 ierr = MatDestroy(&A);CHKERRQ(ierr); 109 ierr = DMDestroy(&udm);CHKERRQ(ierr); 110 ierr = PetscFree(eqns);CHKERRQ(ierr); 111 112 bdLabelName[0] = '\0'; 113 ierr = PetscOptionsGetString(NULL, dm->hdr.prefix, "-dm_plex_coarsen_bd_label", bdLabelName, PETSC_MAX_PATH_LEN-1, NULL);CHKERRQ(ierr); 114 ierr = DMPlexRemesh_Internal(dm, metricVec, bdLabelName, mesh->remeshBd, &dm->coarseMesh);CHKERRQ(ierr); 115 ierr = VecDestroy(&metricVec);CHKERRQ(ierr); 116 } 117 ierr = PetscObjectReference((PetscObject) dm->coarseMesh);CHKERRQ(ierr); 118 *dmCoarsened = dm->coarseMesh; 119 PetscFunctionReturn(0); 120 } 121 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