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