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 /* note the second argument to DMRestoreWorkArray() is always ignored */ 90 ierr = DMRestoreWorkArray(dm,0, PETSC_INT, (void*) rpoints);CHKERRQ(ierr); 91 } 92 PetscFunctionReturn(0); 93 } 94 95 #undef __FUNCT__ 96 #define __FUNCT__ "DMDAGetTransitiveClosure" 97 PetscErrorCode DMDAGetTransitiveClosure(DM dm, PetscInt p, PetscBool useClosure, PetscInt *closureSize, PetscInt **closure) 98 { 99 PetscInt dim = dm->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 PetscValidIntPointer(closureSize,4); 108 ierr = DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);CHKERRQ(ierr); 109 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 110 ierr = DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 111 ierr = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr); 112 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, NULL);CHKERRQ(ierr); 113 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 114 xfStart = fStart; xfEnd = xfStart+nXF; 115 yfStart = xfEnd; 116 if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd); 117 if ((p >= cStart) || (p < cEnd)) { 118 /* Cell */ 119 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 120 else if (dim == 2) { 121 /* 4 faces, 4 vertices 122 Bottom-left vertex follows same order as cells 123 Bottom y-face same order as cells 124 Left x-face follows same order as cells 125 We number the quad: 126 127 8--3--7 128 | | 129 4 0 2 130 | | 131 5--1--6 132 */ 133 PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1); 134 PetscInt v = cy*nVx + cx + vStart; 135 PetscInt xf = cy*nxF + cx + xfStart; 136 PetscInt yf = c + yfStart; 137 138 *closureSize = 9; 139 if (!(*closure)) {ierr = DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);CHKERRQ(ierr);} 140 (*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; 141 } else { 142 /* 6 faces, 12 edges, 8 vertices 143 Bottom-left vertex follows same order as cells 144 Bottom y-face same order as cells 145 Left x-face follows same order as cells 146 We number the quad: 147 148 8--3--7 149 | | 150 4 0 2 151 | | 152 5--1--6 153 */ 154 #if 0 155 PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1), cz = c / ((nVx-1)*(nVy-1)); 156 PetscInt v = cy*nVx + cx + vStart; 157 PetscInt xf = cy*nxF + cx + xfStart; 158 PetscInt yf = c + yfStart; 159 160 *closureSize = 26; 161 if (!(*closure)) {ierr = DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);CHKERRQ(ierr);} 162 #endif 163 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 164 } 165 } else if ((p >= vStart) || (p < vEnd)) { 166 /* Vertex */ 167 *closureSize = 1; 168 if (!(*closure)) {ierr = DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);CHKERRQ(ierr);} 169 (*closure)[0] = p; 170 } else if ((p >= fStart) || (p < fStart + nXF)) { 171 /* X Face */ 172 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 173 else if (dim == 2) { 174 /* 2 vertices: The bottom vertex has the same numbering as the face */ 175 PetscInt f = p - xfStart; 176 177 *closureSize = 3; 178 if (!(*closure)) {ierr = DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);CHKERRQ(ierr);} 179 (*closure)[0] = p; (*closure)[1] = f; (*closure)[2] = f+nVx; 180 } else if (dim == 3) { 181 /* 4 vertices */ 182 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 183 } 184 } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) { 185 /* Y Face */ 186 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 187 else if (dim == 2) { 188 /* 2 vertices: The left vertex has the same numbering as the face */ 189 PetscInt f = p - yfStart; 190 191 *closureSize = 3; 192 if (!(*closure)) {ierr = DMGetWorkArray(dm, *closureSize, PETSC_INT, closure);CHKERRQ(ierr);} 193 (*closure)[0] = p; (*closure)[1] = f; (*closure)[2]= f+1; 194 } else if (dim == 3) { 195 /* 4 vertices */ 196 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 197 } 198 } else { 199 /* Z Face */ 200 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 201 else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D"); 202 else if (dim == 3) { 203 /* 4 vertices */ 204 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 205 } 206 } 207 PetscFunctionReturn(0); 208 } 209 210 #undef __FUNCT__ 211 #define __FUNCT__ "DMDARestoreTransitiveClosure" 212 PetscErrorCode DMDARestoreTransitiveClosure(DM dm, PetscInt p, PetscBool useClosure, PetscInt *closureSize, PetscInt **closure) 213 { 214 PetscErrorCode ierr; 215 216 PetscFunctionBegin; 217 ierr = DMRestoreWorkArray(dm, 0, PETSC_INT, closure);CHKERRQ(ierr); 218 PetscFunctionReturn(0); 219 } 220 221 #undef __FUNCT__ 222 #define __FUNCT__ "DMDAGetClosure" 223 PetscErrorCode DMDAGetClosure(DM dm, PetscSection section, PetscInt p,PetscInt *n,const PetscInt **closure) 224 { 225 PetscInt dim = dm->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 PetscInt dim = dm->dim; 337 PetscInt nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF; 338 PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart; 339 PetscErrorCode ierr; 340 341 PetscFunctionBegin; 342 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 343 PetscValidScalarPointer(vArray, 4); 344 if (values) PetscValidPointer(values, 6); 345 if (!section) {ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr);} 346 if (!section) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection"); 347 ierr = DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);CHKERRQ(ierr); 348 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 349 ierr = DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 350 ierr = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr); 351 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, NULL, NULL);CHKERRQ(ierr); 352 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 353 xfStart = fStart; xfEnd = xfStart+nXF; 354 yfStart = xfEnd; yfEnd = yfStart+nYF; 355 zfStart = yfEnd; 356 if ((p < pStart) || (p >= pEnd)) SETERRQ3(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd); 357 if ((p >= cStart) || (p < cEnd)) { 358 /* Cell */ 359 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 360 else if (dim == 2) { 361 /* 4 faces, 4 vertices 362 Bottom-left vertex follows same order as cells 363 Bottom y-face same order as cells 364 Left x-face follows same order as cells 365 We number the quad: 366 367 8--3--7 368 | | 369 4 0 2 370 | | 371 5--1--6 372 */ 373 PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1); 374 PetscInt v = cy*nVx + cx + vStart; 375 PetscInt xf = cx*nxF + cy + xfStart; 376 PetscInt yf = c + yfStart; 377 PetscInt points[9]; 378 379 points[0] = p; 380 points[1] = yf; points[2] = xf+nxF; points[3] = yf+nyF; points[4] = xf; 381 points[5] = v+0; points[6] = v+1; points[7] = v+nVx+1; points[8] = v+nVx+0; 382 ierr = FillClosureArray_Static(dm, section, 9, points, vArray, csize, values);CHKERRQ(ierr); 383 } else { 384 /* 6 faces, 8 vertices 385 Bottom-left-back vertex follows same order as cells 386 Back z-face follows same order as cells 387 Bottom y-face follows same order as cells 388 Left x-face follows same order as cells 389 390 14-----13 391 /| /| 392 / | 2 / | 393 / 5| / | 394 10-----9 4| 395 | 11-|---12 396 |6 / | / 397 | /1 3| / 398 |/ |/ 399 7-----8 400 */ 401 PetscInt c = p - cStart; 402 PetscInt points[15]; 403 404 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; 405 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; 406 points[12] = c+vStart+nVx*nVy+1; points[13] = c+vStart+nVx*nVy+nVx+1; points[14] = c+vStart+nVx*nVy+nVx+0; 407 408 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken"); 409 ierr = FillClosureArray_Static(dm, section, 15, points, vArray, csize, values);CHKERRQ(ierr); 410 } 411 } else if ((p >= vStart) || (p < vEnd)) { 412 /* Vertex */ 413 ierr = FillClosureArray_Static(dm, section, 1, &p, vArray, csize, values);CHKERRQ(ierr); 414 } else if ((p >= fStart) || (p < fStart + nXF)) { 415 /* X Face */ 416 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 417 else if (dim == 2) { 418 /* 2 vertices: The bottom vertex has the same numbering as the face */ 419 PetscInt f = p - xfStart; 420 PetscInt points[3]; 421 422 points[0] = p; points[1] = f; points[2] = f+nVx; 423 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken"); 424 ierr = FillClosureArray_Static(dm, section, 3, points, vArray, csize, values);CHKERRQ(ierr); 425 } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 426 } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) { 427 /* Y Face */ 428 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 429 else if (dim == 2) { 430 /* 2 vertices: The left vertex has the same numbering as the face */ 431 PetscInt f = p - yfStart; 432 PetscInt points[3]; 433 434 points[0] = p; points[1] = f; points[2] = f+1; 435 SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Broken"); 436 ierr = FillClosureArray_Static(dm, section, 3, points, vArray, csize, values);CHKERRQ(ierr); 437 } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 438 } else { 439 /* Z Face */ 440 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 441 else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D"); 442 else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 443 } 444 PetscFunctionReturn(0); 445 } 446 447 #undef __FUNCT__ 448 #define __FUNCT__ "DMDAVecGetClosure" 449 PetscErrorCode DMDAVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt p, PetscInt *csize, PetscScalar **values) 450 { 451 PetscScalar *vArray; 452 PetscErrorCode ierr; 453 454 PetscFunctionBegin; 455 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 456 PetscValidHeaderSpecific(v, VEC_CLASSID, 3); 457 if (values) PetscValidPointer(values, 6); 458 ierr = VecGetArray(v, &vArray);CHKERRQ(ierr); 459 ierr = DMDAGetClosureScalar(dm, section, p, vArray, csize, values);CHKERRQ(ierr); 460 ierr = VecRestoreArray(v, &vArray);CHKERRQ(ierr); 461 PetscFunctionReturn(0); 462 } 463 464 #undef __FUNCT__ 465 #define __FUNCT__ "DMDARestoreClosureScalar" 466 PetscErrorCode DMDARestoreClosureScalar(DM dm, PetscSection section, PetscInt p, PetscScalar *vArray, PetscInt *csize, PetscScalar **values) 467 { 468 PetscErrorCode ierr; 469 470 PetscFunctionBegin; 471 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 472 PetscValidPointer(values, 6); 473 ierr = DMRestoreWorkArray(dm, csize ? *csize : 0, PETSC_SCALAR, (void*) values);CHKERRQ(ierr); 474 PetscFunctionReturn(0); 475 } 476 477 #undef __FUNCT__ 478 #define __FUNCT__ "DMDAVecRestoreClosure" 479 PetscErrorCode DMDAVecRestoreClosure(DM dm, PetscSection section, Vec v, PetscInt p, PetscInt *csize, PetscScalar **values) 480 { 481 PetscErrorCode ierr; 482 483 PetscFunctionBegin; 484 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 485 PetscValidHeaderSpecific(v, VEC_CLASSID, 3); 486 PetscValidPointer(values, 5); 487 ierr = DMDARestoreClosureScalar(dm, section, p, NULL, csize, values);CHKERRQ(ierr); 488 PetscFunctionReturn(0); 489 } 490 491 #undef __FUNCT__ 492 #define __FUNCT__ "DMDASetClosureScalar" 493 PetscErrorCode DMDASetClosureScalar(DM dm, PetscSection section, PetscInt p,PetscScalar *vArray, const PetscScalar *values, InsertMode mode) 494 { 495 PetscInt dim = dm->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, cx = c % (nVx-1), cy = c / (nVx-1); 534 PetscInt v = cy*nVx + cx + vStart; 535 PetscInt xf = cx*nxF + cy + xfStart; 536 PetscInt yf = c + yfStart; 537 PetscInt points[9]; 538 539 points[0] = p; 540 points[1] = yf; points[2] = xf+nxF; points[3] = yf+nyF; points[4] = xf; 541 points[5] = v+0; points[6] = v+1; points[7] = v+nVx+1; points[8] = v+nVx+0; 542 ierr = FillClosureVec_Private(dm, section, 9, points, vArray, values, mode);CHKERRQ(ierr); 543 } else { 544 /* 6 faces, 8 vertices 545 Bottom-left-back vertex follows same order as cells 546 Back z-face follows same order as cells 547 Bottom y-face follows same order as cells 548 Left x-face follows same order as cells 549 550 14-----13 551 /| /| 552 / | 2 / | 553 / 5| / | 554 10-----9 4| 555 | 11-|---12 556 |6 / | / 557 | /1 3| / 558 |/ |/ 559 7-----8 560 */ 561 PetscInt c = p - cStart; 562 PetscInt points[15]; 563 564 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; 565 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; 566 points[13] = c+vStart+nVx*nVy+nVx+1; points[14] = c+vStart+nVx*nVy+nVx+0; 567 ierr = FillClosureVec_Private(dm, section, 15, points, vArray, values, mode);CHKERRQ(ierr); 568 } 569 } else if ((p >= vStart) || (p < vEnd)) { 570 /* Vertex */ 571 ierr = FillClosureVec_Private(dm, section, 1, &p, vArray, values, mode);CHKERRQ(ierr); 572 } else if ((p >= fStart) || (p < fStart + nXF)) { 573 /* X Face */ 574 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 575 else if (dim == 2) { 576 /* 2 vertices: The bottom vertex has the same numbering as the face */ 577 PetscInt f = p - xfStart; 578 PetscInt points[3]; 579 580 points[0] = p; points[1] = f; points[2] = f+nVx; 581 ierr = FillClosureVec_Private(dm, section, 3, points, vArray, values, mode);CHKERRQ(ierr); 582 } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 583 } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) { 584 /* Y Face */ 585 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 586 else if (dim == 2) { 587 /* 2 vertices: The left vertex has the same numbering as the face */ 588 PetscInt f = p - yfStart; 589 PetscInt points[3]; 590 591 points[0] = p; points[1] = f; points[2] = f+1; 592 ierr = FillClosureVec_Private(dm, section, 3, points, vArray, values, mode);CHKERRQ(ierr); 593 } else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 594 } else { 595 /* Z Face */ 596 if (dim == 1) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no faces in 1D"); 597 else if (dim == 2) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "There are no z-faces in 2D"); 598 else if (dim == 3) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Not implemented"); 599 } 600 PetscFunctionReturn(0); 601 } 602 603 #undef __FUNCT__ 604 #define __FUNCT__ "DMDAVecSetClosure" 605 PetscErrorCode DMDAVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt p, const PetscScalar *values, InsertMode mode) 606 { 607 PetscScalar *vArray; 608 PetscErrorCode ierr; 609 610 PetscFunctionBegin; 611 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 612 PetscValidHeaderSpecific(v, VEC_CLASSID, 3); 613 PetscValidPointer(values, 5); 614 ierr = VecGetArray(v,&vArray);CHKERRQ(ierr); 615 ierr = DMDASetClosureScalar(dm,section,p,vArray,values,mode);CHKERRQ(ierr); 616 ierr = VecRestoreArray(v,&vArray);CHKERRQ(ierr); 617 PetscFunctionReturn(0); 618 } 619 620 #undef __FUNCT__ 621 #define __FUNCT__ "DMDAConvertToCell" 622 /*@ 623 DMDAConvertToCell - Convert (i,j,k) to local cell number 624 625 Not Collective 626 627 Input Parameter: 628 + da - the distributed array 629 - s - A MatStencil giving (i,j,k) 630 631 Output Parameter: 632 . cell - the local cell number 633 634 Level: developer 635 636 .seealso: DMDAVecGetClosure() 637 @*/ 638 PetscErrorCode DMDAConvertToCell(DM dm, MatStencil s, PetscInt *cell) 639 { 640 DM_DA *da = (DM_DA*) dm->data; 641 const PetscInt dim = dm->dim; 642 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys /*, mz = da->Ze - da->Zs*/; 643 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; 644 645 PetscFunctionBegin; 646 *cell = -1; 647 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); 648 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); 649 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); 650 *cell = (kl*my + jl)*mx + il; 651 PetscFunctionReturn(0); 652 } 653 654 #undef __FUNCT__ 655 #define __FUNCT__ "DMDAComputeCellGeometry_2D" 656 PetscErrorCode DMDAComputeCellGeometry_2D(DM dm, const PetscScalar vertices[], const PetscReal refPoint[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 657 { 658 const PetscScalar x0 = vertices[0]; 659 const PetscScalar y0 = vertices[1]; 660 const PetscScalar x1 = vertices[2]; 661 const PetscScalar y1 = vertices[3]; 662 const PetscScalar x2 = vertices[4]; 663 const PetscScalar y2 = vertices[5]; 664 const PetscScalar x3 = vertices[6]; 665 const PetscScalar y3 = vertices[7]; 666 const PetscScalar f_01 = x2 - x1 - x3 + x0; 667 const PetscScalar g_01 = y2 - y1 - y3 + y0; 668 const PetscScalar x = refPoint[0]; 669 const PetscScalar y = refPoint[1]; 670 PetscReal invDet; 671 PetscErrorCode ierr; 672 673 PetscFunctionBegin; 674 #if 0 675 ierr = PetscPrintf(PETSC_COMM_SELF, "Cell (%g,%g)--(%g,%g)--(%g,%g)--(%g,%g)\n", 676 PetscRealPart(x0),PetscRealPart(y0),PetscRealPart(x1),PetscRealPart(y1),PetscRealPart(x2),PetscRealPart(y2),PetscRealPart(x3),PetscRealPart(y3));CHKERRQ(ierr); 677 ierr = PetscPrintf(PETSC_COMM_SELF, "Ref Point (%g,%g)\n", PetscRealPart(x), PetscRealPart(y));CHKERRQ(ierr); 678 #endif 679 J[0] = PetscRealPart(x1 - x0 + f_01*y) * 0.5; J[1] = PetscRealPart(x3 - x0 + f_01*x) * 0.5; 680 J[2] = PetscRealPart(y1 - y0 + g_01*y) * 0.5; J[3] = PetscRealPart(y3 - y0 + g_01*x) * 0.5; 681 *detJ = J[0]*J[3] - J[1]*J[2]; 682 invDet = 1.0/(*detJ); 683 if (invJ) { 684 invJ[0] = invDet*J[3]; invJ[1] = -invDet*J[1]; 685 invJ[2] = -invDet*J[2]; invJ[3] = invDet*J[0]; 686 } 687 ierr = PetscLogFlops(30);CHKERRQ(ierr); 688 PetscFunctionReturn(0); 689 } 690 691 #undef __FUNCT__ 692 #define __FUNCT__ "DMDAComputeCellGeometryFEM" 693 PetscErrorCode DMDAComputeCellGeometryFEM(DM dm, PetscInt cell, PetscQuadrature quad, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal detJ[]) 694 { 695 DM cdm; 696 Vec coordinates; 697 const PetscReal *quadPoints; 698 PetscScalar *vertices = NULL; 699 PetscInt numQuadPoints, csize, dim, d, q; 700 PetscErrorCode ierr; 701 702 PetscFunctionBegin; 703 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 704 ierr = DMDAGetInfo(dm, &dim, 0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 705 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 706 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 707 ierr = DMDAVecGetClosure(cdm, NULL, coordinates, cell, &csize, &vertices);CHKERRQ(ierr); 708 for (d = 0; d < dim; ++d) v0[d] = PetscRealPart(vertices[d]); 709 switch (dim) { 710 case 2: 711 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, NULL);CHKERRQ(ierr); 712 for (q = 0; q < numQuadPoints; ++q) { 713 ierr = DMDAComputeCellGeometry_2D(dm, vertices, &quadPoints[q*dim], J, invJ, detJ);CHKERRQ(ierr); 714 } 715 break; 716 default: 717 SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Dimension %d not supported", dim); 718 } 719 ierr = DMDAVecRestoreClosure(cdm, NULL, coordinates, cell, &csize, &vertices);CHKERRQ(ierr); 720 PetscFunctionReturn(0); 721 } 722