1 2 /* 3 Code for manipulating distributed regular arrays in parallel. 4 */ 5 6 #include <petsc-private/dmdaimpl.h> /*I "petscdmda.h" I*/ 7 #include <petscsf.h> 8 #include <petscfe.h> 9 10 /* 11 This allows the DMDA vectors to properly tell MATLAB their dimensions 12 */ 13 #if defined(PETSC_HAVE_MATLAB_ENGINE) 14 #include <engine.h> /* MATLAB include file */ 15 #include <mex.h> /* MATLAB include file */ 16 #undef __FUNCT__ 17 #define __FUNCT__ "VecMatlabEnginePut_DA2d" 18 static PetscErrorCode VecMatlabEnginePut_DA2d(PetscObject obj,void *mengine) 19 { 20 PetscErrorCode ierr; 21 PetscInt n,m; 22 Vec vec = (Vec)obj; 23 PetscScalar *array; 24 mxArray *mat; 25 DM da; 26 27 PetscFunctionBegin; 28 ierr = VecGetDM(vec, &da);CHKERRQ(ierr); 29 if (!da) SETERRQ(PetscObjectComm((PetscObject)vec),PETSC_ERR_ARG_WRONGSTATE,"Vector not associated with a DMDA"); 30 ierr = DMDAGetGhostCorners(da,0,0,0,&m,&n,0);CHKERRQ(ierr); 31 32 ierr = VecGetArray(vec,&array);CHKERRQ(ierr); 33 #if !defined(PETSC_USE_COMPLEX) 34 mat = mxCreateDoubleMatrix(m,n,mxREAL); 35 #else 36 mat = mxCreateDoubleMatrix(m,n,mxCOMPLEX); 37 #endif 38 ierr = PetscMemcpy(mxGetPr(mat),array,n*m*sizeof(PetscScalar));CHKERRQ(ierr); 39 ierr = PetscObjectName(obj);CHKERRQ(ierr); 40 engPutVariable((Engine*)mengine,obj->name,mat); 41 42 ierr = VecRestoreArray(vec,&array);CHKERRQ(ierr); 43 PetscFunctionReturn(0); 44 } 45 #endif 46 47 48 #undef __FUNCT__ 49 #define __FUNCT__ "DMCreateLocalVector_DA" 50 PetscErrorCode DMCreateLocalVector_DA(DM da,Vec *g) 51 { 52 PetscErrorCode ierr; 53 DM_DA *dd = (DM_DA*)da->data; 54 55 PetscFunctionBegin; 56 PetscValidHeaderSpecific(da,DM_CLASSID,1); 57 PetscValidPointer(g,2); 58 if (da->defaultSection) { 59 ierr = DMCreateLocalVector_Section_Private(da,g);CHKERRQ(ierr); 60 } else { 61 ierr = VecCreate(PETSC_COMM_SELF,g);CHKERRQ(ierr); 62 ierr = VecSetSizes(*g,dd->nlocal,PETSC_DETERMINE);CHKERRQ(ierr); 63 ierr = VecSetBlockSize(*g,dd->w);CHKERRQ(ierr); 64 ierr = VecSetType(*g,da->vectype);CHKERRQ(ierr); 65 ierr = VecSetDM(*g, da);CHKERRQ(ierr); 66 #if defined(PETSC_HAVE_MATLAB_ENGINE) 67 if (dd->w == 1 && dd->dim == 2) { 68 ierr = PetscObjectComposeFunction((PetscObject)*g,"PetscMatlabEnginePut_C",VecMatlabEnginePut_DA2d);CHKERRQ(ierr); 69 } 70 #endif 71 } 72 PetscFunctionReturn(0); 73 } 74 75 #undef __FUNCT__ 76 #define __FUNCT__ "DMDAGetNumCells" 77 PetscErrorCode DMDAGetNumCells(DM dm, PetscInt *numCellsX, PetscInt *numCellsY, PetscInt *numCellsZ, PetscInt *numCells) 78 { 79 DM_DA *da = (DM_DA*) dm->data; 80 const PetscInt dim = da->dim; 81 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 82 const PetscInt nC = (mx)*(dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1); 83 84 PetscFunctionBegin; 85 if (numCellsX) { 86 PetscValidIntPointer(numCellsX,2); 87 *numCellsX = mx; 88 } 89 if (numCellsY) { 90 PetscValidIntPointer(numCellsX,3); 91 *numCellsY = my; 92 } 93 if (numCellsZ) { 94 PetscValidIntPointer(numCellsX,4); 95 *numCellsZ = mz; 96 } 97 if (numCells) { 98 PetscValidIntPointer(numCells,5); 99 *numCells = nC; 100 } 101 PetscFunctionReturn(0); 102 } 103 104 #undef __FUNCT__ 105 #define __FUNCT__ "DMDAGetNumVertices" 106 PetscErrorCode DMDAGetNumVertices(DM dm, PetscInt *numVerticesX, PetscInt *numVerticesY, PetscInt *numVerticesZ, PetscInt *numVertices) 107 { 108 DM_DA *da = (DM_DA*) dm->data; 109 const PetscInt dim = da->dim; 110 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 111 const PetscInt nVx = mx+1; 112 const PetscInt nVy = dim > 1 ? (my+1) : 1; 113 const PetscInt nVz = dim > 2 ? (mz+1) : 1; 114 const PetscInt nV = nVx*nVy*nVz; 115 116 PetscFunctionBegin; 117 if (numVerticesX) { 118 PetscValidIntPointer(numVerticesX,2); 119 *numVerticesX = nVx; 120 } 121 if (numVerticesY) { 122 PetscValidIntPointer(numVerticesY,3); 123 *numVerticesY = nVy; 124 } 125 if (numVerticesZ) { 126 PetscValidIntPointer(numVerticesZ,4); 127 *numVerticesZ = nVz; 128 } 129 if (numVertices) { 130 PetscValidIntPointer(numVertices,5); 131 *numVertices = nV; 132 } 133 PetscFunctionReturn(0); 134 } 135 136 #undef __FUNCT__ 137 #define __FUNCT__ "DMDAGetNumFaces" 138 PetscErrorCode DMDAGetNumFaces(DM dm, PetscInt *numXFacesX, PetscInt *numXFaces, PetscInt *numYFacesY, PetscInt *numYFaces, PetscInt *numZFacesZ, PetscInt *numZFaces) 139 { 140 DM_DA *da = (DM_DA*) dm->data; 141 const PetscInt dim = da->dim; 142 const PetscInt mx = (da->Xe - da->Xs)/da->w, my = da->Ye - da->Ys, mz = da->Ze - da->Zs; 143 const PetscInt nxF = (dim > 1 ? (my)*(dim > 2 ? (mz) : 1) : 1); 144 const PetscInt nXF = (mx+1)*nxF; 145 const PetscInt nyF = mx*(dim > 2 ? mz : 1); 146 const PetscInt nYF = dim > 1 ? (my+1)*nyF : 0; 147 const PetscInt nzF = mx*(dim > 1 ? my : 0); 148 const PetscInt nZF = dim > 2 ? (mz+1)*nzF : 0; 149 150 PetscFunctionBegin; 151 if (numXFacesX) { 152 PetscValidIntPointer(numXFacesX,2); 153 *numXFacesX = nxF; 154 } 155 if (numXFaces) { 156 PetscValidIntPointer(numXFaces,3); 157 *numXFaces = nXF; 158 } 159 if (numYFacesY) { 160 PetscValidIntPointer(numYFacesY,4); 161 *numYFacesY = nyF; 162 } 163 if (numYFaces) { 164 PetscValidIntPointer(numYFaces,5); 165 *numYFaces = nYF; 166 } 167 if (numZFacesZ) { 168 PetscValidIntPointer(numZFacesZ,6); 169 *numZFacesZ = nzF; 170 } 171 if (numZFaces) { 172 PetscValidIntPointer(numZFaces,7); 173 *numZFaces = nZF; 174 } 175 PetscFunctionReturn(0); 176 } 177 178 #undef __FUNCT__ 179 #define __FUNCT__ "DMDAGetHeightStratum" 180 PetscErrorCode DMDAGetHeightStratum(DM dm, PetscInt height, PetscInt *pStart, PetscInt *pEnd) 181 { 182 DM_DA *da = (DM_DA*) dm->data; 183 const PetscInt dim = da->dim; 184 PetscInt nC, nV, nXF, nYF, nZF; 185 PetscErrorCode ierr; 186 187 PetscFunctionBegin; 188 if (pStart) PetscValidIntPointer(pStart,3); 189 if (pEnd) PetscValidIntPointer(pEnd,4); 190 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 191 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr); 192 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr); 193 if (height == 0) { 194 /* Cells */ 195 if (pStart) *pStart = 0; 196 if (pEnd) *pEnd = nC; 197 } else if (height == 1) { 198 /* Faces */ 199 if (pStart) *pStart = nC+nV; 200 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 201 } else if (height == dim) { 202 /* Vertices */ 203 if (pStart) *pStart = nC; 204 if (pEnd) *pEnd = nC+nV; 205 } else if (height < 0) { 206 /* All points */ 207 if (pStart) *pStart = 0; 208 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 209 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of height %d in the DA", height); 210 PetscFunctionReturn(0); 211 } 212 213 #undef __FUNCT__ 214 #define __FUNCT__ "DMDAGetDepthStratum" 215 PetscErrorCode DMDAGetDepthStratum(DM dm, PetscInt depth, PetscInt *pStart, PetscInt *pEnd) 216 { 217 DM_DA *da = (DM_DA*) dm->data; 218 const PetscInt dim = da->dim; 219 PetscInt nC, nV, nXF, nYF, nZF; 220 PetscErrorCode ierr; 221 222 PetscFunctionBegin; 223 if (pStart) PetscValidIntPointer(pStart,3); 224 if (pEnd) PetscValidIntPointer(pEnd,4); 225 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 226 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr); 227 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr); 228 if (depth == dim) { 229 /* Cells */ 230 if (pStart) *pStart = 0; 231 if (pEnd) *pEnd = nC; 232 } else if (depth == dim-1) { 233 /* Faces */ 234 if (pStart) *pStart = nC+nV; 235 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 236 } else if (depth == 0) { 237 /* Vertices */ 238 if (pStart) *pStart = nC; 239 if (pEnd) *pEnd = nC+nV; 240 } else if (depth < 0) { 241 /* All points */ 242 if (pStart) *pStart = 0; 243 if (pEnd) *pEnd = nC+nV+nXF+nYF+nZF; 244 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "No points of depth %d in the DA", depth); 245 PetscFunctionReturn(0); 246 } 247 248 #undef __FUNCT__ 249 #define __FUNCT__ "DMDAGetConeSize" 250 PetscErrorCode DMDAGetConeSize(DM dm, PetscInt p, PetscInt *coneSize) 251 { 252 DM_DA *da = (DM_DA*) dm->data; 253 const PetscInt dim = da->dim; 254 PetscInt nC, nV, nXF, nYF, nZF; 255 PetscErrorCode ierr; 256 257 PetscFunctionBegin; 258 *coneSize = 0; 259 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 260 ierr = DMDAGetNumVertices(dm, NULL, NULL, NULL, &nV);CHKERRQ(ierr); 261 ierr = DMDAGetNumFaces(dm, NULL, &nXF, NULL, &nYF, NULL, &nZF);CHKERRQ(ierr); 262 switch (dim) { 263 case 2: 264 if (p >= 0) { 265 if (p < nC) { 266 *coneSize = 4; 267 } else if (p < nC+nV) { 268 *coneSize = 0; 269 } else if (p < nC+nV+nXF+nYF+nZF) { 270 *coneSize = 2; 271 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF); 272 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p); 273 break; 274 case 3: 275 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D"); 276 break; 277 } 278 PetscFunctionReturn(0); 279 } 280 281 #undef __FUNCT__ 282 #define __FUNCT__ "DMDAGetCone" 283 PetscErrorCode DMDAGetCone(DM dm, PetscInt p, PetscInt *cone[]) 284 { 285 DM_DA *da = (DM_DA*) dm->data; 286 const PetscInt dim = da->dim; 287 PetscInt nCx, nCy, nCz, nC, nVx, nVy, nVz, nV, nxF, nyF, nzF, nXF, nYF, nZF; 288 PetscErrorCode ierr; 289 290 PetscFunctionBegin; 291 if (!cone) {ierr = DMGetWorkArray(dm, 6, PETSC_INT, cone);CHKERRQ(ierr);} 292 ierr = DMDAGetNumCells(dm, &nCx, &nCy, &nCz, &nC);CHKERRQ(ierr); 293 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr); 294 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 295 switch (dim) { 296 case 2: 297 if (p >= 0) { 298 if (p < nC) { 299 const PetscInt cy = p / nCx; 300 const PetscInt cx = p % nCx; 301 302 (*cone)[0] = cy*nxF + cx + nC+nV; 303 (*cone)[1] = cx*nyF + cy + nyF + nC+nV+nXF; 304 (*cone)[2] = cy*nxF + cx + nxF + nC+nV; 305 (*cone)[3] = cx*nyF + cy + nC+nV+nXF; 306 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do cell cones"); 307 } else if (p < nC+nV) { 308 } else if (p < nC+nV+nXF) { 309 const PetscInt fy = (p - nC+nV) / nxF; 310 const PetscInt fx = (p - nC+nV) % nxF; 311 312 (*cone)[0] = fy*nVx + fx + nC; 313 (*cone)[1] = fy*nVx + fx + 1 + nC; 314 } else if (p < nC+nV+nXF+nYF) { 315 const PetscInt fx = (p - nC+nV+nXF) / nyF; 316 const PetscInt fy = (p - nC+nV+nXF) % nyF; 317 318 (*cone)[0] = fy*nVx + fx + nC; 319 (*cone)[1] = fy*nVx + fx + nVx + nC; 320 } else SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d should be in [0, %d)", p, nC+nV+nXF+nYF+nZF); 321 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Negative point %d is invalid", p); 322 break; 323 case 3: 324 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to do 3D"); 325 break; 326 } 327 PetscFunctionReturn(0); 328 } 329 330 #undef __FUNCT__ 331 #define __FUNCT__ "DMDARestoreCone" 332 PetscErrorCode DMDARestoreCone(DM dm, PetscInt p, PetscInt *cone[]) 333 { 334 PetscErrorCode ierr; 335 336 PetscFunctionBegin; 337 ierr = DMGetWorkArray(dm, 6, PETSC_INT, cone);CHKERRQ(ierr); 338 PetscFunctionReturn(0); 339 } 340 341 #undef __FUNCT__ 342 #define __FUNCT__ "DMDACreateSection" 343 /*@C 344 DMDACreateSection - Create a PetscSection inside the DMDA that describes data layout. This allows multiple fields with 345 different numbers of dofs on vertices, cells, and faces in each direction. 346 347 Input Parameters: 348 + dm- The DMDA 349 . numFields - The number of fields 350 . numComp - The number of components in each field 351 . numDof - The number of dofs per dimension for each field 352 . numFaceDof - The number of dofs per face for each field and direction, or NULL 353 354 Level: developer 355 356 Note: 357 The default DMDA numbering is as follows: 358 359 - Cells: [0, nC) 360 - Vertices: [nC, nC+nV) 361 - X-Faces: [nC+nV, nC+nV+nXF) normal is +- x-dir 362 - Y-Faces: [nC+nV+nXF, nC+nV+nXF+nYF) normal is +- y-dir 363 - Z-Faces: [nC+nV+nXF+nYF, nC+nV+nXF+nYF+nZF) normal is +- z-dir 364 365 We interpret the default DMDA partition as a cell partition, and the data assignment as a cell assignment. 366 @*/ 367 PetscErrorCode DMDACreateSection(DM dm, PetscInt numComp[], PetscInt numDof[], PetscInt numFaceDof[], PetscSection *s) 368 { 369 DM_DA *da = (DM_DA*) dm->data; 370 PetscSection section; 371 const PetscInt dim = da->dim; 372 PetscInt numFields, numVertexTotDof = 0, numCellTotDof = 0, numFaceTotDof[3] = {0, 0, 0}; 373 PetscSF sf; 374 PetscMPIInt rank; 375 const PetscMPIInt *neighbors; 376 PetscInt *localPoints; 377 PetscSFNode *remotePoints; 378 PetscInt nleaves = 0, nleavesCheck = 0, nL = 0; 379 PetscInt nC, nVx, nVy, nVz, nV, nxF, nXF, nyF, nYF, nzF, nZF; 380 PetscInt pStart, pEnd, cStart, cEnd, vStart, vEnd, fStart, fEnd, xfStart, xfEnd, yfStart, yfEnd, zfStart, zfEnd; 381 PetscInt f, v, c, xf, yf, zf, xn, yn, zn; 382 PetscErrorCode ierr; 383 384 PetscFunctionBegin; 385 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 386 PetscValidPointer(s, 4); 387 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 388 ierr = DMDAGetNumCells(dm, NULL, NULL, NULL, &nC);CHKERRQ(ierr); 389 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr); 390 ierr = DMDAGetNumFaces(dm, &nxF, &nXF, &nyF, &nYF, &nzF, &nZF);CHKERRQ(ierr); 391 ierr = DMDAGetHeightStratum(dm, -1, &pStart, &pEnd);CHKERRQ(ierr); 392 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 393 ierr = DMDAGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 394 ierr = DMDAGetHeightStratum(dm, dim, &vStart, &vEnd);CHKERRQ(ierr); 395 xfStart = vEnd; xfEnd = xfStart+nXF; 396 yfStart = xfEnd; yfEnd = yfStart+nYF; 397 zfStart = yfEnd; zfEnd = zfStart+nZF; 398 if (zfEnd != fEnd) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid face end %d, should be %d", zfEnd, fEnd); 399 /* Create local section */ 400 ierr = DMDAGetInfo(dm, 0,0,0,0,0,0,0, &numFields, 0,0,0,0,0);CHKERRQ(ierr); 401 for (f = 0; f < numFields; ++f) { 402 numVertexTotDof += numDof[f*(dim+1)+0]; 403 numCellTotDof += numDof[f*(dim+1)+dim]; 404 numFaceTotDof[0] += dim > 0 ? (numFaceDof ? numFaceDof[f*dim+0] : numDof[f*(dim+1)+dim-1]) : 0; 405 numFaceTotDof[1] += dim > 1 ? (numFaceDof ? numFaceDof[f*dim+1] : numDof[f*(dim+1)+dim-1]) : 0; 406 numFaceTotDof[2] += dim > 2 ? (numFaceDof ? numFaceDof[f*dim+2] : numDof[f*(dim+1)+dim-1]) : 0; 407 } 408 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ion);CHKERRQ(ierr); 409 if (numFields > 0) { 410 ierr = PetscSectionSetNumFields(section, numFields);CHKERRQ(ierr); 411 for (f = 0; f < numFields; ++f) { 412 const char *name; 413 414 ierr = DMDAGetFieldName(dm, f, &name);CHKERRQ(ierr); 415 ierr = PetscSectionSetFieldName(section, f, name ? name : "Field");CHKERRQ(ierr); 416 if (numComp) { 417 ierr = PetscSectionSetFieldComponents(section, f, numComp[f]);CHKERRQ(ierr); 418 } 419 } 420 } 421 ierr = PetscSectionSetChart(section, pStart, pEnd);CHKERRQ(ierr); 422 for (v = vStart; v < vEnd; ++v) { 423 for (f = 0; f < numFields; ++f) { 424 ierr = PetscSectionSetFieldDof(section, v, f, numDof[f*(dim+1)+0]);CHKERRQ(ierr); 425 } 426 ierr = PetscSectionSetDof(section, v, numVertexTotDof);CHKERRQ(ierr); 427 } 428 for (xf = xfStart; xf < xfEnd; ++xf) { 429 for (f = 0; f < numFields; ++f) { 430 ierr = PetscSectionSetFieldDof(section, xf, f, numFaceDof ? numFaceDof[f*dim+0] : numDof[f*(dim+1)+dim-1]);CHKERRQ(ierr); 431 } 432 ierr = PetscSectionSetDof(section, xf, numFaceTotDof[0]);CHKERRQ(ierr); 433 } 434 for (yf = yfStart; yf < yfEnd; ++yf) { 435 for (f = 0; f < numFields; ++f) { 436 ierr = PetscSectionSetFieldDof(section, yf, f, numFaceDof ? numFaceDof[f*dim+1] : numDof[f*(dim+1)+dim-1]);CHKERRQ(ierr); 437 } 438 ierr = PetscSectionSetDof(section, yf, numFaceTotDof[1]);CHKERRQ(ierr); 439 } 440 for (zf = zfStart; zf < zfEnd; ++zf) { 441 for (f = 0; f < numFields; ++f) { 442 ierr = PetscSectionSetFieldDof(section, zf, f, numFaceDof ? numFaceDof[f*dim+2] : numDof[f*(dim+1)+dim-1]);CHKERRQ(ierr); 443 } 444 ierr = PetscSectionSetDof(section, zf, numFaceTotDof[2]);CHKERRQ(ierr); 445 } 446 for (c = cStart; c < cEnd; ++c) { 447 for (f = 0; f < numFields; ++f) { 448 ierr = PetscSectionSetFieldDof(section, c, f, numDof[f*(dim+1)+dim]);CHKERRQ(ierr); 449 } 450 ierr = PetscSectionSetDof(section, c, numCellTotDof);CHKERRQ(ierr); 451 } 452 ierr = PetscSectionSetUp(section);CHKERRQ(ierr); 453 /* Create mesh point SF */ 454 ierr = DMDAGetNeighbors(dm, &neighbors);CHKERRQ(ierr); 455 for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) { 456 for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) { 457 for (xn = 0; xn < 3; ++xn) { 458 const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0; 459 const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn]; 460 461 if (neighbor >= 0 && neighbor < rank) { 462 nleaves += (!xp ? nVx : 1) * (!yp ? nVy : 1) * (!zp ? nVz : 1); /* vertices */ 463 if (xp && !yp && !zp) { 464 nleaves += nxF; /* x faces */ 465 } else if (yp && !zp && !xp) { 466 nleaves += nyF; /* y faces */ 467 } else if (zp && !xp && !yp) { 468 nleaves += nzF; /* z faces */ 469 } 470 } 471 } 472 } 473 } 474 ierr = PetscMalloc2(nleaves,&localPoints,nleaves,&remotePoints);CHKERRQ(ierr); 475 for (zn = 0; zn < (dim > 2 ? 3 : 1); ++zn) { 476 for (yn = 0; yn < (dim > 1 ? 3 : 1); ++yn) { 477 for (xn = 0; xn < 3; ++xn) { 478 const PetscInt xp = xn-1, yp = dim > 1 ? yn-1 : 0, zp = dim > 2 ? zn-1 : 0; 479 const PetscInt neighbor = neighbors[(zn*3+yn)*3+xn]; 480 PetscInt xv, yv, zv; 481 482 if (neighbor >= 0 && neighbor < rank) { 483 if (xp < 0) { /* left */ 484 if (yp < 0) { /* bottom */ 485 if (zp < 0) { /* back */ 486 const PetscInt localVertex = ( 0*nVy + 0)*nVx + 0 + nC; 487 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 488 nleavesCheck += 1; /* left bottom back vertex */ 489 490 localPoints[nL] = localVertex; 491 remotePoints[nL].rank = neighbor; 492 remotePoints[nL].index = remoteVertex; 493 ++nL; 494 } else if (zp > 0) { /* front */ 495 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; 496 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 497 nleavesCheck += 1; /* left bottom front vertex */ 498 499 localPoints[nL] = localVertex; 500 remotePoints[nL].rank = neighbor; 501 remotePoints[nL].index = remoteVertex; 502 ++nL; 503 } else { 504 nleavesCheck += nVz; /* left bottom vertices */ 505 for (zv = 0; zv < nVz; ++zv, ++nL) { 506 const PetscInt localVertex = (zv*nVy + 0)*nVx + 0 + nC; 507 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 508 509 localPoints[nL] = localVertex; 510 remotePoints[nL].rank = neighbor; 511 remotePoints[nL].index = remoteVertex; 512 } 513 } 514 } else if (yp > 0) { /* top */ 515 if (zp < 0) { /* back */ 516 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; 517 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 518 nleavesCheck += 1; /* left top back vertex */ 519 520 localPoints[nL] = localVertex; 521 remotePoints[nL].rank = neighbor; 522 remotePoints[nL].index = remoteVertex; 523 ++nL; 524 } else if (zp > 0) { /* front */ 525 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; 526 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 527 nleavesCheck += 1; /* left top front vertex */ 528 529 localPoints[nL] = localVertex; 530 remotePoints[nL].rank = neighbor; 531 remotePoints[nL].index = remoteVertex; 532 ++nL; 533 } else { 534 nleavesCheck += nVz; /* left top vertices */ 535 for (zv = 0; zv < nVz; ++zv, ++nL) { 536 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; 537 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 538 539 localPoints[nL] = localVertex; 540 remotePoints[nL].rank = neighbor; 541 remotePoints[nL].index = remoteVertex; 542 } 543 } 544 } else { 545 if (zp < 0) { /* back */ 546 nleavesCheck += nVy; /* left back vertices */ 547 for (yv = 0; yv < nVy; ++yv, ++nL) { 548 const PetscInt localVertex = ( 0*nVy + yv)*nVx + 0 + nC; 549 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 550 551 localPoints[nL] = localVertex; 552 remotePoints[nL].rank = neighbor; 553 remotePoints[nL].index = remoteVertex; 554 } 555 } else if (zp > 0) { /* front */ 556 nleavesCheck += nVy; /* left front vertices */ 557 for (yv = 0; yv < nVy; ++yv, ++nL) { 558 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; 559 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 560 561 localPoints[nL] = localVertex; 562 remotePoints[nL].rank = neighbor; 563 remotePoints[nL].index = remoteVertex; 564 } 565 } else { 566 nleavesCheck += nVy*nVz; /* left vertices */ 567 for (zv = 0; zv < nVz; ++zv) { 568 for (yv = 0; yv < nVy; ++yv, ++nL) { 569 const PetscInt localVertex = (zv*nVy + yv)*nVx + 0 + nC; 570 const PetscInt remoteVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; /* TODO: Correct this for neighbor sizes */ 571 572 localPoints[nL] = localVertex; 573 remotePoints[nL].rank = neighbor; 574 remotePoints[nL].index = remoteVertex; 575 } 576 } 577 nleavesCheck += nxF; /* left faces */ 578 for (xf = 0; xf < nxF; ++xf, ++nL) { 579 /* THIS IS WRONG */ 580 const PetscInt localFace = 0 + nC+nV; 581 const PetscInt remoteFace = 0 + nC+nV; 582 583 localPoints[nL] = localFace; 584 remotePoints[nL].rank = neighbor; 585 remotePoints[nL].index = remoteFace; 586 } 587 } 588 } 589 } else if (xp > 0) { /* right */ 590 if (yp < 0) { /* bottom */ 591 if (zp < 0) { /* back */ 592 const PetscInt localVertex = ( 0*nVy + 0)*nVx + nVx-1 + nC; 593 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 594 nleavesCheck += 1; /* right bottom back vertex */ 595 596 localPoints[nL] = localVertex; 597 remotePoints[nL].rank = neighbor; 598 remotePoints[nL].index = remoteVertex; 599 ++nL; 600 } else if (zp > 0) { /* front */ 601 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + nVx-1 + nC; 602 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 603 nleavesCheck += 1; /* right bottom front vertex */ 604 605 localPoints[nL] = localVertex; 606 remotePoints[nL].rank = neighbor; 607 remotePoints[nL].index = remoteVertex; 608 ++nL; 609 } else { 610 nleavesCheck += nVz; /* right bottom vertices */ 611 for (zv = 0; zv < nVz; ++zv, ++nL) { 612 const PetscInt localVertex = (zv*nVy + 0)*nVx + nVx-1 + nC; 613 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 614 615 localPoints[nL] = localVertex; 616 remotePoints[nL].rank = neighbor; 617 remotePoints[nL].index = remoteVertex; 618 } 619 } 620 } else if (yp > 0) { /* top */ 621 if (zp < 0) { /* back */ 622 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + nVx-1 + nC; 623 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 624 nleavesCheck += 1; /* right top back vertex */ 625 626 localPoints[nL] = localVertex; 627 remotePoints[nL].rank = neighbor; 628 remotePoints[nL].index = remoteVertex; 629 ++nL; 630 } else if (zp > 0) { /* front */ 631 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + nVx-1 + nC; 632 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 633 nleavesCheck += 1; /* right top front vertex */ 634 635 localPoints[nL] = localVertex; 636 remotePoints[nL].rank = neighbor; 637 remotePoints[nL].index = remoteVertex; 638 ++nL; 639 } else { 640 nleavesCheck += nVz; /* right top vertices */ 641 for (zv = 0; zv < nVz; ++zv, ++nL) { 642 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + nVx-1 + nC; 643 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 644 645 localPoints[nL] = localVertex; 646 remotePoints[nL].rank = neighbor; 647 remotePoints[nL].index = remoteVertex; 648 } 649 } 650 } else { 651 if (zp < 0) { /* back */ 652 nleavesCheck += nVy; /* right back vertices */ 653 for (yv = 0; yv < nVy; ++yv, ++nL) { 654 const PetscInt localVertex = ( 0*nVy + yv)*nVx + nVx-1 + nC; 655 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 656 657 localPoints[nL] = localVertex; 658 remotePoints[nL].rank = neighbor; 659 remotePoints[nL].index = remoteVertex; 660 } 661 } else if (zp > 0) { /* front */ 662 nleavesCheck += nVy; /* right front vertices */ 663 for (yv = 0; yv < nVy; ++yv, ++nL) { 664 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + nVx-1 + nC; 665 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 666 667 localPoints[nL] = localVertex; 668 remotePoints[nL].rank = neighbor; 669 remotePoints[nL].index = remoteVertex; 670 } 671 } else { 672 nleavesCheck += nVy*nVz; /* right vertices */ 673 for (zv = 0; zv < nVz; ++zv) { 674 for (yv = 0; yv < nVy; ++yv, ++nL) { 675 const PetscInt localVertex = (zv*nVy + yv)*nVx + nVx-1 + nC; 676 const PetscInt remoteVertex = (zv*nVy + yv)*nVx + 0 + nC; /* TODO: Correct this for neighbor sizes */ 677 678 localPoints[nL] = localVertex; 679 remotePoints[nL].rank = neighbor; 680 remotePoints[nL].index = remoteVertex; 681 } 682 } 683 nleavesCheck += nxF; /* right faces */ 684 for (xf = 0; xf < nxF; ++xf, ++nL) { 685 /* THIS IS WRONG */ 686 const PetscInt localFace = 0 + nC+nV; 687 const PetscInt remoteFace = 0 + nC+nV; 688 689 localPoints[nL] = localFace; 690 remotePoints[nL].rank = neighbor; 691 remotePoints[nL].index = remoteFace; 692 } 693 } 694 } 695 } else { 696 if (yp < 0) { /* bottom */ 697 if (zp < 0) { /* back */ 698 nleavesCheck += nVx; /* bottom back vertices */ 699 for (xv = 0; xv < nVx; ++xv, ++nL) { 700 const PetscInt localVertex = ( 0*nVy + 0)*nVx + xv + nC; 701 const PetscInt remoteVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 702 703 localPoints[nL] = localVertex; 704 remotePoints[nL].rank = neighbor; 705 remotePoints[nL].index = remoteVertex; 706 } 707 } else if (zp > 0) { /* front */ 708 nleavesCheck += nVx; /* bottom front vertices */ 709 for (xv = 0; xv < nVx; ++xv, ++nL) { 710 const PetscInt localVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; 711 const PetscInt remoteVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 712 713 localPoints[nL] = localVertex; 714 remotePoints[nL].rank = neighbor; 715 remotePoints[nL].index = remoteVertex; 716 } 717 } else { 718 nleavesCheck += nVx*nVz; /* bottom vertices */ 719 for (zv = 0; zv < nVz; ++zv) { 720 for (xv = 0; xv < nVx; ++xv, ++nL) { 721 const PetscInt localVertex = (zv*nVy + 0)*nVx + xv + nC; 722 const PetscInt remoteVertex = (zv*nVy + nVy-1)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 723 724 localPoints[nL] = localVertex; 725 remotePoints[nL].rank = neighbor; 726 remotePoints[nL].index = remoteVertex; 727 } 728 } 729 nleavesCheck += nyF; /* bottom faces */ 730 for (yf = 0; yf < nyF; ++yf, ++nL) { 731 /* THIS IS WRONG */ 732 const PetscInt localFace = 0 + nC+nV; 733 const PetscInt remoteFace = 0 + nC+nV; 734 735 localPoints[nL] = localFace; 736 remotePoints[nL].rank = neighbor; 737 remotePoints[nL].index = remoteFace; 738 } 739 } 740 } else if (yp > 0) { /* top */ 741 if (zp < 0) { /* back */ 742 nleavesCheck += nVx; /* top back vertices */ 743 for (xv = 0; xv < nVx; ++xv, ++nL) { 744 const PetscInt localVertex = ( 0*nVy + nVy-1)*nVx + xv + nC; 745 const PetscInt remoteVertex = ((nVz-1)*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 746 747 localPoints[nL] = localVertex; 748 remotePoints[nL].rank = neighbor; 749 remotePoints[nL].index = remoteVertex; 750 } 751 } else if (zp > 0) { /* front */ 752 nleavesCheck += nVx; /* top front vertices */ 753 for (xv = 0; xv < nVx; ++xv, ++nL) { 754 const PetscInt localVertex = ((nVz-1)*nVy + nVy-1)*nVx + xv + nC; 755 const PetscInt remoteVertex = ( 0*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 756 757 localPoints[nL] = localVertex; 758 remotePoints[nL].rank = neighbor; 759 remotePoints[nL].index = remoteVertex; 760 } 761 } else { 762 nleavesCheck += nVx*nVz; /* top vertices */ 763 for (zv = 0; zv < nVz; ++zv) { 764 for (xv = 0; xv < nVx; ++xv, ++nL) { 765 const PetscInt localVertex = (zv*nVy + nVy-1)*nVx + xv + nC; 766 const PetscInt remoteVertex = (zv*nVy + 0)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 767 768 localPoints[nL] = localVertex; 769 remotePoints[nL].rank = neighbor; 770 remotePoints[nL].index = remoteVertex; 771 } 772 } 773 nleavesCheck += nyF; /* top faces */ 774 for (yf = 0; yf < nyF; ++yf, ++nL) { 775 /* THIS IS WRONG */ 776 const PetscInt localFace = 0 + nC+nV; 777 const PetscInt remoteFace = 0 + nC+nV; 778 779 localPoints[nL] = localFace; 780 remotePoints[nL].rank = neighbor; 781 remotePoints[nL].index = remoteFace; 782 } 783 } 784 } else { 785 if (zp < 0) { /* back */ 786 nleavesCheck += nVx*nVy; /* back vertices */ 787 for (yv = 0; yv < nVy; ++yv) { 788 for (xv = 0; xv < nVx; ++xv, ++nL) { 789 const PetscInt localVertex = ( 0*nVy + yv)*nVx + xv + nC; 790 const PetscInt remoteVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 791 792 localPoints[nL] = localVertex; 793 remotePoints[nL].rank = neighbor; 794 remotePoints[nL].index = remoteVertex; 795 } 796 } 797 nleavesCheck += nzF; /* back faces */ 798 for (zf = 0; zf < nzF; ++zf, ++nL) { 799 /* THIS IS WRONG */ 800 const PetscInt localFace = 0 + nC+nV; 801 const PetscInt remoteFace = 0 + nC+nV; 802 803 localPoints[nL] = localFace; 804 remotePoints[nL].rank = neighbor; 805 remotePoints[nL].index = remoteFace; 806 } 807 } else if (zp > 0) { /* front */ 808 nleavesCheck += nVx*nVy; /* front vertices */ 809 for (yv = 0; yv < nVy; ++yv) { 810 for (xv = 0; xv < nVx; ++xv, ++nL) { 811 const PetscInt localVertex = ((nVz-1)*nVy + yv)*nVx + xv + nC; 812 const PetscInt remoteVertex = ( 0*nVy + yv)*nVx + xv + nC; /* TODO: Correct this for neighbor sizes */ 813 814 localPoints[nL] = localVertex; 815 remotePoints[nL].rank = neighbor; 816 remotePoints[nL].index = remoteVertex; 817 } 818 } 819 nleavesCheck += nzF; /* front faces */ 820 for (zf = 0; zf < nzF; ++zf, ++nL) { 821 /* THIS IS WRONG */ 822 const PetscInt localFace = 0 + nC+nV; 823 const PetscInt remoteFace = 0 + nC+nV; 824 825 localPoints[nL] = localFace; 826 remotePoints[nL].rank = neighbor; 827 remotePoints[nL].index = remoteFace; 828 } 829 } else { 830 /* Nothing is shared from the interior */ 831 } 832 } 833 } 834 } 835 } 836 } 837 } 838 /* TODO: Remove duplication in leaf determination */ 839 if (nleaves != nleavesCheck) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "The number of leaves %d did not match the number of remote leaves %d", nleaves, nleavesCheck); 840 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), &sf);CHKERRQ(ierr); 841 ierr = PetscSFSetGraph(sf, pEnd, nleaves, localPoints, PETSC_OWN_POINTER, remotePoints, PETSC_OWN_POINTER);CHKERRQ(ierr); 842 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr); 843 ierr = PetscSFDestroy(&sf);CHKERRQ(ierr); 844 *s = section; 845 PetscFunctionReturn(0); 846 } 847 848 #undef __FUNCT__ 849 #define __FUNCT__ "DMDASetVertexCoordinates" 850 PetscErrorCode DMDASetVertexCoordinates(DM dm, PetscReal xl, PetscReal xu, PetscReal yl, PetscReal yu, PetscReal zl, PetscReal zu) 851 { 852 DM_DA *da = (DM_DA *) dm->data; 853 Vec coordinates; 854 PetscSection section; 855 PetscScalar *coords; 856 PetscReal h[3]; 857 PetscInt dim, size, M, N, P, nVx, nVy, nVz, nV, vStart, vEnd, v, i, j, k; 858 PetscErrorCode ierr; 859 860 PetscFunctionBegin; 861 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 862 ierr = DMDAGetInfo(dm, &dim, &M, &N, &P, 0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 863 h[0] = (xu - xl)/M; 864 h[1] = (yu - yl)/N; 865 h[2] = (zu - zl)/P; 866 ierr = DMDAGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 867 ierr = DMDAGetNumVertices(dm, &nVx, &nVy, &nVz, &nV);CHKERRQ(ierr); 868 ierr = PetscSectionCreate(PetscObjectComm((PetscObject) dm), §ion);CHKERRQ(ierr); 869 ierr = PetscSectionSetNumFields(section, 1);CHKERRQ(ierr); 870 ierr = PetscSectionSetFieldComponents(section, 0, dim);CHKERRQ(ierr); 871 ierr = PetscSectionSetChart(section, vStart, vEnd);CHKERRQ(ierr); 872 for (v = vStart; v < vEnd; ++v) { 873 ierr = PetscSectionSetDof(section, v, dim);CHKERRQ(ierr); 874 } 875 ierr = PetscSectionSetUp(section);CHKERRQ(ierr); 876 ierr = PetscSectionGetStorageSize(section, &size);CHKERRQ(ierr); 877 ierr = VecCreateSeq(PETSC_COMM_SELF, size, &coordinates);CHKERRQ(ierr); 878 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 879 for (k = 0; k < nVz; ++k) { 880 PetscInt ind[3] = {0, 0, k + da->zs}, d, off; 881 882 for (j = 0; j < nVy; ++j) { 883 ind[1] = j + da->ys; 884 for (i = 0; i < nVx; ++i) { 885 const PetscInt vertex = (k*nVy + j)*nVx + i + vStart; 886 887 ierr = PetscSectionGetOffset(section, vertex, &off);CHKERRQ(ierr); 888 ind[0] = i + da->xs; 889 for (d = 0; d < dim; ++d) { 890 coords[off+d] = h[d]*ind[d]; 891 } 892 } 893 } 894 } 895 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 896 ierr = DMSetCoordinateSection(dm, section);CHKERRQ(ierr); 897 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr); 898 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 899 ierr = VecDestroy(&coordinates);CHKERRQ(ierr); 900 PetscFunctionReturn(0); 901 } 902 903 #undef __FUNCT__ 904 #define __FUNCT__ "DMDAProjectFunctionLocal" 905 PetscErrorCode DMDAProjectFunctionLocal(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *), InsertMode mode, Vec localX) 906 { 907 PetscDualSpace *sp; 908 PetscQuadrature q; 909 PetscSection section; 910 PetscScalar *values; 911 PetscReal *v0, *J, *detJ; 912 PetscInt numFields, numComp, dim, spDim, totDim = 0, numValues, cStart, cEnd, f, c, v, d; 913 PetscErrorCode ierr; 914 915 PetscFunctionBegin; 916 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 917 ierr = PetscFEGetQuadrature(fe[0], &q);CHKERRQ(ierr); 918 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 919 ierr = PetscMalloc(numFields * sizeof(PetscDualSpace), &sp);CHKERRQ(ierr); 920 for (f = 0; f < numFields; ++f) { 921 ierr = PetscFEGetDualSpace(fe[f], &sp[f]);CHKERRQ(ierr); 922 ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr); 923 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 924 totDim += spDim*numComp; 925 } 926 ierr = DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 927 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 928 ierr = DMDAVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 929 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 930 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 931 ierr = PetscMalloc3(dim*q.numPoints,PetscReal,&v0,dim*dim*q.numPoints,PetscReal,&J,q.numPoints,PetscReal,&detJ);CHKERRQ(ierr); 932 for (c = cStart; c < cEnd; ++c) { 933 PetscCellGeometry geom; 934 935 ierr = DMDAComputeCellGeometry(dm, c, &q, v0, J, NULL, detJ);CHKERRQ(ierr); 936 geom.v0 = v0; 937 geom.J = J; 938 geom.detJ = detJ; 939 for (f = 0, v = 0; f < numFields; ++f) { 940 ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr); 941 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 942 for (d = 0; d < spDim; ++d) { 943 ierr = PetscDualSpaceApply(sp[f], d, geom, numComp, funcs[f], &values[v]);CHKERRQ(ierr); 944 v += numComp; 945 } 946 } 947 ierr = DMDAVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr); 948 } 949 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 950 ierr = PetscFree3(v0,J,detJ);CHKERRQ(ierr); 951 ierr = PetscFree(sp);CHKERRQ(ierr); 952 PetscFunctionReturn(0); 953 } 954 955 #undef __FUNCT__ 956 #define __FUNCT__ "DMDAProjectFunction" 957 /*@C 958 DMDAProjectFunction - This projects the given function into the function space provided. 959 960 Input Parameters: 961 + dm - The DM 962 . fe - The PetscFE associated with the field 963 . funcs - The coordinate functions to evaluate 964 - mode - The insertion mode for values 965 966 Output Parameter: 967 . X - vector 968 969 Level: developer 970 971 .seealso: DMDAComputeL2Diff() 972 @*/ 973 PetscErrorCode DMDAProjectFunction(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *), InsertMode mode, Vec X) 974 { 975 Vec localX; 976 PetscErrorCode ierr; 977 978 PetscFunctionBegin; 979 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 980 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 981 ierr = DMDAProjectFunctionLocal(dm, fe, funcs, mode, localX);CHKERRQ(ierr); 982 ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr); 983 ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr); 984 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 985 PetscFunctionReturn(0); 986 } 987 988 #undef __FUNCT__ 989 #define __FUNCT__ "DMDAComputeL2Diff" 990 /*@C 991 DMDAComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h. 992 993 Input Parameters: 994 + dm - The DM 995 . fe - The PetscFE object for each field 996 . funcs - The functions to evaluate for each field component 997 - X - The coefficient vector u_h 998 999 Output Parameter: 1000 . diff - The diff ||u - u_h||_2 1001 1002 Level: developer 1003 1004 .seealso: DMDAProjectFunction() 1005 @*/ 1006 PetscErrorCode DMDAComputeL2Diff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *), Vec X, PetscReal *diff) 1007 { 1008 const PetscInt debug = 0; 1009 PetscSection section; 1010 PetscQuadrature quad; 1011 Vec localX; 1012 PetscScalar *funcVal; 1013 PetscReal *coords, *v0, *J, *invJ, detJ; 1014 PetscReal localDiff = 0.0; 1015 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp; 1016 PetscErrorCode ierr; 1017 1018 PetscFunctionBegin; 1019 ierr = DMDAGetInfo(dm, &dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 1020 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1021 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1022 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 1023 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1024 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 1025 for (field = 0; field < numFields; ++field) { 1026 PetscInt Nc; 1027 1028 ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr); 1029 numComponents += Nc; 1030 } 1031 /* There are no BC values in DAs right now: ierr = DMDAProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 1032 ierr = PetscMalloc5(numComponents,PetscScalar,&funcVal,dim,PetscReal,&coords,dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);CHKERRQ(ierr); 1033 ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1034 ierr = PetscFEGetQuadrature(fe[0], &quad);CHKERRQ(ierr); 1035 for (c = cStart; c < cEnd; ++c) { 1036 PetscScalar *x = NULL; 1037 PetscReal elemDiff = 0.0; 1038 1039 ierr = DMDAComputeCellGeometry(dm, c, &quad, v0, J, invJ, &detJ);CHKERRQ(ierr); 1040 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 1041 ierr = DMDAVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1042 1043 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 1044 const PetscInt numQuadPoints = quad.numPoints; 1045 const PetscReal *quadPoints = quad.points; 1046 const PetscReal *quadWeights = quad.weights; 1047 PetscReal *basis; 1048 PetscInt numBasisFuncs, numBasisComps, q, d, e, fc, f; 1049 1050 ierr = PetscFEGetDimension(fe[field], &numBasisFuncs);CHKERRQ(ierr); 1051 ierr = PetscFEGetNumComponents(fe[field], &numBasisComps);CHKERRQ(ierr); 1052 ierr = PetscFEGetDefaultTabulation(fe[field], &basis, NULL, NULL);CHKERRQ(ierr); 1053 if (debug) { 1054 char title[1024]; 1055 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 1056 ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr); 1057 } 1058 for (q = 0; q < numQuadPoints; ++q) { 1059 for (d = 0; d < dim; d++) { 1060 coords[d] = v0[d]; 1061 for (e = 0; e < dim; e++) { 1062 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 1063 } 1064 } 1065 (*funcs[field])(coords, funcVal); 1066 for (fc = 0; fc < numBasisComps; ++fc) { 1067 PetscScalar interpolant = 0.0; 1068 1069 for (f = 0; f < numBasisFuncs; ++f) { 1070 const PetscInt fidx = f*numBasisComps+fc; 1071 interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx]; 1072 } 1073 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);} 1074 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 1075 } 1076 } 1077 comp += numBasisComps; 1078 fieldOffset += numBasisFuncs*numBasisComps; 1079 } 1080 ierr = DMDAVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 1081 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 1082 localDiff += elemDiff; 1083 } 1084 ierr = PetscFree5(funcVal,coords,v0,J,invJ);CHKERRQ(ierr); 1085 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1086 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 1087 *diff = PetscSqrtReal(*diff); 1088 PetscFunctionReturn(0); 1089 } 1090 1091 /* ------------------------------------------------------------------- */ 1092 1093 #undef __FUNCT__ 1094 #define __FUNCT__ "DMDAGetArray" 1095 /*@C 1096 DMDAGetArray - Gets a work array for a DMDA 1097 1098 Input Parameter: 1099 + da - information about my local patch 1100 - ghosted - do you want arrays for the ghosted or nonghosted patch 1101 1102 Output Parameters: 1103 . vptr - array data structured 1104 1105 Note: The vector values are NOT initialized and may have garbage in them, so you may need 1106 to zero them. 1107 1108 Level: advanced 1109 1110 .seealso: DMDARestoreArray() 1111 1112 @*/ 1113 PetscErrorCode DMDAGetArray(DM da,PetscBool ghosted,void *vptr) 1114 { 1115 PetscErrorCode ierr; 1116 PetscInt j,i,xs,ys,xm,ym,zs,zm; 1117 char *iarray_start; 1118 void **iptr = (void**)vptr; 1119 DM_DA *dd = (DM_DA*)da->data; 1120 1121 PetscFunctionBegin; 1122 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1123 if (ghosted) { 1124 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1125 if (dd->arrayghostedin[i]) { 1126 *iptr = dd->arrayghostedin[i]; 1127 iarray_start = (char*)dd->startghostedin[i]; 1128 dd->arrayghostedin[i] = NULL; 1129 dd->startghostedin[i] = NULL; 1130 1131 goto done; 1132 } 1133 } 1134 xs = dd->Xs; 1135 ys = dd->Ys; 1136 zs = dd->Zs; 1137 xm = dd->Xe-dd->Xs; 1138 ym = dd->Ye-dd->Ys; 1139 zm = dd->Ze-dd->Zs; 1140 } else { 1141 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1142 if (dd->arrayin[i]) { 1143 *iptr = dd->arrayin[i]; 1144 iarray_start = (char*)dd->startin[i]; 1145 dd->arrayin[i] = NULL; 1146 dd->startin[i] = NULL; 1147 1148 goto done; 1149 } 1150 } 1151 xs = dd->xs; 1152 ys = dd->ys; 1153 zs = dd->zs; 1154 xm = dd->xe-dd->xs; 1155 ym = dd->ye-dd->ys; 1156 zm = dd->ze-dd->zs; 1157 } 1158 1159 switch (dd->dim) { 1160 case 1: { 1161 void *ptr; 1162 1163 ierr = PetscMalloc(xm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1164 1165 ptr = (void*)(iarray_start - xs*sizeof(PetscScalar)); 1166 *iptr = (void*)ptr; 1167 break; 1168 } 1169 case 2: { 1170 void **ptr; 1171 1172 ierr = PetscMalloc((ym+1)*sizeof(void*)+xm*ym*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1173 1174 ptr = (void**)(iarray_start + xm*ym*sizeof(PetscScalar) - ys*sizeof(void*)); 1175 for (j=ys; j<ys+ym; j++) ptr[j] = iarray_start + sizeof(PetscScalar)*(xm*(j-ys) - xs); 1176 *iptr = (void*)ptr; 1177 break; 1178 } 1179 case 3: { 1180 void ***ptr,**bptr; 1181 1182 ierr = PetscMalloc((zm+1)*sizeof(void**)+(ym*zm+1)*sizeof(void*)+xm*ym*zm*sizeof(PetscScalar),&iarray_start);CHKERRQ(ierr); 1183 1184 ptr = (void***)(iarray_start + xm*ym*zm*sizeof(PetscScalar) - zs*sizeof(void*)); 1185 bptr = (void**)(iarray_start + xm*ym*zm*sizeof(PetscScalar) + zm*sizeof(void**)); 1186 for (i=zs; i<zs+zm; i++) ptr[i] = bptr + ((i-zs)*ym - ys); 1187 for (i=zs; i<zs+zm; i++) { 1188 for (j=ys; j<ys+ym; j++) { 1189 ptr[i][j] = iarray_start + sizeof(PetscScalar)*(xm*ym*(i-zs) + xm*(j-ys) - xs); 1190 } 1191 } 1192 *iptr = (void*)ptr; 1193 break; 1194 } 1195 default: 1196 SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Dimension %D not supported",dd->dim); 1197 } 1198 1199 done: 1200 /* add arrays to the checked out list */ 1201 if (ghosted) { 1202 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1203 if (!dd->arrayghostedout[i]) { 1204 dd->arrayghostedout[i] = *iptr; 1205 dd->startghostedout[i] = iarray_start; 1206 break; 1207 } 1208 } 1209 } else { 1210 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1211 if (!dd->arrayout[i]) { 1212 dd->arrayout[i] = *iptr; 1213 dd->startout[i] = iarray_start; 1214 break; 1215 } 1216 } 1217 } 1218 PetscFunctionReturn(0); 1219 } 1220 1221 #undef __FUNCT__ 1222 #define __FUNCT__ "DMDARestoreArray" 1223 /*@C 1224 DMDARestoreArray - Restores an array of derivative types for a DMDA 1225 1226 Input Parameter: 1227 + da - information about my local patch 1228 . ghosted - do you want arrays for the ghosted or nonghosted patch 1229 - vptr - array data structured to be passed to ad_FormFunctionLocal() 1230 1231 Level: advanced 1232 1233 .seealso: DMDAGetArray() 1234 1235 @*/ 1236 PetscErrorCode DMDARestoreArray(DM da,PetscBool ghosted,void *vptr) 1237 { 1238 PetscInt i; 1239 void **iptr = (void**)vptr,*iarray_start = 0; 1240 DM_DA *dd = (DM_DA*)da->data; 1241 1242 PetscFunctionBegin; 1243 PetscValidHeaderSpecific(da,DM_CLASSID,1); 1244 if (ghosted) { 1245 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1246 if (dd->arrayghostedout[i] == *iptr) { 1247 iarray_start = dd->startghostedout[i]; 1248 dd->arrayghostedout[i] = NULL; 1249 dd->startghostedout[i] = NULL; 1250 break; 1251 } 1252 } 1253 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1254 if (!dd->arrayghostedin[i]) { 1255 dd->arrayghostedin[i] = *iptr; 1256 dd->startghostedin[i] = iarray_start; 1257 break; 1258 } 1259 } 1260 } else { 1261 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1262 if (dd->arrayout[i] == *iptr) { 1263 iarray_start = dd->startout[i]; 1264 dd->arrayout[i] = NULL; 1265 dd->startout[i] = NULL; 1266 break; 1267 } 1268 } 1269 for (i=0; i<DMDA_MAX_WORK_ARRAYS; i++) { 1270 if (!dd->arrayin[i]) { 1271 dd->arrayin[i] = *iptr; 1272 dd->startin[i] = iarray_start; 1273 break; 1274 } 1275 } 1276 } 1277 PetscFunctionReturn(0); 1278 } 1279 1280