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