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 "" (empty string) if there is no such label, and 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, bdyLabelExt, 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 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_NULL, "Argument bdyLabelName cannot be NULL", dim); 99 } else { 100 ierr = PetscStrcmp(bdyLabelName, "", &flag);CHKERRQ(ierr); 101 if (flag) bdyLabelExt = PETSC_FALSE; 102 else bdyLabelExt = PETSC_TRUE; 103 } 104 if (bdyLabelExt) { 105 ierr = PetscStrcmp(bdyLabelName, "boundary", &flag);CHKERRQ(ierr); 106 if (flag) { 107 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdyLabelName); 108 } 109 ierr = DMGetLabel(dm, bdyLabelName, &bdTags);CHKERRQ(ierr); 110 if (!bdTags) { 111 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Label %s does not exist in DM", bdyLabelName); 112 } 113 } 114 /* TODO To avoid marking bdy facets again if bdyLabelName exists, could/should we loop over bdyLabelName stratas ? */ 115 ierr = DMLabelCreate("boundary", &bd);CHKERRQ(ierr); 116 ierr = DMPlexMarkBoundaryFaces(dm, bd);CHKERRQ(ierr); 117 ierr = DMLabelGetStratumIS(bd, 1, &bdIS);CHKERRQ(ierr); 118 ierr = DMLabelGetStratumSize(bd, 1, &numBdFaces);CHKERRQ(ierr); 119 ierr = ISGetIndices(bdIS, &faces);CHKERRQ(ierr); 120 /* TODO why not assume that we are considering simplicial meshes, in which case bdSize = dim*numBdFaces ? */ 121 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 122 PetscInt *closure = NULL; 123 PetscInt closureSize, cl; 124 125 ierr = DMPlexGetTransitiveClosure(dm, faces[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 126 for (cl = 0; cl < closureSize*2; cl += 2) { 127 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize; 128 } 129 ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 130 } 131 ierr = PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);CHKERRQ(ierr); 132 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 133 PetscInt *closure = NULL; 134 PetscInt closureSize, cl; 135 136 ierr = DMPlexGetTransitiveClosure(dm, faces[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 137 for (cl = 0; cl < closureSize*2; cl += 2) { 138 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart; 139 } 140 if (bdyLabelExt) { 141 ierr = DMLabelGetValue(bdTags, faces[f], &bdFaceIds[f]);CHKERRQ(ierr); 142 } else { 143 bdFaceIds[f] = 1; 144 } 145 ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 146 } 147 ierr = ISDestroy(&bdIS);CHKERRQ(ierr); 148 ierr = DMLabelDestroy(&bd);CHKERRQ(ierr); 149 pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds); 150 pragmatic_set_metric(met); 151 pragmatic_adapt(); 152 ierr = PetscFree5(x, y, z, cells, met);CHKERRQ(ierr); 153 ierr = PetscFree2(bdFaces, bdFaceIds);CHKERRQ(ierr); 154 ierr = PetscFree(coords);CHKERRQ(ierr); 155 ierr = DMDestroy(&udm);CHKERRQ(ierr); 156 ierr = DMDestroy(&coordDM);CHKERRQ(ierr); 157 /* Read out mesh */ 158 pragmatic_get_info(&numVerticesAdp, &numCellsAdp); 159 ierr = PetscMalloc1(numVerticesAdp*dim, &coordsAdp);CHKERRQ(ierr); 160 switch (dim) { 161 case 2: 162 ierr = PetscMalloc2(numVerticesAdp, &xAdp, numVerticesAdp, &yAdp);CHKERRQ(ierr); 163 zAdp = NULL; 164 pragmatic_get_coords_2d(xAdp, yAdp); 165 numCornersAdp = 3; 166 for (v = 0; v < numVerticesAdp; ++v) {coordsAdp[v*2+0] = xAdp[v]; coordsAdp[v*2+1] = yAdp[v];} 167 break; 168 case 3: 169 ierr = PetscMalloc3(numVerticesAdp, &xAdp, numVerticesAdp, &yAdp, numVerticesAdp, &zAdp);CHKERRQ(ierr); 170 pragmatic_get_coords_3d(xAdp, yAdp, zAdp); 171 numCornersAdp = 4; 172 for (v = 0; v < numVerticesAdp; ++v) {coordsAdp[v*3+0] = xAdp[v]; coordsAdp[v*3+1] = yAdp[v]; coordsAdp[v*3+2] = zAdp[v];} 173 break; 174 default: 175 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %d", dim); 176 } 177 ierr = PetscMalloc1(numCellsAdp*(dim+1), &cellsAdp);CHKERRQ(ierr); // only for simplicial meshes 178 pragmatic_get_elements(cellsAdp); 179 ierr = DMPlexCreate(comm, dmAdapted);CHKERRQ(ierr); 180 ierr = DMPlexCreateFromCellList(PetscObjectComm((PetscObject) dm), dim, numCellsAdp, numVerticesAdp, numCornersAdp, PETSC_TRUE, cellsAdp, dim, coordsAdp, dmAdapted);CHKERRQ(ierr); 181 /* Read out boundary tags */ 182 pragmatic_get_boundaryTags(&boundaryTags); 183 if (!bdyLabelExt) bdyLabelName = "boundary"; 184 ierr = DMCreateLabel(*dmAdapted, bdyLabelName);CHKERRQ(ierr); 185 ierr = DMGetLabel(*dmAdapted, bdyLabelName, &bdTagsAdp);CHKERRQ(ierr); 186 ierr = DMPlexGetHeightStratum(*dmAdapted, 0, &cStart, &cEnd);CHKERRQ(ierr); 187 ierr = DMPlexGetDepthStratum(*dmAdapted, 0, &vStart, &vEnd);CHKERRQ(ierr); 188 for (c = cStart; c < cEnd; ++c) { 189 PetscInt *cellClosure = NULL; 190 PetscInt cellClosureSize, cl; 191 PetscInt *facetClosure = NULL; 192 PetscInt facetClosureSize, cl2; 193 const PetscInt *facetList; 194 PetscInt facetListSize, f; 195 196 cellOff = (c-cStart)*(dim+1); // gives the offset corresponding to the cell in pragmatic boundary arrays. Only for simplicial meshes 197 hasBdyFacet = PETSC_FALSE; 198 for (i = 0; i < dim+1; ++i) { // only for simplicial meshes 199 if (boundaryTags[cellOff+i]) { 200 hasBdyFacet = PETSC_TRUE; 201 break; 202 } 203 } 204 if (!hasBdyFacet) continue; // The cell has no boundary edge/facet => no boundary tagging 205 206 cell = &cellsAdp[cellOff]; // pointing to the current cell in the cellsAdp table 207 ierr = DMPlexGetTransitiveClosure(*dmAdapted, c, PETSC_TRUE, &cellClosureSize, &cellClosure);CHKERRQ(ierr); 208 /* first, encode the permutation of the vertices indices betwenn the cell closure and pragmatic ordering */ 209 j = 0; 210 for (cl = 0; cl < cellClosureSize*2; cl += 2) { 211 if ((cellClosure[cl] < vStart) || (cellClosure[cl] >= vEnd)) continue; 212 iVer = cellClosure[cl] - vStart; 213 for (i = 0; i < dim+1; ++i) { // only for simplicial meshes 214 if (iVer == cell[i]) { 215 perm[j] = i; // the cl-th element of the closure is the i-th vertex of the cell 216 break; 217 } 218 } 219 j++; 220 } 221 ierr = DMPlexGetCone(*dmAdapted, c, &facetList);CHKERRQ(ierr); // list of edges/facets of the cell 222 ierr = DMPlexGetConeSize(*dmAdapted, c, &facetListSize);CHKERRQ(ierr); 223 /* 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) */ 224 for (f = 0; f < facetListSize; ++f){ 225 facet = facetList[f]; 226 ierr = DMPlexGetTransitiveClosure(*dmAdapted, facet, PETSC_TRUE, &facetClosureSize, &facetClosure);CHKERRQ(ierr); 227 /* loop over the vertices of the cell closure, and check if each vertex belongs to the edge/facet closure */ 228 /* TODO should we use bitmaps to mark vertices instead of a static array ? */ 229 PetscMemzero(isInFacetClosure, sizeof(isInFacetClosure)); 230 j = 0; 231 for (cl = 0; cl < cellClosureSize*2; cl += 2) { 232 if ((cellClosure[cl] < vStart) || (cellClosure[cl] >= vEnd)) continue; 233 for (cl2 = 0; cl2 < facetClosureSize*2; cl2 += 2){ 234 if ((facetClosure[cl2] < vStart) || (facetClosure[cl2] >= vEnd)) continue; 235 if (cellClosure[cl] == facetClosure[cl2]) { 236 isInFacetClosure[j] = 1; 237 } 238 } 239 j++; 240 } 241 /* 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 */ 242 j = 0; 243 for (cl = 0; cl < cellClosureSize*2; cl += 2) { 244 if ((cellClosure[cl] < vStart) || (cellClosure[cl] >= vEnd)) continue; 245 if (!isInFacetClosure[j]) { 246 idx = cellOff + perm[j]; 247 if (boundaryTags[idx]) { 248 ierr = DMLabelSetValue(bdTagsAdp, facet, boundaryTags[idx]);CHKERRQ(ierr); 249 } 250 break; 251 } 252 j++; 253 } 254 ierr = DMPlexRestoreTransitiveClosure(*dmAdapted, facet, PETSC_TRUE, &facetClosureSize, &facetClosure);CHKERRQ(ierr); 255 } 256 ierr = DMPlexRestoreTransitiveClosure(*dmAdapted, c, PETSC_TRUE, &cellClosureSize, &cellClosure);CHKERRQ(ierr); 257 } 258 pragmatic_finalize(); 259 ierr = PetscFree5(xAdp, yAdp, zAdp, cellsAdp, coordsAdp);CHKERRQ(ierr); 260 #endif 261 PetscFunctionReturn(0); 262 } 263