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 = cx*nxF + cy + xfStart; 377 PetscInt yf = c + yfStart; 378 PetscInt points[9]; 379 380 points[0] = p; 381 points[1] = yf; points[2] = xf+nxF; points[3] = yf+nyF; points[4] = xf; 382 points[5] = v+0; points[6] = v+1; points[7] = v+nVx+1; points[8] = v+nVx+0; 383 ierr = FillClosureArray_Static(dm, section, 9, points, vArray, csize, values);CHKERRQ(ierr); 384 } else { 385 /* 6 faces, 8 vertices 386 Bottom-left-back vertex follows same order as cells 387 Back z-face follows same order as cells 388 Bottom y-face follows same order as cells 389 Left x-face follows same order as cells 390 391 14-----13 392 /| /| 393 / | 2 / | 394 / 5| / | 395 10-----9 4| 396 | 11-|---12 397 |6 / | / 398 | /1 3| / 399 |/ |/ 400 7-----8 401 */ 402 PetscInt c = p - cStart; 403 PetscInt points[15]; 404 405 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; 406 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; 407 points[12] = c+vStart+nVx*nVy+1; points[13] = c+vStart+nVx*nVy+nVx+1; points[14] = c+vStart+nVx*nVy+nVx+0; 408 409 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken"); 410 ierr = FillClosureArray_Static(dm, section, 15, points, vArray, csize, values);CHKERRQ(ierr); 411 } 412 } else if ((p >= vStart) || (p < vEnd)) { 413 /* Vertex */ 414 ierr = FillClosureArray_Static(dm, section, 1, &p, vArray, csize, values);CHKERRQ(ierr); 415 } else if ((p >= fStart) || (p < fStart + nXF)) { 416 /* X Face */ 417 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 418 else if (dim == 2) { 419 /* 2 vertices: The bottom vertex has the same numbering as the face */ 420 PetscInt f = p - xfStart; 421 PetscInt points[3]; 422 423 points[0] = p; points[1] = f; points[2] = f+nVx; 424 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken"); 425 ierr = FillClosureArray_Static(dm, section, 3, points, vArray, csize, values);CHKERRQ(ierr); 426 } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 427 } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) { 428 /* Y Face */ 429 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 430 else if (dim == 2) { 431 /* 2 vertices: The left vertex has the same numbering as the face */ 432 PetscInt f = p - yfStart; 433 PetscInt points[3]; 434 435 points[0] = p; points[1] = f; points[2] = f+1; 436 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken"); 437 ierr = FillClosureArray_Static(dm, section, 3, points, vArray, csize, values);CHKERRQ(ierr); 438 } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 439 } else { 440 /* Z Face */ 441 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 442 else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D"); 443 else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 444 } 445 PetscFunctionReturn(0); 446 } 447 448 #undef __FUNCT__ 449 #define __FUNCT__ "DMDAVecGetClosure" 450 PetscErrorCode DMDAVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt p, PetscInt *csize, PetscScalar **values) 451 { 452 PetscScalar *vArray; 453 PetscErrorCode ierr; 454 455 PetscFunctionBegin; 456 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 457 PetscValidHeaderSpecific(v, VEC_CLASSID, 3); 458 if (values) PetscValidPointer(values, 6); 459 ierr = VecGetArray(v, &vArray);CHKERRQ(ierr); 460 ierr = DMDAGetClosureScalar(dm, section, p, vArray, csize, values);CHKERRQ(ierr); 461 ierr = VecRestoreArray(v, &vArray);CHKERRQ(ierr); 462 PetscFunctionReturn(0); 463 } 464 465 #undef __FUNCT__ 466 #define __FUNCT__ "DMDARestoreClosureScalar" 467 PetscErrorCode DMDARestoreClosureScalar(DM dm, PetscSection section, PetscInt p, PetscScalar *vArray, PetscInt *csize, PetscScalar **values) 468 { 469 PetscErrorCode ierr; 470 471 PetscFunctionBegin; 472 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 473 PetscValidPointer(values, 6); 474 ierr = DMRestoreWorkArray(dm, csize ? *csize : 0, PETSC_SCALAR, (void*) values);CHKERRQ(ierr); 475 PetscFunctionReturn(0); 476 } 477 478 #undef __FUNCT__ 479 #define __FUNCT__ "DMDAVecRestoreClosure" 480 PetscErrorCode DMDAVecRestoreClosure(DM dm, PetscSection section, Vec v, PetscInt p, PetscInt *csize, PetscScalar **values) 481 { 482 PetscErrorCode ierr; 483 484 PetscFunctionBegin; 485 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 486 PetscValidHeaderSpecific(v, VEC_CLASSID, 3); 487 PetscValidPointer(values, 5); 488 ierr = DMDARestoreClosureScalar(dm, section, p, NULL, csize, values);CHKERRQ(ierr); 489 PetscFunctionReturn(0); 490 } 491 492 #undef __FUNCT__ 493 #define __FUNCT__ "DMDASetClosureScalar" 494 PetscErrorCode DMDASetClosureScalar(DM dm, PetscSection section, PetscInt p,PetscScalar *vArray, const PetscScalar *values, InsertMode mode) 495 { 496 DM_DA *da = (DM_DA*) dm->data; 497 PetscInt dim = da->dim; 498 PetscInt nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF, nCx, nCy; 499 PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart; 500 PetscErrorCode ierr; 501 502 PetscFunctionBegin; 503 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 504 PetscValidScalarPointer(values, 4); 505 PetscValidPointer(values, 5); 506 if (!section) {ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr);} 507 if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection"); 508 ierr = DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);CHKERRQ(ierr); 509 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 510 ierr = DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 511 ierr = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr); 512 ierr = DMDAGetNumCells(dm, &nCx, &nCy, NULL, NULL);CHKERRQ(ierr); 513 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);CHKERRQ(ierr); 514 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 515 xfStart = fStart; xfEnd = xfStart+nXF; 516 yfStart = xfEnd; yfEnd = yfStart+nYF; 517 zfStart = yfEnd; 518 if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd); 519 if ((p >= cStart) || (p < cEnd)) { 520 /* Cell */ 521 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 522 else if (dim == 2) { 523 /* 4 faces, 4 vertices 524 Bottom-left vertex follows same order as cells 525 Bottom y-face same order as cells 526 Left x-face follows same order as cells 527 We number the quad: 528 529 8--3--7 530 | | 531 4 0 2 532 | | 533 5--1--6 534 */ 535 PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1); 536 PetscInt v = cy*nVx + cx + vStart; 537 PetscInt xf = cx*nxF + cy + xfStart; 538 PetscInt yf = c + yfStart; 539 PetscInt points[9]; 540 541 points[0] = p; 542 points[1] = yf; points[2] = xf+nxF; points[3] = yf+nyF; points[4] = xf; 543 points[5] = v+0; points[6] = v+1; points[7] = v+nVx+1; points[8] = v+nVx+0; 544 ierr = FillClosureVec_Private(dm, section, 9, points, vArray, values, mode);CHKERRQ(ierr); 545 } else { 546 /* 6 faces, 8 vertices 547 Bottom-left-back vertex follows same order as cells 548 Back z-face follows same order as cells 549 Bottom y-face follows same order as cells 550 Left x-face follows same order as cells 551 552 14-----13 553 /| /| 554 / | 2 / | 555 / 5| / | 556 10-----9 4| 557 | 11-|---12 558 |6 / | / 559 | /1 3| / 560 |/ |/ 561 7-----8 562 */ 563 PetscInt c = p - cStart; 564 PetscInt points[15]; 565 566 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; 567 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; 568 points[13] = c+vStart+nVx*nVy+nVx+1; points[14] = c+vStart+nVx*nVy+nVx+0; 569 ierr = FillClosureVec_Private(dm, section, 15, points, vArray, values, mode);CHKERRQ(ierr); 570 } 571 } else if ((p >= vStart) || (p < vEnd)) { 572 /* Vertex */ 573 ierr = FillClosureVec_Private(dm, section, 1, &p, vArray, values, mode);CHKERRQ(ierr); 574 } else if ((p >= fStart) || (p < fStart + nXF)) { 575 /* X Face */ 576 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 577 else if (dim == 2) { 578 /* 2 vertices: The bottom vertex has the same numbering as the face */ 579 PetscInt f = p - xfStart; 580 PetscInt points[3]; 581 582 points[0] = p; points[1] = f; points[2] = f+nVx; 583 ierr = FillClosureVec_Private(dm, section, 3, points, vArray, values, mode);CHKERRQ(ierr); 584 } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 585 } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) { 586 /* Y Face */ 587 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 588 else if (dim == 2) { 589 /* 2 vertices: The left vertex has the same numbering as the face */ 590 PetscInt f = p - yfStart; 591 PetscInt points[3]; 592 593 points[0] = p; points[1] = f; points[2] = f+1; 594 ierr = FillClosureVec_Private(dm, section, 3, points, vArray, values, mode);CHKERRQ(ierr); 595 } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 596 } else { 597 /* Z Face */ 598 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 599 else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D"); 600 else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 601 } 602 PetscFunctionReturn(0); 603 } 604 605 #undef __FUNCT__ 606 #define __FUNCT__ "DMDAVecSetClosure" 607 PetscErrorCode DMDAVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt p, const PetscScalar *values, InsertMode mode) 608 { 609 PetscScalar *vArray; 610 PetscErrorCode ierr; 611 612 PetscFunctionBegin; 613 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 614 PetscValidHeaderSpecific(v, VEC_CLASSID, 3); 615 PetscValidPointer(values, 5); 616 ierr = VecGetArray(v,&vArray);CHKERRQ(ierr); 617 ierr = DMDASetClosureScalar(dm,section,p,vArray,values,mode);CHKERRQ(ierr); 618 ierr = VecRestoreArray(v,&vArray);CHKERRQ(ierr); 619 PetscFunctionReturn(0); 620 } 621 622 #undef __FUNCT__ 623 #define __FUNCT__ "DMDACnvertToCell" 624 /*@ 625 DMDAConvertToCell - Convert (i,j,k) to local cell number 626 627 Not Collective 628 629 Input Parameter: 630 + da - the distributed array 631 = s - A MatStencil giving (i,j,k) 632 633 Output Parameter: 634 . cell - the local cell number 635 636 Level: developer 637 638 .seealso: DMDAVecGetClosure() 639 @*/ 640 PetscErrorCode DMDAConvertToCell(DM dm, MatStencil s, PetscInt *cell) 641 { 642 DM_DA *da = (DM_DA*) dm->data; 643 const PetscInt dim = da->dim; 644 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys /*, mz = da->Ze - da->Zs*/; 645 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; 646 647 PetscFunctionBegin; 648 *cell = -1; 649 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); 650 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); 651 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); 652 *cell = (kl*my + jl)*mx + il; 653 PetscFunctionReturn(0); 654 } 655 656 #undef __FUNCT__ 657 #define __FUNCT__ "DMDAComputeCellGeometry_2D" 658 PetscErrorCode DMDAComputeCellGeometry_2D(DM dm, const PetscScalar vertices[], const PetscReal refPoint[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 659 { 660 const PetscScalar x0 = vertices[0]; 661 const PetscScalar y0 = vertices[1]; 662 const PetscScalar x1 = vertices[2]; 663 const PetscScalar y1 = vertices[3]; 664 const PetscScalar x2 = vertices[4]; 665 const PetscScalar y2 = vertices[5]; 666 const PetscScalar x3 = vertices[6]; 667 const PetscScalar y3 = vertices[7]; 668 const PetscScalar f_01 = x2 - x1 - x3 + x0; 669 const PetscScalar g_01 = y2 - y1 - y3 + y0; 670 const PetscScalar x = refPoint[0]; 671 const PetscScalar y = refPoint[1]; 672 PetscReal invDet; 673 PetscErrorCode ierr; 674 675 PetscFunctionBegin; 676 #if 0 677 ierr = PetscPrintf(PETSC_COMM_SELF, "Cell (%g,%g)--(%g,%g)--(%g,%g)--(%g,%g)\n", 678 PetscRealPart(x0),PetscRealPart(y0),PetscRealPart(x1),PetscRealPart(y1),PetscRealPart(x2),PetscRealPart(y2),PetscRealPart(x3),PetscRealPart(y3));CHKERRQ(ierr); 679 ierr = PetscPrintf(PETSC_COMM_SELF, "Ref Point (%g,%g)\n", PetscRealPart(x), PetscRealPart(y));CHKERRQ(ierr); 680 #endif 681 J[0] = PetscRealPart(x1 - x0 + f_01*y) * 0.5; J[1] = PetscRealPart(x3 - x0 + f_01*x) * 0.5; 682 J[2] = PetscRealPart(y1 - y0 + g_01*y) * 0.5; J[3] = PetscRealPart(y3 - y0 + g_01*x) * 0.5; 683 *detJ = J[0]*J[3] - J[1]*J[2]; 684 invDet = 1.0/(*detJ); 685 if (invJ) { 686 invJ[0] = invDet*J[3]; invJ[1] = -invDet*J[1]; 687 invJ[2] = -invDet*J[2]; invJ[3] = invDet*J[0]; 688 } 689 ierr = PetscLogFlops(30);CHKERRQ(ierr); 690 PetscFunctionReturn(0); 691 } 692 693 #undef __FUNCT__ 694 #define __FUNCT__ "DMDAComputeCellGeometry" 695 PetscErrorCode DMDAComputeCellGeometry(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal detJ[]) 696 { 697 DM cdm; 698 Vec coordinates; 699 const PetscReal *quadPoints; 700 PetscScalar *vertices = NULL; 701 PetscInt numQuadPoints, csize, dim, d, q; 702 PetscErrorCode ierr; 703 704 PetscFunctionBegin; 705 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 706 ierr = DMDAGetInfo(dm, &dim, 0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 707 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 708 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 709 ierr = DMDAVecGetClosure(cdm, NULL, coordinates, cell, &csize, &vertices);CHKERRQ(ierr); 710 for (d = 0; d < dim; ++d) v0[d] = PetscRealPart(vertices[d]); 711 switch (dim) { 712 case 2: 713 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, NULL);CHKERRQ(ierr); 714 for (q = 0; q < numQuadPoints; ++q) { 715 ierr = DMDAComputeCellGeometry_2D(dm, vertices, &quadPoints[q*dim], J, invJ, detJ);CHKERRQ(ierr); 716 } 717 break; 718 default: 719 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Dimension %d not supported", dim); 720 } 721 ierr = DMDAVecRestoreClosure(cdm, NULL, coordinates, cell, &csize, &vertices);CHKERRQ(ierr); 722 PetscFunctionReturn(0); 723 } 724