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 /* Read out mesh */ 156 pragmatic_get_info(&numVerticesAdp, &numCellsAdp); 157 ierr = PetscMalloc1(numVerticesAdp*dim, &coordsAdp);CHKERRQ(ierr); 158 switch (dim) { 159 case 2: 160 ierr = PetscMalloc2(numVerticesAdp, &xAdp, numVerticesAdp, &yAdp);CHKERRQ(ierr); 161 zAdp = NULL; 162 pragmatic_get_coords_2d(xAdp, yAdp); 163 numCornersAdp = 3; 164 for (v = 0; v < numVerticesAdp; ++v) {coordsAdp[v*2+0] = xAdp[v]; coordsAdp[v*2+1] = yAdp[v];} 165 break; 166 case 3: 167 ierr = PetscMalloc3(numVerticesAdp, &xAdp, numVerticesAdp, &yAdp, numVerticesAdp, &zAdp);CHKERRQ(ierr); 168 pragmatic_get_coords_3d(xAdp, yAdp, zAdp); 169 numCornersAdp = 4; 170 for (v = 0; v < numVerticesAdp; ++v) {coordsAdp[v*3+0] = xAdp[v]; coordsAdp[v*3+1] = yAdp[v]; coordsAdp[v*3+2] = zAdp[v];} 171 break; 172 default: 173 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %d", dim); 174 } 175 ierr = PetscMalloc1(numCellsAdp*(dim+1), &cellsAdp);CHKERRQ(ierr); // only for simplicial meshes 176 pragmatic_get_elements(cellsAdp); 177 ierr = DMPlexCreate(comm, dmAdapted);CHKERRQ(ierr); 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