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