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 ierr = PetscFree5(x, y, z, cells, met);CHKERRQ(ierr); 152 ierr = PetscFree2(bdFaces, bdFaceIds);CHKERRQ(ierr); 153 ierr = PetscFree(coords);CHKERRQ(ierr); 154 ierr = DMDestroy(&udm);CHKERRQ(ierr); 155 pragmatic_adapt(); 156 /* Read out mesh */ 157 pragmatic_get_info(&numVerticesAdp, &numCellsAdp); 158 ierr = PetscMalloc1(numVerticesAdp*dim, &coordsAdp);CHKERRQ(ierr); 159 switch (dim) { 160 case 2: 161 ierr = PetscMalloc2(numVerticesAdp, &xAdp, numVerticesAdp, &yAdp);CHKERRQ(ierr); 162 zAdp = NULL; 163 pragmatic_get_coords_2d(xAdp, yAdp); 164 numCornersAdp = 3; 165 for (v = 0; v < numVerticesAdp; ++v) {coordsAdp[v*2+0] = xAdp[v]; coordsAdp[v*2+1] = yAdp[v];} 166 break; 167 case 3: 168 ierr = PetscMalloc3(numVerticesAdp, &xAdp, numVerticesAdp, &yAdp, numVerticesAdp, &zAdp);CHKERRQ(ierr); 169 pragmatic_get_coords_3d(xAdp, yAdp, zAdp); 170 numCornersAdp = 4; 171 for (v = 0; v < numVerticesAdp; ++v) {coordsAdp[v*3+0] = xAdp[v]; coordsAdp[v*3+1] = yAdp[v]; coordsAdp[v*3+2] = zAdp[v];} 172 break; 173 default: 174 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %d", dim); 175 } 176 ierr = PetscMalloc1(numCellsAdp*(dim+1), &cellsAdp);CHKERRQ(ierr); // only for simplicial meshes 177 pragmatic_get_elements(cellsAdp); 178 ierr = DMPlexCreateFromCellList(PetscObjectComm((PetscObject) dm), dim, numCellsAdp, numVerticesAdp, numCornersAdp, PETSC_TRUE, cellsAdp, dim, coordsAdp, dmAdapted);CHKERRQ(ierr); 179 /* Read out boundary tags */ 180 pragmatic_get_boundaryTags(&boundaryTags); 181 if (!bdyLabelExt) bdyLabelName = "boundary"; 182 ierr = DMCreateLabel(*dmAdapted, bdyLabelName);CHKERRQ(ierr); 183 ierr = DMGetLabel(*dmAdapted, bdyLabelName, &bdTagsAdp);CHKERRQ(ierr); 184 ierr = DMPlexGetHeightStratum(*dmAdapted, 0, &cStart, &cEnd);CHKERRQ(ierr); 185 ierr = DMPlexGetDepthStratum(*dmAdapted, 0, &vStart, &vEnd);CHKERRQ(ierr); 186 for (c = cStart; c < cEnd; ++c) { 187 PetscInt *cellClosure = NULL; 188 PetscInt cellClosureSize, cl; 189 PetscInt *facetClosure = NULL; 190 PetscInt facetClosureSize, cl2; 191 const PetscInt *facetList; 192 PetscInt facetListSize, f; 193 194 cellOff = (c-cStart)*(dim+1); // gives the offset corresponding to the cell in pragmatic boundary arrays. Only for simplicial meshes 195 hasBdyFacet = PETSC_FALSE; 196 for (i = 0; i < dim+1; ++i) { // only for simplicial meshes 197 if (boundaryTags[cellOff+i]) { 198 hasBdyFacet = PETSC_TRUE; 199 break; 200 } 201 } 202 if (!hasBdyFacet) continue; // The cell has no boundary edge/facet => no boundary tagging 203 204 cell = &cellsAdp[cellOff]; // pointing to the current cell in the cellsAdp table 205 ierr = DMPlexGetTransitiveClosure(*dmAdapted, c, PETSC_TRUE, &cellClosureSize, &cellClosure);CHKERRQ(ierr); 206 /* first, encode the permutation of the vertices indices betwenn the cell closure and pragmatic ordering */ 207 j = 0; 208 for (cl = 0; cl < cellClosureSize*2; cl += 2) { 209 if ((cellClosure[cl] < vStart) || (cellClosure[cl] >= vEnd)) continue; 210 iVer = cellClosure[cl] - vStart; 211 for (i = 0; i < dim+1; ++i) { // only for simplicial meshes 212 if (iVer == cell[i]) { 213 perm[j] = i; // the cl-th element of the closure is the i-th vertex of the cell 214 break; 215 } 216 } 217 j++; 218 } 219 ierr = DMPlexGetCone(*dmAdapted, c, &facetList);CHKERRQ(ierr); // list of edges/facets of the cell 220 ierr = DMPlexGetConeSize(*dmAdapted, c, &facetListSize);CHKERRQ(ierr); 221 /* 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) */ 222 for (f = 0; f < facetListSize; ++f){ 223 facet = facetList[f]; 224 ierr = DMPlexGetTransitiveClosure(*dmAdapted, facet, PETSC_TRUE, &facetClosureSize, &facetClosure);CHKERRQ(ierr); 225 /* loop over the vertices of the cell closure, and check if each vertex belongs to the edge/facet closure */ 226 /* TODO should we use bitmaps to mark vertices instead of a static array ? */ 227 PetscMemzero(isInFacetClosure, sizeof(isInFacetClosure)); 228 j = 0; 229 for (cl = 0; cl < cellClosureSize*2; cl += 2) { 230 if ((cellClosure[cl] < vStart) || (cellClosure[cl] >= vEnd)) continue; 231 for (cl2 = 0; cl2 < facetClosureSize*2; cl2 += 2){ 232 if ((facetClosure[cl2] < vStart) || (facetClosure[cl2] >= vEnd)) continue; 233 if (cellClosure[cl] == facetClosure[cl2]) { 234 isInFacetClosure[j] = 1; 235 } 236 } 237 j++; 238 } 239 /* 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 */ 240 j = 0; 241 for (cl = 0; cl < cellClosureSize*2; cl += 2) { 242 if ((cellClosure[cl] < vStart) || (cellClosure[cl] >= vEnd)) continue; 243 if (!isInFacetClosure[j]) { 244 idx = cellOff + perm[j]; 245 if (boundaryTags[idx]) { 246 ierr = DMLabelSetValue(bdTagsAdp, facet, boundaryTags[idx]);CHKERRQ(ierr); 247 } 248 break; 249 } 250 j++; 251 } 252 ierr = DMPlexRestoreTransitiveClosure(*dmAdapted, facet, PETSC_TRUE, &facetClosureSize, &facetClosure);CHKERRQ(ierr); 253 } 254 ierr = DMPlexRestoreTransitiveClosure(*dmAdapted, c, PETSC_TRUE, &cellClosureSize, &cellClosure);CHKERRQ(ierr); 255 } 256 pragmatic_finalize(); 257 ierr = PetscFree5(xAdp, yAdp, zAdp, cellsAdp, coordsAdp);CHKERRQ(ierr); 258 #endif 259 PetscFunctionReturn(0); 260 } 261