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