1 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 2 3 PETSC_STATIC_INLINE PetscErrorCode FillClosureArray_Static(DM dm, PetscSection section, PetscInt nP, const PetscInt points[], PetscScalar *vArray, PetscInt *csize, PetscScalar **array) 4 { 5 PetscScalar *a; 6 PetscInt pStart, pEnd, size = 0, dof, off, d, k, i; 7 PetscErrorCode ierr; 8 9 PetscFunctionBegin; 10 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr); 11 for (i = 0; i < nP; ++i) { 12 const PetscInt p = points[i]; 13 14 if ((p < pStart) || (p >= pEnd)) continue; 15 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 16 size += dof; 17 } 18 if (csize) *csize = size; 19 if (array) { 20 ierr = DMGetWorkArray(dm, size, PETSC_SCALAR, &a);CHKERRQ(ierr); 21 for (i = 0, k = 0; i < nP; ++i) { 22 const PetscInt p = points[i]; 23 24 if ((p < pStart) || (p >= pEnd)) continue; 25 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 26 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 27 28 for (d = 0; d < dof; ++d, ++k) a[k] = vArray[off+d]; 29 } 30 *array = a; 31 } 32 PetscFunctionReturn(0); 33 } 34 35 PETSC_STATIC_INLINE PetscErrorCode FillClosureVec_Private(DM dm, PetscSection section, PetscInt nP, const PetscInt points[], PetscScalar *vArray, const PetscScalar *array, InsertMode mode) 36 { 37 PetscInt dof, off, d, k, i; 38 PetscErrorCode ierr; 39 40 PetscFunctionBegin; 41 if ((mode == INSERT_VALUES) || (mode == INSERT_ALL_VALUES)) { 42 for (i = 0, k = 0; i < nP; ++i) { 43 ierr = PetscSectionGetDof(section, points[i], &dof);CHKERRQ(ierr); 44 ierr = PetscSectionGetOffset(section, points[i], &off);CHKERRQ(ierr); 45 46 for (d = 0; d < dof; ++d, ++k) vArray[off+d] = array[k]; 47 } 48 } else { 49 for (i = 0, k = 0; i < nP; ++i) { 50 ierr = PetscSectionGetDof(section, points[i], &dof);CHKERRQ(ierr); 51 ierr = PetscSectionGetOffset(section, points[i], &off);CHKERRQ(ierr); 52 53 for (d = 0; d < dof; ++d, ++k) vArray[off+d] += array[k]; 54 } 55 } 56 PetscFunctionReturn(0); 57 } 58 59 PETSC_STATIC_INLINE PetscErrorCode GetPointArray_Private(DM dm,PetscInt n,PetscInt *points,PetscInt *rn,const PetscInt **rpoints) 60 { 61 PetscErrorCode ierr; 62 PetscInt *work; 63 64 PetscFunctionBegin; 65 if (rn) *rn = n; 66 if (rpoints) { 67 ierr = DMGetWorkArray(dm,n,PETSC_INT,&work);CHKERRQ(ierr); 68 ierr = PetscMemcpy(work,points,n*sizeof(PetscInt));CHKERRQ(ierr); 69 *rpoints = work; 70 } 71 PetscFunctionReturn(0); 72 } 73 74 PETSC_STATIC_INLINE PetscErrorCode RestorePointArray_Private(DM dm,PetscInt *rn,const PetscInt **rpoints) 75 { 76 PetscErrorCode ierr; 77 78 PetscFunctionBegin; 79 if (rn) *rn = 0; 80 if (rpoints) { 81 /* note the second argument to DMRestoreWorkArray() is always ignored */ 82 ierr = DMRestoreWorkArray(dm,0, PETSC_INT, (void*) rpoints);CHKERRQ(ierr); 83 } 84 PetscFunctionReturn(0); 85 } 86 87 PetscErrorCode DMDAGetTransitiveClosure(DM dm, PetscInt p, PetscBool useClosure, PetscInt *closureSize, PetscInt **closure) 88 { 89 PetscInt dim = dm->dim; 90 PetscInt nVx, nVy, nVz, nxF, nXF, nyF, nYF, nzF, nZF; 91 PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart; 92 PetscErrorCode ierr; 93 94 PetscFunctionBegin; 95 if (!useClosure) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Star operation is not yet supported"); 96 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 97 PetscValidIntPointer(closureSize,4); 98 ierr = DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);CHKERRQ(ierr); 99 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 100 ierr = DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 101 ierr = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr); 102 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, NULL);CHKERRQ(ierr); 103 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 104 xfStart = fStart; xfEnd = xfStart+nXF; 105 yfStart = xfEnd; 106 if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd); 107 if ((p >= cStart) || (p < cEnd)) { 108 /* Cell */ 109 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 110 else if (dim == 2) { 111 /* 4 faces, 4 vertices 112 Bottom-left vertex follows same order as cells 113 Bottom y-face same order as cells 114 Left x-face follows same order as cells 115 We number the quad: 116 117 8--3--7 118 | | 119 4 0 2 120 | | 121 5--1--6 122 */ 123 PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1); 124 PetscInt v = cy*nVx + cx + vStart; 125 PetscInt xf = cy*nxF + cx + xfStart; 126 PetscInt yf = c + yfStart; 127 128 *closureSize = 9; 129 if (!(*closure)) {ierr = DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);CHKERRQ(ierr);} 130 (*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; 131 } else { 132 /* 6 faces, 12 edges, 8 vertices 133 Bottom-left vertex follows same order as cells 134 Bottom y-face same order as cells 135 Left x-face follows same order as cells 136 We number the quad: 137 138 8--3--7 139 | | 140 4 0 2 141 | | 142 5--1--6 143 */ 144 #if 0 145 PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1), cz = c / ((nVx-1)*(nVy-1)); 146 PetscInt v = cy*nVx + cx + vStart; 147 PetscInt xf = cy*nxF + cx + xfStart; 148 PetscInt yf = c + yfStart; 149 150 *closureSize = 26; 151 if (!(*closure)) {ierr = DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);CHKERRQ(ierr);} 152 #endif 153 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 154 } 155 } else if ((p >= vStart) || (p < vEnd)) { 156 /* Vertex */ 157 *closureSize = 1; 158 if (!(*closure)) {ierr = DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);CHKERRQ(ierr);} 159 (*closure)[0] = p; 160 } else if ((p >= fStart) || (p < fStart + nXF)) { 161 /* X Face */ 162 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 163 else if (dim == 2) { 164 /* 2 vertices: The bottom vertex has the same numbering as the face */ 165 PetscInt f = p - xfStart; 166 167 *closureSize = 3; 168 if (!(*closure)) {ierr = DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);CHKERRQ(ierr);} 169 (*closure)[0] = p; (*closure)[1] = f; (*closure)[2] = f+nVx; 170 } else if (dim == 3) { 171 /* 4 vertices */ 172 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 173 } 174 } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) { 175 /* Y Face */ 176 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 177 else if (dim == 2) { 178 /* 2 vertices: The left vertex has the same numbering as the face */ 179 PetscInt f = p - yfStart; 180 181 *closureSize = 3; 182 if (!(*closure)) {ierr = DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);CHKERRQ(ierr);} 183 (*closure)[0] = p; (*closure)[1] = f; (*closure)[2]= f+1; 184 } else if (dim == 3) { 185 /* 4 vertices */ 186 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 187 } 188 } else { 189 /* Z Face */ 190 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 191 else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D"); 192 else if (dim == 3) { 193 /* 4 vertices */ 194 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 195 } 196 } 197 PetscFunctionReturn(0); 198 } 199 200 PetscErrorCode DMDARestoreTransitiveClosure(DM dm, PetscInt p, PetscBool useClosure, PetscInt *closureSize, PetscInt **closure) 201 { 202 PetscErrorCode ierr; 203 204 PetscFunctionBegin; 205 ierr = DMRestoreWorkArray(dm, 0, PETSC_INT, closure);CHKERRQ(ierr); 206 PetscFunctionReturn(0); 207 } 208 209 PetscErrorCode DMDAGetClosure(DM dm, PetscSection section, PetscInt p,PetscInt *n,const PetscInt **closure) 210 { 211 PetscInt dim = dm->dim; 212 PetscInt nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF; 213 PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart; 214 PetscErrorCode ierr; 215 216 PetscFunctionBegin; 217 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 218 if (n) PetscValidIntPointer(n,4); 219 if (closure) PetscValidPointer(closure, 5); 220 if (!section) {ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr);} 221 if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection"); 222 ierr = DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);CHKERRQ(ierr); 223 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 224 ierr = DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 225 ierr = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr); 226 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);CHKERRQ(ierr); 227 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 228 xfStart = fStart; xfEnd = xfStart+nXF; 229 yfStart = xfEnd; 230 if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd); 231 if ((p >= cStart) || (p < cEnd)) { 232 /* Cell */ 233 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 234 else if (dim == 2) { 235 /* 4 faces, 4 vertices 236 Bottom-left vertex follows same order as cells 237 Bottom y-face same order as cells 238 Left x-face follows same order as cells 239 We number the quad: 240 241 8--3--7 242 | | 243 4 0 2 244 | | 245 5--1--6 246 */ 247 PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1); 248 PetscInt v = cy*nVx + cx + vStart; 249 PetscInt xf = cy*nxF + cx + xfStart; 250 PetscInt yf = c + yfStart; 251 PetscInt points[9]; 252 253 /* Note: initializer list is not computable at compile time, hence we fill the array manually */ 254 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; 255 256 ierr = GetPointArray_Private(dm,9,points,n,closure);CHKERRQ(ierr); 257 } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 258 } else if ((p >= vStart) || (p < vEnd)) { 259 /* Vertex */ 260 ierr = GetPointArray_Private(dm,1,&p,n,closure);CHKERRQ(ierr); 261 } else if ((p >= fStart) || (p < fStart + nXF)) { 262 /* X Face */ 263 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 264 else if (dim == 2) { 265 /* 2 vertices: The bottom vertex has the same numbering as the face */ 266 PetscInt f = p - xfStart; 267 PetscInt points[3]; 268 269 points[0] = p; points[1] = f; points[2] = f+nVx; 270 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken"); 271 ierr = GetPointArray_Private(dm,3,points,n,closure);CHKERRQ(ierr); 272 } else if (dim == 3) { 273 /* 4 vertices */ 274 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 275 } 276 } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) { 277 /* Y Face */ 278 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 279 else if (dim == 2) { 280 /* 2 vertices: The left vertex has the same numbering as the face */ 281 PetscInt f = p - yfStart; 282 PetscInt points[3]; 283 284 points[0] = p; points[1] = f; points[2]= f+1; 285 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken"); 286 ierr = GetPointArray_Private(dm, 3, points, n, closure);CHKERRQ(ierr); 287 } else if (dim == 3) { 288 /* 4 vertices */ 289 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 290 } 291 } else { 292 /* Z Face */ 293 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 294 else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D"); 295 else if (dim == 3) { 296 /* 4 vertices */ 297 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 298 } 299 } 300 PetscFunctionReturn(0); 301 } 302 303 /* Zeros n and closure. */ 304 PetscErrorCode DMDARestoreClosure(DM dm, PetscSection section, PetscInt p,PetscInt *n,const PetscInt **closure) 305 { 306 PetscErrorCode ierr; 307 308 PetscFunctionBegin; 309 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 310 if (n) PetscValidIntPointer(n,4); 311 if (closure) PetscValidPointer(closure, 5); 312 ierr = RestorePointArray_Private(dm,n,closure);CHKERRQ(ierr); 313 PetscFunctionReturn(0); 314 } 315 316 PetscErrorCode DMDAGetClosureScalar(DM dm, PetscSection section, PetscInt p, PetscScalar *vArray, PetscInt *csize, PetscScalar **values) 317 { 318 PetscInt dim = dm->dim; 319 PetscInt nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF; 320 PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart; 321 PetscErrorCode ierr; 322 323 PetscFunctionBegin; 324 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 325 PetscValidScalarPointer(vArray, 4); 326 if (values) PetscValidPointer(values, 6); 327 if (!section) {ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr);} 328 if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection"); 329 ierr = DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);CHKERRQ(ierr); 330 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 331 ierr = DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 332 ierr = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr); 333 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);CHKERRQ(ierr); 334 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 335 xfStart = fStart; xfEnd = xfStart+nXF; 336 yfStart = xfEnd; yfEnd = yfStart+nYF; 337 zfStart = yfEnd; 338 if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd); 339 if ((p >= cStart) || (p < cEnd)) { 340 /* Cell */ 341 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 342 else if (dim == 2) { 343 /* 4 faces, 4 vertices 344 Bottom-left vertex follows same order as cells 345 Bottom y-face same order as cells 346 Left x-face follows same order as cells 347 We number the quad: 348 349 8--3--7 350 | | 351 4 0 2 352 | | 353 5--1--6 354 */ 355 PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1); 356 PetscInt v = cy*nVx + cx + vStart; 357 PetscInt xf = cx*nxF + cy + xfStart; 358 PetscInt yf = c + yfStart; 359 PetscInt points[9]; 360 361 points[0] = p; 362 points[1] = yf; points[2] = xf+nxF; points[3] = yf+nyF; points[4] = xf; 363 points[5] = v+0; points[6] = v+1; points[7] = v+nVx+1; points[8] = v+nVx+0; 364 ierr = FillClosureArray_Static(dm, section, 9, points, vArray, csize, values);CHKERRQ(ierr); 365 } else { 366 /* 6 faces, 8 vertices 367 Bottom-left-back vertex follows same order as cells 368 Back z-face follows same order as cells 369 Bottom y-face follows same order as cells 370 Left x-face follows same order as cells 371 372 14-----13 373 /| /| 374 / | 2 / | 375 / 5| / | 376 10-----9 4| 377 | 11-|---12 378 |6 / | / 379 | /1 3| / 380 |/ |/ 381 7-----8 382 */ 383 PetscInt c = p - cStart; 384 PetscInt points[15]; 385 386 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; 387 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; 388 points[12] = c+vStart+nVx*nVy+1; points[13] = c+vStart+nVx*nVy+nVx+1; points[14] = c+vStart+nVx*nVy+nVx+0; 389 390 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken"); 391 ierr = FillClosureArray_Static(dm, section, 15, points, vArray, csize, values);CHKERRQ(ierr); 392 } 393 } else if ((p >= vStart) || (p < vEnd)) { 394 /* Vertex */ 395 ierr = FillClosureArray_Static(dm, section, 1, &p, vArray, csize, values);CHKERRQ(ierr); 396 } else if ((p >= fStart) || (p < fStart + nXF)) { 397 /* X Face */ 398 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 399 else if (dim == 2) { 400 /* 2 vertices: The bottom vertex has the same numbering as the face */ 401 PetscInt f = p - xfStart; 402 PetscInt points[3]; 403 404 points[0] = p; points[1] = f; points[2] = f+nVx; 405 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken"); 406 ierr = FillClosureArray_Static(dm, section, 3, points, vArray, csize, values);CHKERRQ(ierr); 407 } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 408 } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) { 409 /* Y Face */ 410 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 411 else if (dim == 2) { 412 /* 2 vertices: The left vertex has the same numbering as the face */ 413 PetscInt f = p - yfStart; 414 PetscInt points[3]; 415 416 points[0] = p; points[1] = f; points[2] = f+1; 417 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken"); 418 ierr = FillClosureArray_Static(dm, section, 3, points, vArray, csize, values);CHKERRQ(ierr); 419 } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 420 } else { 421 /* Z Face */ 422 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 423 else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D"); 424 else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 425 } 426 PetscFunctionReturn(0); 427 } 428 429 PetscErrorCode DMDAVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt p, PetscInt *csize, PetscScalar **values) 430 { 431 PetscScalar *vArray; 432 PetscErrorCode ierr; 433 434 PetscFunctionBegin; 435 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 436 PetscValidHeaderSpecific(v, VEC_CLASSID, 3); 437 if (values) PetscValidPointer(values, 6); 438 ierr = VecGetArray(v, &vArray);CHKERRQ(ierr); 439 ierr = DMDAGetClosureScalar(dm, section, p, vArray, csize, values);CHKERRQ(ierr); 440 ierr = VecRestoreArray(v, &vArray);CHKERRQ(ierr); 441 PetscFunctionReturn(0); 442 } 443 444 PetscErrorCode DMDARestoreClosureScalar(DM dm, PetscSection section, PetscInt p, PetscScalar *vArray, PetscInt *csize, PetscScalar **values) 445 { 446 PetscErrorCode ierr; 447 448 PetscFunctionBegin; 449 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 450 PetscValidPointer(values, 6); 451 ierr = DMRestoreWorkArray(dm, csize ? *csize : 0, PETSC_SCALAR, (void*) values);CHKERRQ(ierr); 452 PetscFunctionReturn(0); 453 } 454 455 PetscErrorCode DMDAVecRestoreClosure(DM dm, PetscSection section, Vec v, PetscInt p, PetscInt *csize, PetscScalar **values) 456 { 457 PetscErrorCode ierr; 458 459 PetscFunctionBegin; 460 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 461 PetscValidHeaderSpecific(v, VEC_CLASSID, 3); 462 PetscValidPointer(values, 5); 463 ierr = DMDARestoreClosureScalar(dm, section, p, NULL, csize, values);CHKERRQ(ierr); 464 PetscFunctionReturn(0); 465 } 466 467 PetscErrorCode DMDASetClosureScalar(DM dm, PetscSection section, PetscInt p,PetscScalar *vArray, const PetscScalar *values, InsertMode mode) 468 { 469 PetscInt dim = dm->dim; 470 PetscInt nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF, nCx, nCy; 471 PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart; 472 PetscErrorCode ierr; 473 474 PetscFunctionBegin; 475 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 476 PetscValidScalarPointer(values, 4); 477 PetscValidPointer(values, 5); 478 if (!section) {ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr);} 479 if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection"); 480 ierr = DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);CHKERRQ(ierr); 481 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 482 ierr = DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 483 ierr = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr); 484 ierr = DMDAGetNumCells(dm, &nCx, &nCy, NULL, NULL);CHKERRQ(ierr); 485 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);CHKERRQ(ierr); 486 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 487 xfStart = fStart; xfEnd = xfStart+nXF; 488 yfStart = xfEnd; yfEnd = yfStart+nYF; 489 zfStart = yfEnd; 490 if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd); 491 if ((p >= cStart) || (p < cEnd)) { 492 /* Cell */ 493 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 494 else if (dim == 2) { 495 /* 4 faces, 4 vertices 496 Bottom-left vertex follows same order as cells 497 Bottom y-face same order as cells 498 Left x-face follows same order as cells 499 We number the quad: 500 501 8--3--7 502 | | 503 4 0 2 504 | | 505 5--1--6 506 */ 507 PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1); 508 PetscInt v = cy*nVx + cx + vStart; 509 PetscInt xf = cx*nxF + cy + xfStart; 510 PetscInt yf = c + yfStart; 511 PetscInt points[9]; 512 513 points[0] = p; 514 points[1] = yf; points[2] = xf+nxF; points[3] = yf+nyF; points[4] = xf; 515 points[5] = v+0; points[6] = v+1; points[7] = v+nVx+1; points[8] = v+nVx+0; 516 ierr = FillClosureVec_Private(dm, section, 9, points, vArray, values, mode);CHKERRQ(ierr); 517 } else { 518 /* 6 faces, 8 vertices 519 Bottom-left-back vertex follows same order as cells 520 Back z-face follows same order as cells 521 Bottom y-face follows same order as cells 522 Left x-face follows same order as cells 523 524 14-----13 525 /| /| 526 / | 2 / | 527 / 5| / | 528 10-----9 4| 529 | 11-|---12 530 |6 / | / 531 | /1 3| / 532 |/ |/ 533 7-----8 534 */ 535 PetscInt c = p - cStart; 536 PetscInt points[15]; 537 538 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; 539 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; 540 points[13] = c+vStart+nVx*nVy+nVx+1; points[14] = c+vStart+nVx*nVy+nVx+0; 541 ierr = FillClosureVec_Private(dm, section, 15, points, vArray, values, mode);CHKERRQ(ierr); 542 } 543 } else if ((p >= vStart) || (p < vEnd)) { 544 /* Vertex */ 545 ierr = FillClosureVec_Private(dm, section, 1, &p, vArray, values, mode);CHKERRQ(ierr); 546 } else if ((p >= fStart) || (p < fStart + nXF)) { 547 /* X Face */ 548 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 549 else if (dim == 2) { 550 /* 2 vertices: The bottom vertex has the same numbering as the face */ 551 PetscInt f = p - xfStart; 552 PetscInt points[3]; 553 554 points[0] = p; points[1] = f; points[2] = f+nVx; 555 ierr = FillClosureVec_Private(dm, section, 3, points, vArray, values, mode);CHKERRQ(ierr); 556 } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 557 } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) { 558 /* Y Face */ 559 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 560 else if (dim == 2) { 561 /* 2 vertices: The left vertex has the same numbering as the face */ 562 PetscInt f = p - yfStart; 563 PetscInt points[3]; 564 565 points[0] = p; points[1] = f; points[2] = f+1; 566 ierr = FillClosureVec_Private(dm, section, 3, points, vArray, values, mode);CHKERRQ(ierr); 567 } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 568 } else { 569 /* Z Face */ 570 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 571 else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D"); 572 else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 573 } 574 PetscFunctionReturn(0); 575 } 576 577 PetscErrorCode DMDAVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt p, const PetscScalar *values, InsertMode mode) 578 { 579 PetscScalar *vArray; 580 PetscErrorCode ierr; 581 582 PetscFunctionBegin; 583 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 584 PetscValidHeaderSpecific(v, VEC_CLASSID, 3); 585 PetscValidPointer(values, 5); 586 ierr = VecGetArray(v,&vArray);CHKERRQ(ierr); 587 ierr = DMDASetClosureScalar(dm,section,p,vArray,values,mode);CHKERRQ(ierr); 588 ierr = VecRestoreArray(v,&vArray);CHKERRQ(ierr); 589 PetscFunctionReturn(0); 590 } 591 592 /*@ 593 DMDAConvertToCell - Convert (i,j,k) to local cell number 594 595 Not Collective 596 597 Input Parameter: 598 + da - the distributed array 599 - s - A MatStencil giving (i,j,k) 600 601 Output Parameter: 602 . cell - the local cell number 603 604 Level: developer 605 606 .seealso: DMDAVecGetClosure() 607 @*/ 608 PetscErrorCode DMDAConvertToCell(DM dm, MatStencil s, PetscInt *cell) 609 { 610 DM_DA *da = (DM_DA*) dm->data; 611 const PetscInt dim = dm->dim; 612 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys /*, mz = da->Ze - da->Zs*/; 613 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; 614 615 PetscFunctionBegin; 616 *cell = -1; 617 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); 618 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); 619 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); 620 *cell = (kl*my + jl)*mx + il; 621 PetscFunctionReturn(0); 622 } 623 624 PetscErrorCode DMDAComputeCellGeometry_2D(DM dm, const PetscScalar vertices[], const PetscReal refPoint[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 625 { 626 const PetscScalar x0 = vertices[0]; 627 const PetscScalar y0 = vertices[1]; 628 const PetscScalar x1 = vertices[2]; 629 const PetscScalar y1 = vertices[3]; 630 const PetscScalar x2 = vertices[4]; 631 const PetscScalar y2 = vertices[5]; 632 const PetscScalar x3 = vertices[6]; 633 const PetscScalar y3 = vertices[7]; 634 const PetscScalar f_01 = x2 - x1 - x3 + x0; 635 const PetscScalar g_01 = y2 - y1 - y3 + y0; 636 const PetscScalar x = refPoint[0]; 637 const PetscScalar y = refPoint[1]; 638 PetscReal invDet; 639 PetscErrorCode ierr; 640 641 PetscFunctionBegin; 642 #if 0 643 ierr = PetscPrintf(PETSC_COMM_SELF, "Cell (%g,%g)--(%g,%g)--(%g,%g)--(%g,%g)\n", 644 PetscRealPart(x0),PetscRealPart(y0),PetscRealPart(x1),PetscRealPart(y1),PetscRealPart(x2),PetscRealPart(y2),PetscRealPart(x3),PetscRealPart(y3));CHKERRQ(ierr); 645 ierr = PetscPrintf(PETSC_COMM_SELF, "Ref Point (%g,%g)\n", PetscRealPart(x), PetscRealPart(y));CHKERRQ(ierr); 646 #endif 647 J[0] = PetscRealPart(x1 - x0 + f_01*y) * 0.5; J[1] = PetscRealPart(x3 - x0 + f_01*x) * 0.5; 648 J[2] = PetscRealPart(y1 - y0 + g_01*y) * 0.5; J[3] = PetscRealPart(y3 - y0 + g_01*x) * 0.5; 649 *detJ = J[0]*J[3] - J[1]*J[2]; 650 invDet = 1.0/(*detJ); 651 if (invJ) { 652 invJ[0] = invDet*J[3]; invJ[1] = -invDet*J[1]; 653 invJ[2] = -invDet*J[2]; invJ[3] = invDet*J[0]; 654 } 655 ierr = PetscLogFlops(30);CHKERRQ(ierr); 656 PetscFunctionReturn(0); 657 } 658 659 PetscErrorCode DMDAComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal detJ[]) 660 { 661 DM cdm; 662 Vec coordinates; 663 const PetscReal *quadPoints; 664 PetscScalar *vertices = NULL; 665 PetscInt numQuadPoints, csize, dim, d, q; 666 PetscErrorCode ierr; 667 668 PetscFunctionBegin; 669 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 670 ierr = DMDAGetInfo(dm, &dim, 0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 671 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 672 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 673 ierr = DMDAVecGetClosure(cdm, NULL, coordinates, cell, &csize, &vertices);CHKERRQ(ierr); 674 for (d = 0; d < dim; ++d) v0[d] = PetscRealPart(vertices[d]); 675 switch (dim) { 676 case 2: 677 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, NULL);CHKERRQ(ierr); 678 for (q = 0; q < numQuadPoints; ++q) { 679 ierr = DMDAComputeCellGeometry_2D(dm, vertices, &quadPoints[q*dim], J, invJ, detJ);CHKERRQ(ierr); 680 } 681 break; 682 default: 683 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Dimension %d not supported", dim); 684 } 685 ierr = DMDAVecRestoreClosure(cdm, NULL, coordinates, cell, &csize, &vertices);CHKERRQ(ierr); 686 PetscFunctionReturn(0); 687 } 688