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