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, const char bdyLabelName[], DM *dmCoarsened) 11 { 12 #ifdef PETSC_HAVE_PRAGMATIC 13 DM udm, coordDM; 14 DMLabel bd, bdtags; 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 if ( bdyLabelName ) { 79 ierr = DMGetLabel(dm, bdyLabelName, &bdtags);CHKERRQ(ierr); 80 } 81 // TODO fix if bdyLabelName = "boundary" 82 // TODO a way to use only bdyLabelName would be to loop over its strata, if I can't get all the strata at once 83 ierr = DMLabelCreate("boundary", &bd);CHKERRQ(ierr); 84 ierr = DMPlexMarkBoundaryFaces(dm, bd);CHKERRQ(ierr); 85 ierr = DMLabelGetStratumIS(bd, 1, &bdIS);CHKERRQ(ierr); 86 ierr = DMLabelGetStratumSize(bd, 1, &numBdFaces);CHKERRQ(ierr); 87 ierr = ISGetIndices(bdIS, &faces);CHKERRQ(ierr); 88 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 89 PetscInt *closure = NULL; 90 PetscInt closureSize, cl; 91 92 ierr = DMPlexGetTransitiveClosure(dm, faces[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 93 for (cl = 0; cl < closureSize*2; cl += 2) { 94 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize; 95 } 96 ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 97 } 98 ierr = PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);CHKERRQ(ierr); 99 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 100 PetscInt *closure = NULL; 101 PetscInt closureSize, cl; 102 103 ierr = DMPlexGetTransitiveClosure(dm, faces[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 104 for (cl = 0; cl < closureSize*2; cl += 2) { 105 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart; 106 } 107 if ( bdyLabelName ) 108 DMLabelGetValue(bdtags, faces[f], &bdFaceIds[f]); 109 else 110 bdFaceIds[f] = 1; 111 ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 112 } 113 ierr = ISDestroy(&bdIS);CHKERRQ(ierr); 114 ierr = DMLabelDestroy(&bd);CHKERRQ(ierr); 115 pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds); 116 pragmatic_set_metric(met); 117 pragmatic_adapt(); 118 /* Read out mesh */ 119 pragmatic_get_info(&numCoarseVertices, &numCoarseCells); 120 ierr = PetscMalloc1(numCoarseVertices*dim, &coarseCoords);CHKERRQ(ierr); 121 switch (dim) { 122 case 2: 123 ierr = PetscMalloc2(numCoarseVertices, &x, numCoarseVertices, &y);CHKERRQ(ierr); 124 pragmatic_get_coords_2d(x, y); 125 numCorners = 3; 126 for (v = 0; v < numCoarseVertices; ++v) {coarseCoords[v*2+0] = x[v]; coarseCoords[v*2+1] = y[v];} 127 break; 128 case 3: 129 ierr = PetscMalloc3(numCoarseVertices, &x, numCoarseVertices, &y, numCoarseVertices, &z);CHKERRQ(ierr); 130 pragmatic_get_coords_3d(x, y, z); 131 numCorners = 4; 132 for (v = 0; v < numCoarseVertices; ++v) {coarseCoords[v*3+0] = x[v]; coarseCoords[v*3+1] = y[v]; coarseCoords[v*3+2] = z[v];} 133 break; 134 default: 135 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic coarsening defined for dimension %d", dim); 136 } 137 ierr = PetscMalloc1(numCoarseCells*(dim+1), &ccells);CHKERRQ(ierr); // only for simplicial meshes 138 pragmatic_get_elements(ccells); 139 /* TODO Read out markers for boundary */ 140 ierr = DMPlexCreateFromCellList(PetscObjectComm((PetscObject) dm), dim, numCoarseCells, numCoarseVertices, numCorners, PETSC_TRUE, ccells, dim, coarseCoords, &dm->coarseMesh);CHKERRQ(ierr); 141 pragmatic_finalize(); 142 ierr = PetscFree4(x, y, z, cells);CHKERRQ(ierr); 143 ierr = PetscFree2(bdFaces, bdFaceIds);CHKERRQ(ierr); 144 ierr = PetscFree(coarseCoords);CHKERRQ(ierr); 145 } 146 #endif 147 ierr = PetscObjectReference((PetscObject) dm->coarseMesh);CHKERRQ(ierr); 148 *dmCoarsened = dm->coarseMesh; 149 PetscFunctionReturn(0); 150 } 151