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