1 #include <petscsf.h> 2 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 3 4 PETSC_STATIC_INLINE PetscErrorCode FillClosureArray_Static(DM dm, PetscSection section, PetscInt nP, const PetscInt points[], PetscScalar *vArray, PetscInt *csize, PetscScalar **array) 5 { 6 PetscScalar *a; 7 PetscInt pStart, pEnd, size = 0, dof, off, d, k, i; 8 PetscErrorCode ierr; 9 10 PetscFunctionBegin; 11 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr); 12 for (i = 0; i < nP; ++i) { 13 const PetscInt p = points[i]; 14 15 if ((p < pStart) || (p >= pEnd)) continue; 16 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 17 size += dof; 18 } 19 if (csize) *csize = size; 20 if (array) { 21 ierr = DMGetWorkArray(dm, size, PETSC_SCALAR, &a);CHKERRQ(ierr); 22 for (i = 0, k = 0; i < nP; ++i) { 23 const PetscInt p = points[i]; 24 25 if ((p < pStart) || (p >= pEnd)) continue; 26 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 27 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 28 29 for (d = 0; d < dof; ++d, ++k) a[k] = vArray[off+d]; 30 } 31 *array = a; 32 } 33 PetscFunctionReturn(0); 34 } 35 36 PETSC_STATIC_INLINE PetscErrorCode FillClosureVec_Private(DM dm, PetscSection section, PetscInt nP, const PetscInt points[], PetscScalar *vArray, const PetscScalar *array, InsertMode mode) 37 { 38 PetscInt dof, off, d, k, i; 39 PetscErrorCode ierr; 40 41 PetscFunctionBegin; 42 if ((mode == INSERT_VALUES) || (mode == INSERT_ALL_VALUES)) { 43 for (i = 0, k = 0; i < nP; ++i) { 44 ierr = PetscSectionGetDof(section, points[i], &dof);CHKERRQ(ierr); 45 ierr = PetscSectionGetOffset(section, points[i], &off);CHKERRQ(ierr); 46 47 for (d = 0; d < dof; ++d, ++k) vArray[off+d] = array[k]; 48 } 49 } else { 50 for (i = 0, k = 0; i < nP; ++i) { 51 ierr = PetscSectionGetDof(section, points[i], &dof);CHKERRQ(ierr); 52 ierr = PetscSectionGetOffset(section, points[i], &off);CHKERRQ(ierr); 53 54 for (d = 0; d < dof; ++d, ++k) vArray[off+d] += array[k]; 55 } 56 } 57 PetscFunctionReturn(0); 58 } 59 60 PETSC_STATIC_INLINE PetscErrorCode GetPointArray_Private(DM dm,PetscInt n,PetscInt *points,PetscInt *rn,const PetscInt **rpoints) 61 { 62 PetscErrorCode ierr; 63 PetscInt *work; 64 65 PetscFunctionBegin; 66 if (rn) *rn = n; 67 if (rpoints) { 68 ierr = DMGetWorkArray(dm,n,PETSC_INT,&work);CHKERRQ(ierr); 69 ierr = PetscMemcpy(work,points,n*sizeof(PetscInt));CHKERRQ(ierr); 70 *rpoints = work; 71 } 72 PetscFunctionReturn(0); 73 } 74 75 PETSC_STATIC_INLINE PetscErrorCode RestorePointArray_Private(DM dm,PetscInt *rn,const PetscInt **rpoints) 76 { 77 PetscErrorCode ierr; 78 79 PetscFunctionBegin; 80 if (rn) *rn = 0; 81 if (rpoints) { 82 /* note the second argument to DMRestoreWorkArray() is always ignored */ 83 ierr = DMRestoreWorkArray(dm,0, PETSC_INT, (void*) rpoints);CHKERRQ(ierr); 84 } 85 PetscFunctionReturn(0); 86 } 87 88 PetscErrorCode DMDAGetTransitiveClosure(DM dm, PetscInt p, PetscBool useClosure, PetscInt *closureSize, PetscInt **closure) 89 { 90 PetscInt dim = dm->dim; 91 PetscInt nVx, nVy, nVz, nxF, nXF, nyF, nYF, nzF, nZF; 92 PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart; 93 PetscErrorCode ierr; 94 95 PetscFunctionBegin; 96 if (!useClosure) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Star operation is not yet supported"); 97 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 98 PetscValidIntPointer(closureSize,4); 99 ierr = DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);CHKERRQ(ierr); 100 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 101 ierr = DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 102 ierr = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr); 103 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, NULL);CHKERRQ(ierr); 104 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 105 xfStart = fStart; xfEnd = xfStart+nXF; 106 yfStart = xfEnd; 107 if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd); 108 if ((p >= cStart) || (p < cEnd)) { 109 /* Cell */ 110 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 111 else if (dim == 2) { 112 /* 4 faces, 4 vertices 113 Bottom-left vertex follows same order as cells 114 Bottom y-face same order as cells 115 Left x-face follows same order as cells 116 We number the quad: 117 118 8--3--7 119 | | 120 4 0 2 121 | | 122 5--1--6 123 */ 124 PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1); 125 PetscInt v = cy*nVx + cx + vStart; 126 PetscInt xf = cy*nxF + cx + xfStart; 127 PetscInt yf = c + yfStart; 128 129 *closureSize = 9; 130 if (!(*closure)) {ierr = DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);CHKERRQ(ierr);} 131 (*closure)[0] = p; (*closure)[1] = yf; (*closure)[2] = xf+1; (*closure)[3] = yf+nyF; (*closure)[4] = xf+0; (*closure)[5] = v+0; (*closure)[6]= v+1; (*closure)[7] = v+nVx+1; (*closure)[8] = v+nVx+0; 132 } else { 133 /* 6 faces, 12 edges, 8 vertices 134 Bottom-left vertex follows same order as cells 135 Bottom y-face same order as cells 136 Left x-face follows same order as cells 137 We number the quad: 138 139 8--3--7 140 | | 141 4 0 2 142 | | 143 5--1--6 144 */ 145 #if 0 146 PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1), cz = c / ((nVx-1)*(nVy-1)); 147 PetscInt v = cy*nVx + cx + vStart; 148 PetscInt xf = cy*nxF + cx + xfStart; 149 PetscInt yf = c + yfStart; 150 151 *closureSize = 26; 152 if (!(*closure)) {ierr = DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);CHKERRQ(ierr);} 153 #endif 154 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 155 } 156 } else if ((p >= vStart) || (p < vEnd)) { 157 /* Vertex */ 158 *closureSize = 1; 159 if (!(*closure)) {ierr = DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);CHKERRQ(ierr);} 160 (*closure)[0] = p; 161 } else if ((p >= fStart) || (p < fStart + nXF)) { 162 /* X Face */ 163 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 164 else if (dim == 2) { 165 /* 2 vertices: The bottom vertex has the same numbering as the face */ 166 PetscInt f = p - xfStart; 167 168 *closureSize = 3; 169 if (!(*closure)) {ierr = DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);CHKERRQ(ierr);} 170 (*closure)[0] = p; (*closure)[1] = f; (*closure)[2] = f+nVx; 171 } else if (dim == 3) { 172 /* 4 vertices */ 173 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 174 } 175 } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) { 176 /* Y Face */ 177 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 178 else if (dim == 2) { 179 /* 2 vertices: The left vertex has the same numbering as the face */ 180 PetscInt f = p - yfStart; 181 182 *closureSize = 3; 183 if (!(*closure)) {ierr = DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);CHKERRQ(ierr);} 184 (*closure)[0] = p; (*closure)[1] = f; (*closure)[2]= f+1; 185 } else if (dim == 3) { 186 /* 4 vertices */ 187 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 188 } 189 } else { 190 /* Z Face */ 191 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 192 else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D"); 193 else if (dim == 3) { 194 /* 4 vertices */ 195 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 196 } 197 } 198 PetscFunctionReturn(0); 199 } 200 201 PetscErrorCode DMDARestoreTransitiveClosure(DM dm, PetscInt p, PetscBool useClosure, PetscInt *closureSize, PetscInt **closure) 202 { 203 PetscErrorCode ierr; 204 205 PetscFunctionBegin; 206 ierr = DMRestoreWorkArray(dm, 0, PETSC_INT, closure);CHKERRQ(ierr); 207 PetscFunctionReturn(0); 208 } 209 210 PetscErrorCode DMDAGetClosure(DM dm, PetscSection section, PetscInt p,PetscInt *n,const PetscInt **closure) 211 { 212 PetscInt dim = dm->dim; 213 PetscInt nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF; 214 PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart; 215 PetscErrorCode ierr; 216 217 PetscFunctionBegin; 218 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 219 if (n) PetscValidIntPointer(n,4); 220 if (closure) PetscValidPointer(closure, 5); 221 if (!section) {ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr);} 222 if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection"); 223 ierr = DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);CHKERRQ(ierr); 224 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 225 ierr = DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 226 ierr = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr); 227 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);CHKERRQ(ierr); 228 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 229 xfStart = fStart; xfEnd = xfStart+nXF; 230 yfStart = xfEnd; 231 if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd); 232 if ((p >= cStart) || (p < cEnd)) { 233 /* Cell */ 234 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 235 else if (dim == 2) { 236 /* 4 faces, 4 vertices 237 Bottom-left vertex follows same order as cells 238 Bottom y-face same order as cells 239 Left x-face follows same order as cells 240 We number the quad: 241 242 8--3--7 243 | | 244 4 0 2 245 | | 246 5--1--6 247 */ 248 PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1); 249 PetscInt v = cy*nVx + cx + vStart; 250 PetscInt xf = cy*nxF + cx + xfStart; 251 PetscInt yf = c + yfStart; 252 PetscInt points[9]; 253 254 /* Note: initializer list is not computable at compile time, hence we fill the array manually */ 255 points[0] = p; points[1] = yf; points[2] = xf+1; points[3] = yf+nyF; points[4] = xf+0; points[5] = v+0; points[6]= v+1; points[7] = v+nVx+1; points[8] = v+nVx+0; 256 257 ierr = GetPointArray_Private(dm,9,points,n,closure);CHKERRQ(ierr); 258 } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 259 } else if ((p >= vStart) || (p < vEnd)) { 260 /* Vertex */ 261 ierr = GetPointArray_Private(dm,1,&p,n,closure);CHKERRQ(ierr); 262 } else if ((p >= fStart) || (p < fStart + nXF)) { 263 /* X Face */ 264 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 265 else if (dim == 2) { 266 /* 2 vertices: The bottom vertex has the same numbering as the face */ 267 PetscInt f = p - xfStart; 268 PetscInt points[3]; 269 270 points[0] = p; points[1] = f; points[2] = f+nVx; 271 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken"); 272 ierr = GetPointArray_Private(dm,3,points,n,closure);CHKERRQ(ierr); 273 } else if (dim == 3) { 274 /* 4 vertices */ 275 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 276 } 277 } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) { 278 /* Y Face */ 279 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 280 else if (dim == 2) { 281 /* 2 vertices: The left vertex has the same numbering as the face */ 282 PetscInt f = p - yfStart; 283 PetscInt points[3]; 284 285 points[0] = p; points[1] = f; points[2]= f+1; 286 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken"); 287 ierr = GetPointArray_Private(dm, 3, points, n, closure);CHKERRQ(ierr); 288 } else if (dim == 3) { 289 /* 4 vertices */ 290 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 291 } 292 } else { 293 /* Z Face */ 294 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 295 else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D"); 296 else if (dim == 3) { 297 /* 4 vertices */ 298 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 299 } 300 } 301 PetscFunctionReturn(0); 302 } 303 304 /* Zeros n and closure. */ 305 PetscErrorCode DMDARestoreClosure(DM dm, PetscSection section, PetscInt p,PetscInt *n,const PetscInt **closure) 306 { 307 PetscErrorCode ierr; 308 309 PetscFunctionBegin; 310 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 311 if (n) PetscValidIntPointer(n,4); 312 if (closure) PetscValidPointer(closure, 5); 313 ierr = RestorePointArray_Private(dm,n,closure);CHKERRQ(ierr); 314 PetscFunctionReturn(0); 315 } 316 317 PetscErrorCode DMDAGetClosureScalar(DM dm, PetscSection section, PetscInt p, PetscScalar *vArray, PetscInt *csize, PetscScalar **values) 318 { 319 PetscInt dim = dm->dim; 320 PetscInt nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF; 321 PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart; 322 PetscErrorCode ierr; 323 324 PetscFunctionBegin; 325 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 326 PetscValidScalarPointer(vArray, 4); 327 if (values) PetscValidPointer(values, 6); 328 if (!section) {ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr);} 329 if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection"); 330 ierr = DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);CHKERRQ(ierr); 331 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 332 ierr = DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 333 ierr = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr); 334 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);CHKERRQ(ierr); 335 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 336 xfStart = fStart; xfEnd = xfStart+nXF; 337 yfStart = xfEnd; yfEnd = yfStart+nYF; 338 zfStart = yfEnd; 339 if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd); 340 if ((p >= cStart) || (p < cEnd)) { 341 /* Cell */ 342 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 343 else if (dim == 2) { 344 /* 4 faces, 4 vertices 345 Bottom-left vertex follows same order as cells 346 Bottom y-face same order as cells 347 Left x-face follows same order as cells 348 We number the quad: 349 350 8--3--7 351 | | 352 4 0 2 353 | | 354 5--1--6 355 */ 356 PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1); 357 PetscInt v = cy*nVx + cx + vStart; 358 PetscInt xf = cx*nxF + cy + xfStart; 359 PetscInt yf = c + yfStart; 360 PetscInt points[9]; 361 362 points[0] = p; 363 points[1] = yf; points[2] = xf+nxF; points[3] = yf+nyF; points[4] = xf; 364 points[5] = v+0; points[6] = v+1; points[7] = v+nVx+1; points[8] = v+nVx+0; 365 ierr = FillClosureArray_Static(dm, section, 9, points, vArray, csize, values);CHKERRQ(ierr); 366 } else { 367 /* 6 faces, 8 vertices 368 Bottom-left-back vertex follows same order as cells 369 Back z-face follows same order as cells 370 Bottom y-face follows same order as cells 371 Left x-face follows same order as cells 372 373 14-----13 374 /| /| 375 / | 2 / | 376 / 5| / | 377 10-----9 4| 378 | 11-|---12 379 |6 / | / 380 | /1 3| / 381 |/ |/ 382 7-----8 383 */ 384 PetscInt c = p - cStart; 385 PetscInt points[15]; 386 387 points[0] = p; points[1] = c+zfStart; points[2] = c+zfStart+nzF; points[3] = c+yfStart; points[4] = c+xfStart+nxF; points[5] = c+yfStart+nyF; points[6] = c+xfStart; 388 points[7] = c+vStart+0; points[8] = c+vStart+1; points[9] = c+vStart+nVx+1; points[10] = c+vStart+nVx+0; points[11] = c+vStart+nVx*nVy+0; 389 points[12] = c+vStart+nVx*nVy+1; points[13] = c+vStart+nVx*nVy+nVx+1; points[14] = c+vStart+nVx*nVy+nVx+0; 390 391 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken"); 392 ierr = FillClosureArray_Static(dm, section, 15, points, vArray, csize, values);CHKERRQ(ierr); 393 } 394 } else if ((p >= vStart) || (p < vEnd)) { 395 /* Vertex */ 396 ierr = FillClosureArray_Static(dm, section, 1, &p, vArray, csize, values);CHKERRQ(ierr); 397 } else if ((p >= fStart) || (p < fStart + nXF)) { 398 /* X Face */ 399 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 400 else if (dim == 2) { 401 /* 2 vertices: The bottom vertex has the same numbering as the face */ 402 PetscInt f = p - xfStart; 403 PetscInt points[3]; 404 405 points[0] = p; points[1] = f; points[2] = f+nVx; 406 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken"); 407 ierr = FillClosureArray_Static(dm, section, 3, points, vArray, csize, values);CHKERRQ(ierr); 408 } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 409 } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) { 410 /* Y Face */ 411 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 412 else if (dim == 2) { 413 /* 2 vertices: The left vertex has the same numbering as the face */ 414 PetscInt f = p - yfStart; 415 PetscInt points[3]; 416 417 points[0] = p; points[1] = f; points[2] = f+1; 418 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken"); 419 ierr = FillClosureArray_Static(dm, section, 3, points, vArray, csize, values);CHKERRQ(ierr); 420 } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 421 } else { 422 /* Z Face */ 423 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 424 else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D"); 425 else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 426 } 427 PetscFunctionReturn(0); 428 } 429 430 PetscErrorCode DMDAVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt p, PetscInt *csize, PetscScalar **values) 431 { 432 PetscScalar *vArray; 433 PetscErrorCode ierr; 434 435 PetscFunctionBegin; 436 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 437 PetscValidHeaderSpecific(v, VEC_CLASSID, 3); 438 if (values) PetscValidPointer(values, 6); 439 ierr = VecGetArray(v, &vArray);CHKERRQ(ierr); 440 ierr = DMDAGetClosureScalar(dm, section, p, vArray, csize, values);CHKERRQ(ierr); 441 ierr = VecRestoreArray(v, &vArray);CHKERRQ(ierr); 442 PetscFunctionReturn(0); 443 } 444 445 PetscErrorCode DMDARestoreClosureScalar(DM dm, PetscSection section, PetscInt p, PetscScalar *vArray, PetscInt *csize, PetscScalar **values) 446 { 447 PetscErrorCode ierr; 448 449 PetscFunctionBegin; 450 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 451 PetscValidPointer(values, 6); 452 ierr = DMRestoreWorkArray(dm, csize ? *csize : 0, PETSC_SCALAR, (void*) values);CHKERRQ(ierr); 453 PetscFunctionReturn(0); 454 } 455 456 PetscErrorCode DMDAVecRestoreClosure(DM dm, PetscSection section, Vec v, PetscInt p, PetscInt *csize, PetscScalar **values) 457 { 458 PetscErrorCode ierr; 459 460 PetscFunctionBegin; 461 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 462 PetscValidHeaderSpecific(v, VEC_CLASSID, 3); 463 PetscValidPointer(values, 5); 464 ierr = DMDARestoreClosureScalar(dm, section, p, NULL, csize, values);CHKERRQ(ierr); 465 PetscFunctionReturn(0); 466 } 467 468 PetscErrorCode DMDASetClosureScalar(DM dm, PetscSection section, PetscInt p,PetscScalar *vArray, const PetscScalar *values, InsertMode mode) 469 { 470 PetscInt dim = dm->dim; 471 PetscInt nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF, nCx, nCy; 472 PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart; 473 PetscErrorCode ierr; 474 475 PetscFunctionBegin; 476 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 477 PetscValidScalarPointer(values, 4); 478 PetscValidPointer(values, 5); 479 if (!section) {ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr);} 480 if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection"); 481 ierr = DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);CHKERRQ(ierr); 482 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 483 ierr = DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 484 ierr = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr); 485 ierr = DMDAGetNumCells(dm, &nCx, &nCy, NULL, NULL);CHKERRQ(ierr); 486 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);CHKERRQ(ierr); 487 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 488 xfStart = fStart; xfEnd = xfStart+nXF; 489 yfStart = xfEnd; yfEnd = yfStart+nYF; 490 zfStart = yfEnd; 491 if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd); 492 if ((p >= cStart) || (p < cEnd)) { 493 /* Cell */ 494 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 495 else if (dim == 2) { 496 /* 4 faces, 4 vertices 497 Bottom-left vertex follows same order as cells 498 Bottom y-face same order as cells 499 Left x-face follows same order as cells 500 We number the quad: 501 502 8--3--7 503 | | 504 4 0 2 505 | | 506 5--1--6 507 */ 508 PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1); 509 PetscInt v = cy*nVx + cx + vStart; 510 PetscInt xf = cx*nxF + cy + xfStart; 511 PetscInt yf = c + yfStart; 512 PetscInt points[9]; 513 514 points[0] = p; 515 points[1] = yf; points[2] = xf+nxF; points[3] = yf+nyF; points[4] = xf; 516 points[5] = v+0; points[6] = v+1; points[7] = v+nVx+1; points[8] = v+nVx+0; 517 ierr = FillClosureVec_Private(dm, section, 9, points, vArray, values, mode);CHKERRQ(ierr); 518 } else { 519 /* 6 faces, 8 vertices 520 Bottom-left-back vertex follows same order as cells 521 Back z-face follows same order as cells 522 Bottom y-face follows same order as cells 523 Left x-face follows same order as cells 524 525 14-----13 526 /| /| 527 / | 2 / | 528 / 5| / | 529 10-----9 4| 530 | 11-|---12 531 |6 / | / 532 | /1 3| / 533 |/ |/ 534 7-----8 535 */ 536 PetscInt c = p - cStart; 537 PetscInt points[15]; 538 539 points[0] = p; points[1] = c+zfStart; points[2] = c+zfStart+nzF; points[3] = c+yfStart; points[4] = c+xfStart+nxF; points[5] = c+yfStart+nyF; points[6] = c+xfStart; 540 points[7] = c+vStart+0; points[8] = c+vStart+1; points[9] = c+vStart+nVx+1; points[10] = c+vStart+nVx+0; points[11] = c+vStart+nVx*nVy+0; points[12] = c+vStart+nVx*nVy+1; 541 points[13] = c+vStart+nVx*nVy+nVx+1; points[14] = c+vStart+nVx*nVy+nVx+0; 542 ierr = FillClosureVec_Private(dm, section, 15, points, vArray, values, mode);CHKERRQ(ierr); 543 } 544 } else if ((p >= vStart) || (p < vEnd)) { 545 /* Vertex */ 546 ierr = FillClosureVec_Private(dm, section, 1, &p, vArray, values, mode);CHKERRQ(ierr); 547 } else if ((p >= fStart) || (p < fStart + nXF)) { 548 /* X Face */ 549 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 550 else if (dim == 2) { 551 /* 2 vertices: The bottom vertex has the same numbering as the face */ 552 PetscInt f = p - xfStart; 553 PetscInt points[3]; 554 555 points[0] = p; points[1] = f; points[2] = f+nVx; 556 ierr = FillClosureVec_Private(dm, section, 3, points, vArray, values, mode);CHKERRQ(ierr); 557 } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 558 } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) { 559 /* Y Face */ 560 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 561 else if (dim == 2) { 562 /* 2 vertices: The left vertex has the same numbering as the face */ 563 PetscInt f = p - yfStart; 564 PetscInt points[3]; 565 566 points[0] = p; points[1] = f; points[2] = f+1; 567 ierr = FillClosureVec_Private(dm, section, 3, points, vArray, values, mode);CHKERRQ(ierr); 568 } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 569 } else { 570 /* Z Face */ 571 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 572 else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D"); 573 else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 574 } 575 PetscFunctionReturn(0); 576 } 577 578 PetscErrorCode DMDAVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt p, const PetscScalar *values, InsertMode mode) 579 { 580 PetscScalar *vArray; 581 PetscErrorCode ierr; 582 583 PetscFunctionBegin; 584 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 585 PetscValidHeaderSpecific(v, VEC_CLASSID, 3); 586 PetscValidPointer(values, 5); 587 ierr = VecGetArray(v,&vArray);CHKERRQ(ierr); 588 ierr = DMDASetClosureScalar(dm,section,p,vArray,values,mode);CHKERRQ(ierr); 589 ierr = VecRestoreArray(v,&vArray);CHKERRQ(ierr); 590 PetscFunctionReturn(0); 591 } 592 593 /*@ 594 DMDAConvertToCell - Convert (i,j,k) to local cell number 595 596 Not Collective 597 598 Input Parameter: 599 + da - the distributed array 600 - s - A MatStencil giving (i,j,k) 601 602 Output Parameter: 603 . cell - the local cell number 604 605 Level: developer 606 607 .seealso: DMDAVecGetClosure() 608 @*/ 609 PetscErrorCode DMDAConvertToCell(DM dm, MatStencil s, PetscInt *cell) 610 { 611 DM_DA *da = (DM_DA*) dm->data; 612 const PetscInt dim = dm->dim; 613 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys /*, mz = da->Ze - da->Zs*/; 614 const PetscInt il = s.i - da->Xs/da->w, jl = dim > 1 ? s.j - da->Ys : 0, kl = dim > 2 ? s.k - da->Zs : 0; 615 616 PetscFunctionBegin; 617 *cell = -1; 618 if ((s.i < da->Xs/da->w) || (s.i >= da->Xe/da->w)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Stencil i %d should be in [%d, %d)", s.i, da->Xs/da->w, da->Xe/da->w); 619 if ((dim > 1) && ((s.j < da->Ys) || (s.j >= da->Ye))) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Stencil j %d should be in [%d, %d)", s.j, da->Ys, da->Ye); 620 if ((dim > 2) && ((s.k < da->Zs) || (s.k >= da->Ze))) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Stencil k %d should be in [%d, %d)", s.k, da->Zs, da->Ze); 621 *cell = (kl*my + jl)*mx + il; 622 PetscFunctionReturn(0); 623 } 624 625 PetscErrorCode DMDAComputeCellGeometry_2D(DM dm, const PetscScalar vertices[], const PetscReal refPoint[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 626 { 627 const PetscScalar x0 = vertices[0]; 628 const PetscScalar y0 = vertices[1]; 629 const PetscScalar x1 = vertices[2]; 630 const PetscScalar y1 = vertices[3]; 631 const PetscScalar x2 = vertices[4]; 632 const PetscScalar y2 = vertices[5]; 633 const PetscScalar x3 = vertices[6]; 634 const PetscScalar y3 = vertices[7]; 635 const PetscScalar f_01 = x2 - x1 - x3 + x0; 636 const PetscScalar g_01 = y2 - y1 - y3 + y0; 637 const PetscScalar x = refPoint[0]; 638 const PetscScalar y = refPoint[1]; 639 PetscReal invDet; 640 PetscErrorCode ierr; 641 642 PetscFunctionBegin; 643 #if 0 644 ierr = PetscPrintf(PETSC_COMM_SELF, "Cell (%g,%g)--(%g,%g)--(%g,%g)--(%g,%g)\n", 645 PetscRealPart(x0),PetscRealPart(y0),PetscRealPart(x1),PetscRealPart(y1),PetscRealPart(x2),PetscRealPart(y2),PetscRealPart(x3),PetscRealPart(y3));CHKERRQ(ierr); 646 ierr = PetscPrintf(PETSC_COMM_SELF, "Ref Point (%g,%g)\n", PetscRealPart(x), PetscRealPart(y));CHKERRQ(ierr); 647 #endif 648 J[0] = PetscRealPart(x1 - x0 + f_01*y) * 0.5; J[1] = PetscRealPart(x3 - x0 + f_01*x) * 0.5; 649 J[2] = PetscRealPart(y1 - y0 + g_01*y) * 0.5; J[3] = PetscRealPart(y3 - y0 + g_01*x) * 0.5; 650 *detJ = J[0]*J[3] - J[1]*J[2]; 651 invDet = 1.0/(*detJ); 652 if (invJ) { 653 invJ[0] = invDet*J[3]; invJ[1] = -invDet*J[1]; 654 invJ[2] = -invDet*J[2]; invJ[3] = invDet*J[0]; 655 } 656 ierr = PetscLogFlops(30);CHKERRQ(ierr); 657 PetscFunctionReturn(0); 658 } 659 660 PetscErrorCode DMDAComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal detJ[]) 661 { 662 DM cdm; 663 Vec coordinates; 664 const PetscReal *quadPoints; 665 PetscScalar *vertices = NULL; 666 PetscInt Nq, csize, dim, d, q; 667 PetscErrorCode ierr; 668 669 PetscFunctionBegin; 670 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 671 ierr = DMDAGetInfo(dm, &dim, 0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 672 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 673 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 674 ierr = DMDAVecGetClosure(cdm, NULL, coordinates, cell, &csize, &vertices);CHKERRQ(ierr); 675 for (d = 0; d < dim; ++d) v0[d] = PetscRealPart(vertices[d]); 676 switch (dim) { 677 case 2: 678 ierr = PetscQuadratureGetData(quad, NULL, NULL, &Nq, &quadPoints, NULL);CHKERRQ(ierr); 679 for (q = 0; q < Nq; ++q) { 680 ierr = DMDAComputeCellGeometry_2D(dm, vertices, &quadPoints[q*dim], J, invJ, detJ);CHKERRQ(ierr); 681 } 682 break; 683 default: 684 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Dimension %d not supported", dim); 685 } 686 ierr = DMDAVecRestoreClosure(cdm, NULL, coordinates, cell, &csize, &vertices);CHKERRQ(ierr); 687 PetscFunctionReturn(0); 688 } 689 690 PetscErrorCode private_DMDALocatePointsIS_2D_Regular(DM dmregular,Vec pos,IS *iscell) 691 { 692 PetscInt n,bs,p,npoints; 693 PetscInt xs,xe,Xs,Xe,mxlocal; 694 PetscInt ys,ye,Ys,Ye,mylocal; 695 PetscInt d,c0,c1; 696 PetscReal gmin_l[2],gmax_l[2],dx[2]; 697 PetscReal gmin[2],gmax[2]; 698 PetscInt *cellidx; 699 Vec coor; 700 const PetscScalar *_coor; 701 PetscErrorCode ierr; 702 703 PetscFunctionBegin; 704 ierr = DMDAGetCorners(dmregular,&xs,&ys,0,&xe,&ye,0);CHKERRQ(ierr); 705 ierr = DMDAGetGhostCorners(dmregular,&Xs,&Ys,0,&Xe,&Ye,0);CHKERRQ(ierr); 706 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 707 ye += ys; Ye += Ys; if (ys != Ys) ys -= 1; 708 709 mxlocal = xe - xs - 1; 710 mylocal = ye - ys - 1; 711 712 ierr = DMGetCoordinatesLocal(dmregular,&coor);CHKERRQ(ierr); 713 ierr = VecGetArrayRead(coor,&_coor);CHKERRQ(ierr); 714 c0 = (xs-Xs) + (ys-Ys)*(Xe-Xs); 715 c1 = (xe-2-Xs+1) + (ye-2-Ys+1)*(Xe-Xs); 716 717 gmin_l[0] = PetscRealPart(_coor[2*c0+0]); 718 gmin_l[1] = PetscRealPart(_coor[2*c0+1]); 719 720 gmax_l[0] = PetscRealPart(_coor[2*c1+0]); 721 gmax_l[1] = PetscRealPart(_coor[2*c1+1]); 722 723 dx[0] = (gmax_l[0]-gmin_l[0])/((PetscReal)mxlocal); 724 dx[1] = (gmax_l[1]-gmin_l[1])/((PetscReal)mylocal); 725 726 ierr = VecRestoreArrayRead(coor,&_coor);CHKERRQ(ierr); 727 728 ierr = DMDAGetBoundingBox(dmregular,gmin,gmax);CHKERRQ(ierr); 729 730 ierr = VecGetLocalSize(pos,&n);CHKERRQ(ierr); 731 ierr = VecGetBlockSize(pos,&bs);CHKERRQ(ierr); 732 npoints = n/bs; 733 734 ierr = PetscMalloc1(npoints,&cellidx);CHKERRQ(ierr); 735 ierr = VecGetArrayRead(pos,&_coor);CHKERRQ(ierr); 736 for (p=0; p<npoints; p++) { 737 PetscReal coor_p[2]; 738 PetscInt mi[2]; 739 740 coor_p[0] = PetscRealPart(_coor[2*p]); 741 coor_p[1] = PetscRealPart(_coor[2*p+1]); 742 743 cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND; 744 745 if (coor_p[0] < gmin_l[0]) { continue; } 746 if (coor_p[0] > gmax_l[0]) { continue; } 747 if (coor_p[1] < gmin_l[1]) { continue; } 748 if (coor_p[1] > gmax_l[1]) { continue; } 749 750 for (d=0; d<2; d++) { 751 mi[d] = (PetscInt)( (coor_p[d] - gmin[d])/dx[d] ); 752 } 753 754 if (mi[0] < xs) { continue; } 755 if (mi[0] > (xe-1)) { continue; } 756 if (mi[1] < ys) { continue; } 757 if (mi[1] > (ye-1)) { continue; } 758 759 if (mi[0] == (xe-1)) { mi[0]--; } 760 if (mi[1] == (ye-1)) { mi[1]--; } 761 762 cellidx[p] = (mi[0]-xs) + (mi[1]-ys) * mxlocal; 763 } 764 ierr = VecRestoreArrayRead(pos,&_coor);CHKERRQ(ierr); 765 ierr = ISCreateGeneral(PETSC_COMM_SELF,npoints,cellidx,PETSC_OWN_POINTER,iscell);CHKERRQ(ierr); 766 PetscFunctionReturn(0); 767 } 768 769 PetscErrorCode private_DMDALocatePointsIS_3D_Regular(DM dmregular,Vec pos,IS *iscell) 770 { 771 PetscInt n,bs,p,npoints; 772 PetscInt xs,xe,Xs,Xe,mxlocal; 773 PetscInt ys,ye,Ys,Ye,mylocal; 774 PetscInt zs,ze,Zs,Ze,mzlocal; 775 PetscInt d,c0,c1; 776 PetscReal gmin_l[3],gmax_l[3],dx[3]; 777 PetscReal gmin[3],gmax[3]; 778 PetscInt *cellidx; 779 Vec coor; 780 const PetscScalar *_coor; 781 PetscErrorCode ierr; 782 783 PetscFunctionBegin; 784 ierr = DMDAGetCorners(dmregular,&xs,&ys,&zs,&xe,&ye,&ze);CHKERRQ(ierr); 785 ierr = DMDAGetGhostCorners(dmregular,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);CHKERRQ(ierr); 786 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 787 ye += ys; Ye += Ys; if (ys != Ys) ys -= 1; 788 ze += zs; Ze += Zs; if (zs != Zs) zs -= 1; 789 790 mxlocal = xe - xs - 1; 791 mylocal = ye - ys - 1; 792 mzlocal = ze - zs - 1; 793 794 ierr = DMGetCoordinatesLocal(dmregular,&coor);CHKERRQ(ierr); 795 ierr = VecGetArrayRead(coor,&_coor);CHKERRQ(ierr); 796 c0 = (xs-Xs) + (ys-Ys) *(Xe-Xs) + (zs-Zs) *(Xe-Xs)*(Ye-Ys); 797 c1 = (xe-2-Xs+1) + (ye-2-Ys+1)*(Xe-Xs) + (ze-2-Zs+1)*(Xe-Xs)*(Ye-Ys); 798 799 gmin_l[0] = PetscRealPart(_coor[3*c0+0]); 800 gmin_l[1] = PetscRealPart(_coor[3*c0+1]); 801 gmin_l[2] = PetscRealPart(_coor[3*c0+2]); 802 803 gmax_l[0] = PetscRealPart(_coor[3*c1+0]); 804 gmax_l[1] = PetscRealPart(_coor[3*c1+1]); 805 gmax_l[2] = PetscRealPart(_coor[3*c1+2]); 806 807 dx[0] = (gmax_l[0]-gmin_l[0])/((PetscReal)mxlocal); 808 dx[1] = (gmax_l[1]-gmin_l[1])/((PetscReal)mylocal); 809 dx[2] = (gmax_l[2]-gmin_l[2])/((PetscReal)mzlocal); 810 811 ierr = VecRestoreArrayRead(coor,&_coor);CHKERRQ(ierr); 812 813 ierr = DMDAGetBoundingBox(dmregular,gmin,gmax);CHKERRQ(ierr); 814 815 ierr = VecGetLocalSize(pos,&n);CHKERRQ(ierr); 816 ierr = VecGetBlockSize(pos,&bs);CHKERRQ(ierr); 817 npoints = n/bs; 818 819 ierr = PetscMalloc1(npoints,&cellidx);CHKERRQ(ierr); 820 ierr = VecGetArrayRead(pos,&_coor);CHKERRQ(ierr); 821 for (p=0; p<npoints; p++) { 822 PetscReal coor_p[3]; 823 PetscInt mi[3]; 824 825 coor_p[0] = PetscRealPart(_coor[3*p]); 826 coor_p[1] = PetscRealPart(_coor[3*p+1]); 827 coor_p[2] = PetscRealPart(_coor[3*p+2]); 828 829 cellidx[p] = DMLOCATEPOINT_POINT_NOT_FOUND; 830 831 if (coor_p[0] < gmin_l[0]) { continue; } 832 if (coor_p[0] > gmax_l[0]) { continue; } 833 if (coor_p[1] < gmin_l[1]) { continue; } 834 if (coor_p[1] > gmax_l[1]) { continue; } 835 if (coor_p[2] < gmin_l[2]) { continue; } 836 if (coor_p[2] > gmax_l[2]) { continue; } 837 838 for (d=0; d<3; d++) { 839 mi[d] = (PetscInt)( (coor_p[d] - gmin[d])/dx[d] ); 840 } 841 842 if (mi[0] < xs) { continue; } 843 if (mi[0] > (xe-1)) { continue; } 844 if (mi[1] < ys) { continue; } 845 if (mi[1] > (ye-1)) { continue; } 846 if (mi[2] < zs) { continue; } 847 if (mi[2] > (ze-1)) { continue; } 848 849 if (mi[0] == (xe-1)) { mi[0]--; } 850 if (mi[1] == (ye-1)) { mi[1]--; } 851 if (mi[2] == (ze-1)) { mi[2]--; } 852 853 cellidx[p] = (mi[0]-xs) + (mi[1]-ys) * mxlocal + (mi[2]-zs) * mxlocal * mylocal; 854 } 855 ierr = VecRestoreArrayRead(pos,&_coor);CHKERRQ(ierr); 856 ierr = ISCreateGeneral(PETSC_COMM_SELF,npoints,cellidx,PETSC_OWN_POINTER,iscell);CHKERRQ(ierr); 857 PetscFunctionReturn(0); 858 } 859 860 PetscErrorCode DMLocatePoints_DA_Regular(DM dm,Vec pos,DMPointLocationType ltype,PetscSF cellSF) 861 { 862 IS iscell; 863 PetscSFNode *cells; 864 PetscInt p,bs,dim,npoints,nfound; 865 const PetscInt *boxCells; 866 PetscErrorCode ierr; 867 868 PetscFunctionBegin; 869 ierr = VecGetBlockSize(pos,&dim);CHKERRQ(ierr); 870 switch (dim) { 871 case 1: 872 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Support not provided for 1D"); 873 break; 874 case 2: 875 ierr = private_DMDALocatePointsIS_2D_Regular(dm,pos,&iscell);CHKERRQ(ierr); 876 break; 877 case 3: 878 ierr = private_DMDALocatePointsIS_3D_Regular(dm,pos,&iscell);CHKERRQ(ierr); 879 break; 880 default: 881 SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupport spatial dimension"); 882 break; 883 } 884 885 ierr = VecGetLocalSize(pos,&npoints);CHKERRQ(ierr); 886 ierr = VecGetBlockSize(pos,&bs);CHKERRQ(ierr); 887 npoints = npoints / bs; 888 889 ierr = PetscMalloc1(npoints, &cells);CHKERRQ(ierr); 890 ierr = ISGetIndices(iscell, &boxCells);CHKERRQ(ierr); 891 892 for (p=0; p<npoints; p++) { 893 cells[p].rank = 0; 894 cells[p].index = DMLOCATEPOINT_POINT_NOT_FOUND; 895 896 cells[p].index = boxCells[p]; 897 } 898 ierr = ISRestoreIndices(iscell, &boxCells);CHKERRQ(ierr); 899 900 nfound = npoints; 901 ierr = PetscSFSetGraph(cellSF, npoints, nfound, NULL, PETSC_OWN_POINTER, cells, PETSC_OWN_POINTER);CHKERRQ(ierr); 902 ierr = ISDestroy(&iscell);CHKERRQ(ierr); 903 904 PetscFunctionReturn(0); 905 } 906