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