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