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__ "DMPlexAdaptdm->coarseMesh" 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[28], isInFacetClosure[28], iVer, i, idx, facet; 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 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 92 PetscInt *closure = NULL; 93 PetscInt closureSize, cl; 94 95 ierr = DMPlexGetTransitiveClosure(dm, faces[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 96 for (cl = 0; cl < closureSize*2; cl += 2) { 97 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize; 98 } 99 ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 100 } 101 ierr = PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);CHKERRQ(ierr); 102 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 103 PetscInt *closure = NULL; 104 PetscInt closureSize, cl; 105 106 ierr = DMPlexGetTransitiveClosure(dm, faces[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 107 for (cl = 0; cl < closureSize*2; cl += 2) { 108 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart; 109 } 110 if ( bdyLabelName ) 111 DMLabelGetValue(bdtags, faces[f], &bdFaceIds[f]); 112 else 113 bdFaceIds[f] = 1; 114 ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 115 } 116 ierr = ISDestroy(&bdIS);CHKERRQ(ierr); 117 ierr = DMLabelDestroy(&bd);CHKERRQ(ierr); 118 pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds); 119 pragmatic_set_metric(met); 120 pragmatic_adapt(); 121 /* Read out mesh */ 122 pragmatic_get_info(&numCoarseVertices, &numCoarseCells); 123 ierr = PetscMalloc1(numCoarseVertices*dim, &coarseCoords);CHKERRQ(ierr); 124 switch (dim) { 125 case 2: 126 ierr = PetscMalloc2(numCoarseVertices, &x, numCoarseVertices, &y);CHKERRQ(ierr); 127 pragmatic_get_coords_2d(x, y); 128 numCorners = 3; 129 for (v = 0; v < numCoarseVertices; ++v) {coarseCoords[v*2+0] = x[v]; coarseCoords[v*2+1] = y[v];} 130 break; 131 case 3: 132 ierr = PetscMalloc3(numCoarseVertices, &x, numCoarseVertices, &y, numCoarseVertices, &z);CHKERRQ(ierr); 133 pragmatic_get_coords_3d(x, y, z); 134 numCorners = 4; 135 for (v = 0; v < numCoarseVertices; ++v) {coarseCoords[v*3+0] = x[v]; coarseCoords[v*3+1] = y[v]; coarseCoords[v*3+2] = z[v];} 136 break; 137 default: 138 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic coarsening defined for dimension %d", dim); 139 } 140 ierr = PetscMalloc1(numCoarseCells*(dim+1), &ccells);CHKERRQ(ierr); // only for simplicial meshes 141 pragmatic_get_elements(ccells); 142 ierr = DMPlexCreateFromCellList(PetscObjectComm((PetscObject) dm), dim, numCoarseCells, numCoarseVertices, numCorners, PETSC_TRUE, ccells, dim, coarseCoords, &dm->coarseMesh);CHKERRQ(ierr); 143 /* Read out boundary tags */ 144 pragmatic_get_boundaryTags(&boundaryTags); 145 if ( !bdyLabelName ) bdyLabelName = "boundary"; 146 ierr = DMCreateLabel(dm->coarseMesh, bdyLabelName);;CHKERRQ(ierr); 147 ierr = DMPlexGetHeightStratum(dm->coarseMesh, 0, &cStart, &cEnd);CHKERRQ(ierr); 148 ierr = DMPlexGetDepthStratum(dm->coarseMesh, 0, &vStart, &vEnd);CHKERRQ(ierr); 149 for (c = cStart; c < cEnd; ++c) { 150 PetscInt *cellClosure = NULL; 151 PetscInt cellClosureSize, cl; 152 PetscInt *facetClosure = NULL; 153 PetscInt facetClosureSize, cl2; 154 const PetscInt *facetList; 155 PetscInt facetListSize, f; 156 157 cell = &ccells[(c-cStart)*(dim+1)]; // pointing to the current cell in the ccell table 158 ierr = DMPlexGetTransitiveClosure(dm->coarseMesh, c, PETSC_TRUE, &cellClosureSize, &cellClosure);CHKERRQ(ierr); 159 // first, encode the permutation of the vertices indices betwenn the cell closure and pragmatic ordering 160 for (cl = 0; cl < cellClosureSize*2; cl+=2 ) { 161 if ((cellClosure[cl] < vStart) || (cellClosure[cl] >= vEnd)) continue; 162 iVer = cellClosure[cl] - vStart; 163 for (i=0; i<dim+1; ++i) { 164 if ( iVer == cell[i] ) { 165 perm[cl] = i; // the cl-th vertex of the closure is the i-th vertex of the ccell 166 break; 167 } 168 } 169 } 170 ierr = DMPlexGetCone(dm->coarseMesh, c, &facetList);CHKERRQ(ierr); // list of edges/facets of the cell 171 ierr = DMPlexGetConeSize(dm->coarseMesh, c, &facetListSize);CHKERRQ(ierr); 172 // then, for each edge/facet of the cell, find the opposite vertex (ie the one not in the closure of the facet/edge) 173 for (f=0; f<facetListSize; ++f){ 174 facet = facetList[f]; 175 ierr = DMPlexGetTransitiveClosure(dm->coarseMesh, facet, PETSC_TRUE, &facetClosureSize, &facetClosure);CHKERRQ(ierr); 176 // loop over the vertices of the cell closure, and check if each vertex belongs to the edge/facet closure 177 PetscMemzero(isInFacetClosure, sizeof(isInFacetClosure)); 178 for (cl = 0; cl < cellClosureSize*2; cl+=2) { 179 if ((cellClosure[cl] < vStart) || (cellClosure[cl] >= vEnd)) isInFacetClosure[cl] = 1; 180 for (cl2 =0; cl2<facetClosureSize*2; cl2+=2){ 181 if ((facetClosure[cl2] < vStart) || (facetClosure[cl2] >= vEnd)) continue; 182 if ( cellClosure[cl] == facetClosure[cl2] ) { 183 isInFacetClosure[cl] = 1; 184 } 185 } 186 } 187 // 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 188 for (cl = 0; cl < cellClosureSize*2; cl+=2 ) { 189 if ( !isInFacetClosure[cl] ) { 190 idx = (c-cStart)*(dim+1) + perm[cl]; 191 if ( boundaryTags[idx] ) 192 ierr = DMSetLabelValue(dm->coarseMesh, bdyLabelName, facet, boundaryTags[idx]);CHKERRQ(ierr); 193 break; 194 } 195 } 196 ierr = DMPlexRestoreTransitiveClosure(dm->coarseMesh, facet, PETSC_TRUE, &facetClosureSize, &facetClosure);CHKERRQ(ierr); 197 } 198 ierr = DMPlexRestoreTransitiveClosure(dm->coarseMesh, c, PETSC_TRUE, &cellClosureSize, &cellClosure);CHKERRQ(ierr); 199 } 200 pragmatic_finalize(); 201 ierr = PetscFree4(x, y, z, cells);CHKERRQ(ierr); 202 ierr = PetscFree2(bdFaces, bdFaceIds);CHKERRQ(ierr); 203 ierr = PetscFree(coarseCoords);CHKERRQ(ierr); 204 } 205 #endif 206 ierr = PetscObjectReference((PetscObject) dm->coarseMesh);CHKERRQ(ierr); 207 *dmCoarsened = dm->coarseMesh; 208 PetscFunctionReturn(0); 209 } 210