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, &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__ "DMDAVecGetClosure" 61 PetscErrorCode DMDAVecGetClosure(DM dm, PetscSection section, Vec v, PetscInt p, const PetscScalar **values) 62 { 63 DM_DA *da = (DM_DA *) dm->data; 64 PetscInt dim = da->dim; 65 PetscScalar *vArray; 66 PetscInt nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF; 67 PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart; 68 PetscErrorCode ierr; 69 70 PetscFunctionBegin; 71 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 72 PetscValidHeaderSpecific(v, VEC_CLASSID, 3); 73 PetscValidPointer(values, 5); 74 if (!section) {ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr);} 75 if (!section) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection"); 76 ierr = DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);CHKERRQ(ierr); 77 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 78 ierr = DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 79 ierr = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr); 80 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 81 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 82 xfStart = fStart; xfEnd = xfStart+nXF; 83 yfStart = xfEnd; yfEnd = yfStart+nYF; 84 zfStart = yfEnd; 85 if ((p < pStart) || (p >= pEnd)) SETERRQ3(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd); 86 ierr = VecGetArray(v, &vArray);CHKERRQ(ierr); 87 if ((p >= cStart) || (p < cEnd)) { 88 /* Cell */ 89 if (dim == 1) { 90 SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Not implemented"); 91 } else if (dim == 2) { 92 /* 4 faces, 4 vertices 93 Bottom-left vertex follows same order as cells 94 Bottom y-face same order as cells 95 Left x-face follows same order as cells 96 We number the quad: 97 98 8--3--7 99 | | 100 4 0 2 101 | | 102 5--1--6 103 */ 104 PetscInt c = p - cStart, cx = c % (nVx-1), cy = c / (nVx-1); 105 PetscInt v = cy*nVx + cx + vStart; 106 PetscInt xf = cy*nxF + cx + xfStart; 107 PetscInt yf = c + yfStart; 108 PetscInt points[9] = {p, yf, xf+1, yf+nyF, xf+0, v+0, v+1, v+nVx+1, v+nVx+0}; 109 110 ierr = FillClosureArray_Private(dm, section, 9, points, vArray, values);CHKERRQ(ierr); 111 } else { 112 /* 6 faces, 8 vertices 113 Bottom-left-back vertex follows same order as cells 114 Back z-face follows same order as cells 115 Bottom y-face follows same order as cells 116 Left x-face follows same order as cells 117 118 14-----13 119 /| /| 120 / | 2 / | 121 / 5| / | 122 10-----9 4| 123 | 11-|---12 124 |6 / | / 125 | /1 3| / 126 |/ |/ 127 7-----8 128 */ 129 PetscInt c = p - cStart; 130 PetscInt points[15] = {p, c+zfStart, c+zfStart+nzF, c+yfStart, c+xfStart+nxF, c+yfStart+nyF, c+xfStart, 131 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}; 132 133 SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Broken"); 134 ierr = FillClosureArray_Private(dm, section, 15, points, vArray, values);CHKERRQ(ierr); 135 } 136 } else if ((p >= vStart) || (p < vEnd)) { 137 /* Vertex */ 138 ierr = FillClosureArray_Private(dm, section, 1, &p, vArray, values);CHKERRQ(ierr); 139 } else if ((p >= fStart) || (p < fStart + nXF)) { 140 /* X Face */ 141 if (dim == 1) { 142 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 = FillClosureArray_Private(dm, section, 3, points, vArray, values);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) { 157 SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_PLIB, "There are no faces in 1D"); 158 } else if (dim == 2) { 159 /* 2 vertices: The left vertex has the same numbering as the face */ 160 PetscInt f = p - yfStart; 161 PetscInt points[3] = {p, f, f+1}; 162 163 SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Broken"); 164 ierr = FillClosureArray_Private(dm, section, 3, points, vArray, values);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) { 172 SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_PLIB, "There are no faces in 1D"); 173 } else if (dim == 2) { 174 SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_PLIB, "There are no z-faces in 2D"); 175 } else if (dim == 3) { 176 /* 4 vertices */ 177 SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Not implemented"); 178 } 179 } 180 ierr = VecRestoreArray(v, &vArray);CHKERRQ(ierr); 181 PetscFunctionReturn(0); 182 } 183 184 #undef __FUNCT__ 185 #define __FUNCT__ "DMDAVecSetClosure" 186 PetscErrorCode DMDAVecSetClosure(DM dm, PetscSection section, Vec v, PetscInt p, const PetscScalar *values, InsertMode mode) 187 { 188 DM_DA *da = (DM_DA *) dm->data; 189 PetscInt dim = da->dim; 190 PetscScalar *vArray; 191 PetscInt nVx, nVy, nxF, nXF, nyF, nYF, nzF, nZF; 192 PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart; 193 PetscErrorCode ierr; 194 195 PetscFunctionBegin; 196 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 197 PetscValidHeaderSpecific(v, VEC_CLASSID, 3); 198 PetscValidPointer(values, 5); 199 if (!section) {ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr);} 200 if (!section) SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_ARG_WRONG, "This DM has not default PetscSection"); 201 ierr = DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);CHKERRQ(ierr); 202 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 203 ierr = DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 204 ierr = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr); 205 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 206 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 207 xfStart = fStart; xfEnd = xfStart+nXF; 208 yfStart = xfEnd; yfEnd = yfStart+nYF; 209 zfStart = yfEnd; 210 if ((p < pStart) || (p >= pEnd)) SETERRQ3(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid point %d should be in [%d, %d)", p, pStart, pEnd); 211 ierr = VecGetArray(v, &vArray);CHKERRQ(ierr); 212 if ((p >= cStart) || (p < cEnd)) { 213 /* Cell */ 214 if (dim == 1) { 215 SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Not implemented"); 216 } else if (dim == 2) { 217 /* 4 faces, 4 vertices 218 Bottom-left vertex follows same order as cells 219 Bottom y-face same order as cells 220 Left x-face follows same order as cells 221 We number the quad: 222 223 8--3--7 224 | | 225 4 0 2 226 | | 227 5--1--6 228 */ 229 PetscInt c = p - cStart; 230 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}; 231 232 ierr = FillClosureVec_Private(dm, section, 9, points, vArray, values, mode);CHKERRQ(ierr); 233 } else { 234 /* 6 faces, 8 vertices 235 Bottom-left-back vertex follows same order as cells 236 Back z-face follows same order as cells 237 Bottom y-face follows same order as cells 238 Left x-face follows same order as cells 239 240 14-----13 241 /| /| 242 / | 2 / | 243 / 5| / | 244 10-----9 4| 245 | 11-|---12 246 |6 / | / 247 | /1 3| / 248 |/ |/ 249 7-----8 250 */ 251 PetscInt c = p - cStart; 252 PetscInt points[15] = {p, c+zfStart, c+zfStart+nzF, c+yfStart, c+xfStart+nxF, c+yfStart+nyF, c+xfStart, 253 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}; 254 255 ierr = FillClosureVec_Private(dm, section, 15, points, vArray, values, mode);CHKERRQ(ierr); 256 } 257 } else if ((p >= vStart) || (p < vEnd)) { 258 /* Vertex */ 259 ierr = FillClosureVec_Private(dm, section, 1, &p, vArray, values, mode);CHKERRQ(ierr); 260 } else if ((p >= fStart) || (p < fStart + nXF)) { 261 /* X Face */ 262 if (dim == 1) { 263 SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_PLIB, "There are no faces in 1D"); 264 } else if (dim == 2) { 265 /* 2 vertices: The bottom vertex has the same numbering as the face */ 266 PetscInt f = p - xfStart; 267 PetscInt points[3] = {p, f, f+nVx}; 268 269 ierr = FillClosureVec_Private(dm, section, 3, points, vArray, values, mode);CHKERRQ(ierr); 270 } else if (dim == 3) { 271 /* 4 vertices */ 272 SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Not implemented"); 273 } 274 } else if ((p >= fStart + nXF) || (p < fStart + nXF + nYF)) { 275 /* Y Face */ 276 if (dim == 1) { 277 SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_PLIB, "There are no faces in 1D"); 278 } else if (dim == 2) { 279 /* 2 vertices: The left vertex has the same numbering as the face */ 280 PetscInt f = p - yfStart; 281 PetscInt points[3] = {p, f, f+1}; 282 283 ierr = FillClosureVec_Private(dm, section, 3, points, vArray, values, mode);CHKERRQ(ierr); 284 } else if (dim == 3) { 285 /* 4 vertices */ 286 SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Not implemented"); 287 } 288 } else { 289 /* Z Face */ 290 if (dim == 1) { 291 SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_PLIB, "There are no faces in 1D"); 292 } else if (dim == 2) { 293 SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_PLIB, "There are no z-faces in 2D"); 294 } else if (dim == 3) { 295 /* 4 vertices */ 296 SETERRQ(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Not implemented"); 297 } 298 } 299 ierr = VecRestoreArray(v, &vArray);CHKERRQ(ierr); 300 PetscFunctionReturn(0); 301 } 302 303 #undef __FUNCT__ 304 #define __FUNCT__ "DMDACnvertToCell" 305 /*@ 306 DMDAConvertToCell - Convert (i,j,k) to local cell number 307 308 Not Collective 309 310 Input Parameter: 311 + da - the distributed array 312 = s - A MatStencil giving (i,j,k) 313 314 Output Parameter: 315 . cell - the local cell number 316 317 Level: developer 318 319 .seealso: DMDAVecGetClosure() 320 @*/ 321 PetscErrorCode DMDAConvertToCell(DM dm, MatStencil s, PetscInt *cell) 322 { 323 DM_DA *da = (DM_DA *) dm->data; 324 const PetscInt dim = da->dim; 325 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 326 const PetscInt nC = (mx )*(dim > 1 ? (my )*(dim > 2 ? (mz ) : 1) : 1); 327 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; 328 329 PetscFunctionBegin; 330 *cell = -1; 331 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); 332 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); 333 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); 334 *cell = (kl*my + jl)*mx + il; 335 PetscFunctionReturn(0); 336 } 337 338 #undef __FUNCT__ 339 #define __FUNCT__ "DMDAComputeCellGeometry_2D" 340 PetscErrorCode DMDAComputeCellGeometry_2D(DM dm, const PetscScalar vertices[], const PetscReal refPoint[], PetscReal J[], PetscReal invJ[], PetscReal *detJ) 341 { 342 const PetscScalar x0 = vertices[0]; 343 const PetscScalar y0 = vertices[1]; 344 const PetscScalar x1 = vertices[2]; 345 const PetscScalar y1 = vertices[3]; 346 const PetscScalar x2 = vertices[4]; 347 const PetscScalar y2 = vertices[5]; 348 const PetscScalar x3 = vertices[6]; 349 const PetscScalar y3 = vertices[7]; 350 const PetscScalar f_01 = x2 - x1 - x3 + x0; 351 const PetscScalar g_01 = y2 - y1 - y3 + y0; 352 const PetscScalar x = refPoint[0]; 353 const PetscScalar y = refPoint[1]; 354 PetscReal invDet; 355 PetscErrorCode ierr; 356 357 PetscFunctionBegin; 358 #if defined(PETSC_USE_DEBUG) 359 ierr = PetscPrintf(PETSC_COMM_SELF, "Cell (%g,%g)--(%g,%g)--(%g,%g)--(%g,%g)\n", 360 PetscRealPart(x0),PetscRealPart(y0),PetscRealPart(x1),PetscRealPart(y1),PetscRealPart(x2),PetscRealPart(y2),PetscRealPart(x3),PetscRealPart(y3));CHKERRQ(ierr); 361 ierr = PetscPrintf(PETSC_COMM_SELF, "Ref Point (%g,%g)\n", PetscRealPart(x), PetscRealPart(y));CHKERRQ(ierr); 362 #endif 363 J[0] = PetscRealPart(x1 - x0 + f_01*y) * 0.5; J[1] = PetscRealPart(x3 - x0 + f_01*x) * 0.5; 364 J[2] = PetscRealPart(y1 - y0 + g_01*y) * 0.5; J[3] = PetscRealPart(y3 - y0 + g_01*x) * 0.5; 365 *detJ = J[0]*J[3] - J[1]*J[2]; 366 invDet = 1.0/(*detJ); 367 invJ[0] = invDet*J[3]; invJ[1] = -invDet*J[1]; 368 invJ[2] = -invDet*J[2]; invJ[3] = invDet*J[0]; 369 ierr = PetscLogFlops(30);CHKERRQ(ierr); 370 PetscFunctionReturn(0); 371 } 372 373 #undef __FUNCT__ 374 #define __FUNCT__ "DMDAComputeCellGeometry" 375 PetscErrorCode DMDAComputeCellGeometry(DM dm, PetscInt cell, PetscQuadrature *quad, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal detJ[]) 376 { 377 DM cdm; 378 Vec coordinates; 379 const PetscScalar *vertices; 380 PetscInt dim, d, q; 381 PetscErrorCode ierr; 382 383 PetscFunctionBegin; 384 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 385 ierr = DMDAGetInfo(dm, &dim, 0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 386 ierr = DMDAGetCoordinates(dm, &coordinates);CHKERRQ(ierr); 387 ierr = DMDAGetCoordinateDA(dm, &cdm);CHKERRQ(ierr); 388 ierr = DMDAVecGetClosure(cdm, PETSC_NULL, coordinates, cell, &vertices);CHKERRQ(ierr); 389 for(d = 0; d < dim; ++d) { 390 v0[d] = PetscRealPart(vertices[d]); 391 } 392 switch(dim) { 393 case 2: 394 for(q = 0; q < quad->numQuadPoints; ++q) { 395 ierr = DMDAComputeCellGeometry_2D(dm, vertices, &quad->quadPoints[q*dim], J, invJ, detJ);CHKERRQ(ierr); 396 } 397 break; 398 default: 399 SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_SUP, "Dimension %d not supported", dim); 400 } 401 PetscFunctionReturn(0); 402 } 403