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__ "DMPlexAdapt" 10 /*@ 11 DMPlexAdapt - Generates a mesh adapted to the specified metric field using the pragmatic library. 12 13 Input Parameters: 14 + dm - The DM object 15 . metric - The metric to which the mesh is adapted, defined vertex-wise. 16 - bdyLabelName - Label name for boundary tags. These will be preserved in the output mesh. bdyLabelName should be different from "boundary". 17 18 Output Parameter: 19 . dmAdapted - Pointer to the DM object containing the adapted mesh 20 21 Level: advanced 22 23 .seealso: DMCoarsen(), DMRefine() 24 @*/ 25 PetscErrorCode DMPlexAdapt(DM dm, Vec metric, const char bdyLabelName[], DM *dmAdapted) 26 { 27 #ifdef PETSC_HAVE_PRAGMATIC 28 DM udm, coordDM; 29 DMLabel bd, bdTags, bdTagsAdp; 30 Vec coordinates; 31 PetscSection coordSection; 32 IS bdIS; 33 double *coordsAdp; 34 const PetscScalar *coords; 35 PetscReal *x, *y, *z, *xAdp, *yAdp, *zAdp; 36 const PetscScalar *metricArray; 37 PetscReal *met; 38 const PetscInt *faces; 39 PetscInt *cells, *cellsAdp, *bdFaces, *bdFaceIds; 40 PetscInt *boundaryTags; 41 PetscInt dim, numCornersAdp, cStart, cEnd, numCells, numCellsAdp, c; 42 PetscInt vStart, vEnd, numVertices, numVerticesAdp, v; 43 PetscInt numBdFaces, f, maxConeSize, bdSize, coff; 44 PetscInt *cell, iVer, cellOff, i, j, idx, facet; 45 PetscInt perm[4], isInFacetClosure[4]; // 4 = max number of facets for an element in 2D & 3D. Only for simplicial meshes 46 PetscBool flag, hasBdyFacet; 47 #endif 48 PetscErrorCode ierr; 49 MPI_Comm comm; 50 51 PetscFunctionBegin; 52 ierr = PetscObjectGetComm((PetscObject)dm, &comm);CHKERRQ(ierr); 53 #ifdef PETSC_HAVE_PRAGMATIC 54 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 55 ierr = DMGetCoordinateDM(dm, &coordDM);CHKERRQ(ierr); 56 ierr = DMGetDefaultSection(coordDM, &coordSection);CHKERRQ(ierr); 57 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 58 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 59 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 60 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr); 61 ierr = DMPlexGetMaxSizes(udm, &maxConeSize, NULL);CHKERRQ(ierr); 62 numCells = cEnd - cStart; 63 numVertices = vEnd - vStart; 64 ierr = PetscCalloc5(numVertices, &x, numVertices, &y, numVertices, &z, numCells*maxConeSize, &cells, dim*dim*numVertices, &met);CHKERRQ(ierr); 65 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 66 ierr = VecGetArrayRead(metric, &metricArray);CHKERRQ(ierr); 67 for (v = vStart; v < vEnd; ++v) { 68 PetscInt off; 69 70 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 71 x[v-vStart] = coords[off+0]; 72 y[v-vStart] = coords[off+1]; 73 if (dim > 2) z[v-vStart] = coords[off+2]; 74 ierr = PetscMemcpy(&met[dim*dim*(v-vStart)], &metricArray[dim*off], dim*dim*sizeof(PetscScalar));CHKERRQ(ierr); 75 } 76 ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr); 77 ierr = VecRestoreArrayRead(metric, &metricArray);CHKERRQ(ierr); 78 for (c = 0, coff = 0; c < numCells; ++c) { 79 const PetscInt *cone; 80 PetscInt coneSize, cl; 81 82 ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr); 83 ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr); 84 for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart; 85 } 86 switch (dim) { 87 case 2: 88 pragmatic_2d_init(&numVertices, &numCells, cells, x, y); 89 break; 90 case 3: 91 pragmatic_3d_init(&numVertices, &numCells, cells, x, y, z); 92 break; 93 default: 94 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %d", dim); 95 } 96 /* Create boundary mesh */ 97 if (bdyLabelName) { 98 ierr = PetscStrcmp(bdyLabelName, "boundary", &flag);CHKERRQ(ierr); 99 if (flag) { 100 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdyLabelName); 101 } 102 ierr = DMGetLabel(dm, bdyLabelName, &bdTags);CHKERRQ(ierr); 103 } 104 /* TODO To avoid marking bdy facets again if bdyLabelName exists, could/should we loop over bdyLabelName stratas ? */ 105 ierr = DMLabelCreate("boundary", &bd);CHKERRQ(ierr); 106 ierr = DMPlexMarkBoundaryFaces(dm, bd);CHKERRQ(ierr); 107 ierr = DMLabelGetStratumIS(bd, 1, &bdIS);CHKERRQ(ierr); 108 ierr = DMLabelGetStratumSize(bd, 1, &numBdFaces);CHKERRQ(ierr); 109 ierr = ISGetIndices(bdIS, &faces);CHKERRQ(ierr); 110 /* TODO why not assume that we are considering simplicial meshes, in which case bdSize = dim*numBdFaces ? */ 111 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 112 PetscInt *closure = NULL; 113 PetscInt closureSize, cl; 114 115 ierr = DMPlexGetTransitiveClosure(dm, faces[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 116 for (cl = 0; cl < closureSize*2; cl += 2) { 117 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize; 118 } 119 ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 120 } 121 ierr = PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);CHKERRQ(ierr); 122 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 123 PetscInt *closure = NULL; 124 PetscInt closureSize, cl; 125 126 ierr = DMPlexGetTransitiveClosure(dm, faces[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 127 for (cl = 0; cl < closureSize*2; cl += 2) { 128 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart; 129 } 130 if (bdyLabelName) { 131 ierr = DMLabelGetValue(bdTags, faces[f], &bdFaceIds[f]);CHKERRQ(ierr); 132 } else { 133 bdFaceIds[f] = 1; 134 } 135 ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 136 } 137 ierr = ISDestroy(&bdIS);CHKERRQ(ierr); 138 ierr = DMLabelDestroy(&bd);CHKERRQ(ierr); 139 pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds); 140 pragmatic_set_metric(met); 141 pragmatic_adapt(); 142 ierr = PetscFree5(x, y, z, cells, met);CHKERRQ(ierr); 143 ierr = PetscFree2(bdFaces, bdFaceIds);CHKERRQ(ierr); 144 ierr = PetscFree(coords);CHKERRQ(ierr); 145 /* Read out mesh */ 146 pragmatic_get_info(&numVerticesAdp, &numCellsAdp); 147 ierr = PetscMalloc1(numVerticesAdp*dim, &coordsAdp);CHKERRQ(ierr); 148 switch (dim) { 149 case 2: 150 ierr = PetscMalloc2(numVerticesAdp, &xAdp, numVerticesAdp, &yAdp);CHKERRQ(ierr); 151 zAdp = NULL; 152 pragmatic_get_coords_2d(xAdp, yAdp); 153 numCornersAdp = 3; 154 for (v = 0; v < numVerticesAdp; ++v) {coordsAdp[v*2+0] = xAdp[v]; coordsAdp[v*2+1] = yAdp[v];} 155 break; 156 case 3: 157 ierr = PetscMalloc3(numVerticesAdp, &xAdp, numVerticesAdp, &yAdp, numVerticesAdp, &zAdp);CHKERRQ(ierr); 158 pragmatic_get_coords_3d(xAdp, yAdp, zAdp); 159 numCornersAdp = 4; 160 for (v = 0; v < numVerticesAdp; ++v) {coordsAdp[v*3+0] = xAdp[v]; coordsAdp[v*3+1] = yAdp[v]; coordsAdp[v*3+2] = zAdp[v];} 161 break; 162 default: 163 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %d", dim); 164 } 165 ierr = PetscMalloc1(numCellsAdp*(dim+1), &cellsAdp);CHKERRQ(ierr); // only for simplicial meshes 166 pragmatic_get_elements(cellsAdp); 167 ierr = DMPlexCreate(comm, dmAdapted);CHKERRQ(ierr); 168 ierr = DMPlexCreateFromCellList(PetscObjectComm((PetscObject) dm), dim, numCellsAdp, numVerticesAdp, numCornersAdp, PETSC_TRUE, cellsAdp, dim, coordsAdp, dmAdapted);CHKERRQ(ierr); 169 /* Read out boundary tags */ 170 pragmatic_get_boundaryTags(&boundaryTags); 171 if (!bdyLabelName) bdyLabelName = "boundary"; 172 ierr = DMCreateLabel(*dmAdapted, bdyLabelName);CHKERRQ(ierr); 173 ierr = DMGetLabel(*dmAdapted, bdyLabelName, &bdTagsAdp);CHKERRQ(ierr); 174 ierr = DMPlexGetHeightStratum(*dmAdapted, 0, &cStart, &cEnd);CHKERRQ(ierr); 175 ierr = DMPlexGetDepthStratum(*dmAdapted, 0, &vStart, &vEnd);CHKERRQ(ierr); 176 for (c = cStart; c < cEnd; ++c) { 177 PetscInt *cellClosure = NULL; 178 PetscInt cellClosureSize, cl; 179 PetscInt *facetClosure = NULL; 180 PetscInt facetClosureSize, cl2; 181 const PetscInt *facetList; 182 PetscInt facetListSize, f; 183 184 cellOff = (c-cStart)*(dim+1); // gives the offset corresponding to the cell in pragmatic boundary arrays. Only for simplicial meshes 185 hasBdyFacet = 0; 186 for (i = 0; i < dim+1; ++i) { // only for simplicial meshes 187 if (boundaryTags[cellOff+i]) { 188 hasBdyFacet = 1; 189 break; 190 } 191 } 192 if (!hasBdyFacet) continue; // The cell has no boundary edge/facet => no boundary tagging 193 194 cell = &cellsAdp[cellOff]; // pointing to the current cell in the cellsAdp table 195 ierr = DMPlexGetTransitiveClosure(*dmAdapted, c, PETSC_TRUE, &cellClosureSize, &cellClosure);CHKERRQ(ierr); 196 /* first, encode the permutation of the vertices indices betwenn the cell closure and pragmatic ordering */ 197 j = 0; 198 for (cl = 0; cl < cellClosureSize*2; cl += 2) { 199 if ((cellClosure[cl] < vStart) || (cellClosure[cl] >= vEnd)) continue; 200 iVer = cellClosure[cl] - vStart; 201 for (i = 0; i < dim+1; ++i) { // only for simplicial meshes 202 if (iVer == cell[i]) { 203 perm[j] = i; // the cl-th element of the closure is the i-th vertex of the cell 204 break; 205 } 206 } 207 j++; 208 } 209 ierr = DMPlexGetCone(*dmAdapted, c, &facetList);CHKERRQ(ierr); // list of edges/facets of the cell 210 ierr = DMPlexGetConeSize(*dmAdapted, c, &facetListSize);CHKERRQ(ierr); 211 /* then, for each edge/facet of the cell, find the opposite vertex (ie the one in the closure of the cell and not in the closure of the facet/edge) */ 212 for (f = 0; f < facetListSize; ++f){ 213 facet = facetList[f]; 214 ierr = DMPlexGetTransitiveClosure(*dmAdapted, facet, PETSC_TRUE, &facetClosureSize, &facetClosure);CHKERRQ(ierr); 215 /* loop over the vertices of the cell closure, and check if each vertex belongs to the edge/facet closure */ 216 /* TODO should we use bitmaps to mark vertices instead of a static array ? */ 217 PetscMemzero(isInFacetClosure, sizeof(isInFacetClosure)); 218 j = 0; 219 for (cl = 0; cl < cellClosureSize*2; cl += 2) { 220 if ((cellClosure[cl] < vStart) || (cellClosure[cl] >= vEnd)) continue; 221 for (cl2 = 0; cl2 < facetClosureSize*2; cl2 += 2){ 222 if ((facetClosure[cl2] < vStart) || (facetClosure[cl2] >= vEnd)) continue; 223 if (cellClosure[cl] == facetClosure[cl2]) { 224 isInFacetClosure[j] = 1; 225 } 226 } 227 j++; 228 } 229 /* the vertex that was not marked is the vertex opposite to the edge/facet, ie the one giving the edge/facet boundary tag in pragmatic */ 230 j = 0; 231 for (cl = 0; cl < cellClosureSize*2; cl += 2) { 232 if ((cellClosure[cl] < vStart) || (cellClosure[cl] >= vEnd)) continue; 233 if (!isInFacetClosure[j]) { 234 idx = cellOff + perm[j]; 235 if (boundaryTags[idx]) { 236 ierr = DMLabelSetValue(bdTagsAdp, facet, boundaryTags[idx]);CHKERRQ(ierr); 237 } 238 break; 239 } 240 j++; 241 } 242 ierr = DMPlexRestoreTransitiveClosure(*dmAdapted, facet, PETSC_TRUE, &facetClosureSize, &facetClosure);CHKERRQ(ierr); 243 } 244 ierr = DMPlexRestoreTransitiveClosure(*dmAdapted, c, PETSC_TRUE, &cellClosureSize, &cellClosure);CHKERRQ(ierr); 245 } 246 pragmatic_finalize(); 247 ierr = PetscFree5(xAdp, yAdp, zAdp, cellsAdp, coordsAdp);CHKERRQ(ierr); 248 #endif 249 PetscFunctionReturn(0); 250 } 251