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