1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #ifdef PETSC_HAVE_PRAGMATIC 3 #include <pragmatic/cpragmatic.h> 4 #endif 5 6 #undef __FUNCT__ 7 #define __FUNCT__ "DMCoarsen_Plex" 8 PetscErrorCode DMCoarsen_Plex(DM dm, MPI_Comm comm, DM *dmCoarsened) 9 { 10 DM_Plex *mesh = (DM_Plex*) dm->data; 11 DM udm, coordDM; 12 DMLabel bd; 13 Vec coordinates; 14 PetscSection coordSection; 15 const PetscScalar *coords; 16 double *coarseCoords; 17 IS bdIS; 18 PetscReal *x, *y, *z, *metric; 19 const PetscInt *faces; 20 PetscInt *cells, *bdFaces, *bdFaceIds; 21 PetscInt dim, numCorners, cStart, cEnd, numCells, numCoarseCells, c, vStart, vEnd, numVertices, numCoarseVertices, v, numBdFaces, numCoarseBdFaces, f, maxConeSize, bdSize, coff; 22 PetscErrorCode ierr; 23 24 PetscFunctionBegin; 25 #ifdef PETSC_HAVE_PRAGMATIC 26 if (!mesh->coarseMesh) { 27 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 28 ierr = DMGetCoordinateDM(dm, &coordDM);CHKERRQ(ierr); 29 ierr = DMGetDefaultSection(coordDM, &coordSection);CHKERRQ(ierr); 30 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 31 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 32 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 33 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr); 34 ierr = DMPlexGetMaxSizes(udm, &maxConeSize, NULL);CHKERRQ(ierr); 35 numCells = cEnd - cStart; 36 numVertices = vEnd - vStart; 37 ierr = PetscCalloc5(numVertices, &x, numVertices, &y, numVertices, &z, numVertices, &metric, numCells*maxConeSize, &cells);CHKERRQ(ierr); 38 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 39 for (v = 0; v < numVertices; ++v) { 40 PetscInt off; 41 42 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 43 x[v] = coords[off+0]; 44 y[v] = coords[off+1]; 45 if (dim > 2) z[v] = coords[off+2]; 46 } 47 ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr); 48 for (c = 0, coff = 0; c < numCells; ++c) { 49 const PetscInt *cone; 50 PetscInt coneSize, cl; 51 52 ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr); 53 ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr); 54 for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl]; 55 } 56 switch (dim) { 57 case 2: 58 pragmatic_2d_init(&numVertices, &numCells, cells, x, y); 59 break; 60 case 3: 61 pragmatic_3d_init(&numVertices, &numCells, cells, x, y, z); 62 break; 63 default: 64 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic coarsening defined for dimension %d", dim); 65 } 66 ierr = DMDestroy(&udm);CHKERRQ(ierr); 67 /* Create boundary mesh */ 68 ierr = DMLabelCreate("boundary", &bd);CHKERRQ(ierr); 69 ierr = DMPlexMarkBoundaryFaces(dm, bd);CHKERRQ(ierr); 70 ierr = DMLabelGetStratumIS(bd, 1, &bdIS);CHKERRQ(ierr); 71 ierr = DMLabelGetStratumSize(bd, 1, &numBdFaces);CHKERRQ(ierr); 72 ierr = ISGetIndices(bdIS, &faces);CHKERRQ(ierr); 73 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 74 PetscInt *closure = NULL; 75 PetscInt closureSize, cl; 76 77 ierr = DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 78 for (cl = 0; cl < closureSize*2; cl += 2) { 79 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize; 80 } 81 ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 82 } 83 ierr = PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);CHKERRQ(ierr); 84 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 85 PetscInt *closure = NULL; 86 PetscInt closureSize, cl; 87 88 ierr = DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 89 for (cl = 0; cl < closureSize*2; cl += 2) { 90 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl]; 91 } 92 /* TODO Fix */ 93 bdFaceIds[f] = 1; 94 ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 95 } 96 ierr = ISDestroy(&bdIS);CHKERRQ(ierr); 97 pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds); 98 /* Create metric */ 99 if (dim == 2) {for (v = 0; v < numVertices; ++v) metric[v*4+0] = metric[v*4+3] = 1.0;} 100 else {for (v = 0; v < numVertices; ++v) metric[v*9+0] = metric[v*9+4] = metric[v*9+8] = 1.0;} 101 pragmatic_set_metric(metric); 102 pragmatic_adapt(); 103 /* Read out mesh */ 104 pragmatic_get_info(&numCoarseVertices, &numCoarseCells); 105 switch (dim) { 106 case 2: 107 pragmatic_get_coords_2d(x, y); 108 numCorners = 3; 109 for (v = 0; v < numCoarseVertices; ++v) {coarseCoords[v*2+0] = x[v]; coarseCoords[v*2+1] = y[v];} 110 break; 111 case 3: 112 pragmatic_get_coords_3d(x, y, z); 113 numCorners = 4; 114 for (v = 0; v < numCoarseVertices; ++v) {coarseCoords[v*3+0] = x[v]; coarseCoords[v*3+1] = y[v]; coarseCoords[v*3+2] = z[v];} 115 break; 116 default: 117 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic coarsening defined for dimension %d", dim); 118 } 119 pragmatic_get_elements(cells); 120 /* TODO Read out markers for boundary */ 121 ierr = PetscMalloc1(numCoarseVertices*dim, &coarseCoords);CHKERRQ(ierr); 122 ierr = DMPlexCreateFromCellList(PetscObjectComm((PetscObject) dm), dim, numCoarseCells, numCoarseVertices, numCorners, PETSC_TRUE, cells, dim, coarseCoords, &mesh->coarseMesh);CHKERRQ(ierr); 123 pragmatic_finalize(); 124 ierr = PetscFree4(x, y, z, metric);CHKERRQ(ierr); 125 ierr = PetscFree2(bdFaces, bdFaceIds);CHKERRQ(ierr); 126 ierr = PetscFree(coarseCoords);CHKERRQ(ierr); 127 } 128 #endif 129 ierr = PetscObjectReference((PetscObject) mesh->coarseMesh);CHKERRQ(ierr); 130 *dmCoarsened = mesh->coarseMesh; 131 PetscFunctionReturn(0); 132 } 133