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