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