1 #include <petsc/private/dmfieldimpl.h> /*I "petscdmfield.h" I*/ 2 #include <petsc/private/dmimpl.h> /*I "petscdm.h" I*/ 3 #include <petscdmda.h> 4 5 typedef struct _n_DMField_DA 6 { 7 PetscScalar *cornerVals; 8 PetscScalar *cornerCoeffs; 9 PetscScalar *work; 10 PetscReal coordRange[3][2]; 11 } 12 DMField_DA; 13 14 static PetscErrorCode DMFieldDestroy_DA(DMField field) 15 { 16 DMField_DA *dafield; 17 PetscErrorCode ierr; 18 19 PetscFunctionBegin; 20 dafield = (DMField_DA *) field->data; 21 ierr = PetscFree3(dafield->cornerVals,dafield->cornerCoeffs,dafield->work);CHKERRQ(ierr); 22 ierr = PetscFree(dafield);CHKERRQ(ierr); 23 PetscFunctionReturn(0); 24 } 25 26 static PetscErrorCode DMFieldView_DA(DMField field,PetscViewer viewer) 27 { 28 DMField_DA *dafield = (DMField_DA *) field->data; 29 PetscBool iascii; 30 PetscErrorCode ierr; 31 32 PetscFunctionBegin; 33 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 34 if (iascii) { 35 PetscInt i, c, dim; 36 PetscInt nc; 37 DM dm = field->dm; 38 39 ierr = PetscViewerASCIIPrintf(viewer, "Field corner values:\n");CHKERRQ(ierr); 40 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 41 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 42 nc = field->numComponents; 43 for (i = 0, c = 0; i < (1 << dim); i++) { 44 PetscInt j; 45 46 for (j = 0; j < nc; j++, c++) { 47 PetscScalar val = dafield->cornerVals[nc * i + j]; 48 49 #if !defined(PETSC_USE_COMPLEX) 50 ierr = PetscViewerASCIIPrintf(viewer,"%g ",(double) val);CHKERRQ(ierr); 51 #else 52 ierr = PetscViewerASCIIPrintf(viewer,"%g+i%g ",(double) PetscRealPart(val),(double) PetscImaginaryPart(val));CHKERRQ(ierr); 53 #endif 54 } 55 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 56 } 57 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 58 } 59 PetscFunctionReturn(0); 60 } 61 62 #define MEdot(y,A,x,m,c,cast) \ 63 do { \ 64 PetscInt _k, _l; \ 65 for (_k = 0; _k < (c); _k++) (y)[_k] = 0.; \ 66 for (_l = 0; _l < (m); _l++) { \ 67 for (_k = 0; _k < (c); _k++) { \ 68 (y)[_k] += cast((A)[(c) * _l + _k]) * (x)[_l]; \ 69 } \ 70 } \ 71 } while (0) 72 73 #define MEHess(out,cf,etaB,etaD,dim,nc,cast) \ 74 do { \ 75 PetscInt _m, _j, _k; \ 76 for (_m = 0; _m < (nc) * (dim) * (dim); _m++) (out)[_m] = 0.; \ 77 for (_j = 0; _j < (dim); _j++) { \ 78 for (_k = _j + 1; _k < (dim); _k++) { \ 79 PetscInt _ind = (1 << _j) + (1 << _k); \ 80 for (_m = 0; _m < (nc); _m++) { \ 81 PetscScalar c = (cf)[_m] * (etaB)[_ind] * (etaD)[_ind]; \ 82 (out)[(_m * (dim) + _k) * (dim) + _j] += cast(c); \ 83 (out)[(_m * (dim) + _j) * (dim) + _k] += cast(c); \ 84 } \ 85 } \ 86 } \ 87 } while (0) 88 89 static void MultilinearEvaluate(PetscInt dim, PetscReal (*coordRange)[2], PetscInt nc, PetscScalar *cf, PetscScalar *cfWork, PetscInt nPoints, const PetscScalar *points, PetscDataType datatype, void *B, void *D, void *H) 90 { 91 PetscInt i, j, k, l, m; 92 PetscInt whol = 1 << dim; 93 PetscInt half = whol >> 1; 94 95 PetscFunctionBeginHot; 96 if (!B && !D && !H) PetscFunctionReturnVoid(); 97 for (i = 0; i < nPoints; i++) { 98 const PetscScalar *point = &points[dim * i]; 99 PetscReal deta[3] = {0.}; 100 PetscReal etaB[8] = {1.,1.,1.,1.,1.,1.,1.,1.}; 101 PetscReal etaD[8] = {1.,1.,1.,1.,1.,1.,1.,1.}; 102 PetscReal work[8]; 103 104 for (j = 0; j < dim; j++) { 105 PetscReal e, d; 106 107 e = (PetscRealPart(point[j]) - coordRange[j][0]) / coordRange[j][1]; 108 deta[j] = d = 1. / coordRange[j][1]; 109 for (k = 0; k < whol; k++) {work[k] = etaB[k];} 110 for (k = 0; k < half; k++) { 111 etaB[k] = work[2 * k] * e; 112 etaB[k + half] = work[2 * k + 1]; 113 } 114 if (H) { 115 for (k = 0; k < whol; k++) {work[k] = etaD[k];} 116 for (k = 0; k < half; k++) { 117 etaD[k + half] = work[2 * k]; 118 etaD[k ] = work[2 * k + 1] * d; 119 } 120 } 121 } 122 if (B) { 123 if (datatype == PETSC_SCALAR) { 124 PetscScalar *out = &((PetscScalar *)B)[nc * i]; 125 126 MEdot(out,cf,etaB,(1 << dim),nc,(PetscScalar)); 127 } else { 128 PetscReal *out = &((PetscReal *)B)[nc * i]; 129 130 MEdot(out,cf,etaB,(1 << dim),nc,PetscRealPart); 131 } 132 } 133 if (D) { 134 if (datatype == PETSC_SCALAR) { 135 PetscScalar *out = &((PetscScalar *)D)[nc * dim * i]; 136 137 for (m = 0; m < nc * dim; m++) out[m] = 0.; 138 } else { 139 PetscReal *out = &((PetscReal *)D)[nc * dim * i]; 140 141 for (m = 0; m < nc * dim; m++) out[m] = 0.; 142 } 143 for (j = 0; j < dim; j++) { 144 PetscReal d = deta[j]; 145 146 for (k = 0; k < whol * nc; k++) {cfWork[k] = cf[k];} 147 for (k = 0; k < whol; k++) {work[k] = etaB[k];} 148 for (k = 0; k < half; k++) { 149 PetscReal e; 150 151 etaB[k] = work[2 * k]; 152 etaB[k + half] = e = work[2 * k + 1]; 153 154 for (l = 0; l < nc; l++) { 155 cf[ k * nc + l] = cfWork[ 2 * k * nc + l]; 156 cf[(k + half) * nc + l] = cfWork[(2 * k + 1) * nc + l]; 157 } 158 if (datatype == PETSC_SCALAR) { 159 PetscScalar *out = &((PetscScalar *)D)[nc * dim * i]; 160 161 for (l = 0; l < nc; l++) { 162 out[l * dim + j] += d * e * cf[k * nc + l]; 163 } 164 } else { 165 PetscReal *out = &((PetscReal *)D)[nc * dim * i]; 166 167 for (l = 0; l < nc; l++) { 168 out[l * dim + j] += d * e * PetscRealPart(cf[k * nc + l]); 169 } 170 } 171 } 172 } 173 } 174 if (H) { 175 if (datatype == PETSC_SCALAR) { 176 PetscScalar *out = &((PetscScalar *)H)[nc * dim * dim * i]; 177 178 MEHess(out,cf,etaB,etaD,dim,nc,(PetscScalar)); 179 } else { 180 PetscReal *out = &((PetscReal *)H)[nc * dim * dim * i]; 181 182 MEHess(out,cf,etaB,etaD,dim,nc,PetscRealPart); 183 } 184 } 185 } 186 PetscFunctionReturnVoid(); 187 } 188 189 static PetscErrorCode DMFieldEvaluate_DA(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H) 190 { 191 DM dm; 192 DMField_DA *dafield; 193 PetscInt dim; 194 PetscInt N, n, nc; 195 const PetscScalar *array; 196 PetscReal (*coordRange)[2]; 197 PetscErrorCode ierr; 198 199 PetscFunctionBegin; 200 dm = field->dm; 201 nc = field->numComponents; 202 dafield = (DMField_DA *) field->data; 203 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 204 ierr = VecGetLocalSize(points,&N);CHKERRQ(ierr); 205 if (N % dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Point vector size %D not divisible by coordinate dimension %D\n",N,dim); 206 n = N / dim; 207 coordRange = &(dafield->coordRange[0]); 208 ierr = VecGetArrayRead(points,&array);CHKERRQ(ierr); 209 MultilinearEvaluate(dim,coordRange,nc,dafield->cornerCoeffs,dafield->work,n,array,datatype,B,D,H); 210 ierr = VecRestoreArrayRead(points,&array);CHKERRQ(ierr); 211 PetscFunctionReturn(0); 212 } 213 214 static PetscErrorCode DMFieldEvaluateFE_DA(DMField field, IS cellIS, PetscQuadrature points, PetscDataType datatype, void *B, void *D, void *H) 215 { 216 PetscInt c, i, j, k, dim, cellsPer[3] = {0}, first[3] = {0}, whol, half; 217 PetscReal stepPer[3] = {0.}; 218 PetscReal cellCoordRange[3][2] = {{0.,1.},{0.,1.},{0.,1.}}; 219 PetscScalar *cellCoeffs, *work; 220 DM dm; 221 DMDALocalInfo info; 222 PetscInt cStart, cEnd; 223 PetscInt nq, nc; 224 const PetscReal *q; 225 #if defined(PETSC_USE_COMPLEX) 226 PetscScalar *qs; 227 #else 228 const PetscScalar *qs; 229 #endif 230 DMField_DA *dafield; 231 PetscBool isStride; 232 const PetscInt *cells = NULL; 233 PetscInt sfirst = -1, stride = -1, nCells; 234 PetscErrorCode ierr; 235 236 PetscFunctionBegin; 237 dafield = (DMField_DA *) field->data; 238 dm = field->dm; 239 nc = field->numComponents; 240 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 241 dim = info.dim; 242 work = dafield->work; 243 stepPer[0] = 1./ info.mx; 244 stepPer[1] = 1./ info.my; 245 stepPer[2] = 1./ info.mz; 246 first[0] = info.gxs; 247 first[1] = info.gys; 248 first[2] = info.gzs; 249 cellsPer[0] = info.gxm; 250 cellsPer[1] = info.gym; 251 cellsPer[2] = info.gzm; 252 /* TODO: probably take components into account */ 253 ierr = PetscQuadratureGetData(points, NULL, NULL, &nq, &q, NULL);CHKERRQ(ierr); 254 #if defined(PETSC_USE_COMPLEX) 255 ierr = DMGetWorkArray(dm,nq * dim,MPIU_SCALAR,&qs);CHKERRQ(ierr); 256 for (i = 0; i < nq * dim; i++) qs[i] = q[i]; 257 #else 258 qs = q; 259 #endif 260 ierr = DMDAGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 261 ierr = DMGetWorkArray(dm,(1 << dim) * nc,MPIU_SCALAR,&cellCoeffs);CHKERRQ(ierr); 262 whol = (1 << dim); 263 half = whol >> 1; 264 ierr = ISGetLocalSize(cellIS,&nCells);CHKERRQ(ierr); 265 ierr = PetscObjectTypeCompare((PetscObject)cellIS,ISSTRIDE,&isStride);CHKERRQ(ierr); 266 if (isStride) { 267 ierr = ISStrideGetInfo(cellIS,&sfirst,&stride);CHKERRQ(ierr); 268 } else { 269 ierr = ISGetIndices(cellIS,&cells);CHKERRQ(ierr); 270 } 271 for (c = 0; c < nCells; c++) { 272 PetscInt cell = isStride ? (sfirst + c * stride) : cells[c]; 273 PetscInt rem = cell; 274 PetscInt ijk[3] = {0}; 275 void *cB, *cD, *cH; 276 277 if (datatype == PETSC_SCALAR) { 278 cB = B ? &((PetscScalar *)B)[nc * nq * c] : NULL; 279 cD = D ? &((PetscScalar *)D)[nc * nq * dim * c] : NULL; 280 cH = H ? &((PetscScalar *)H)[nc * nq * dim * dim * c] : NULL; 281 } else { 282 cB = B ? &((PetscReal *)B)[nc * nq * c] : NULL; 283 cD = D ? &((PetscReal *)D)[nc * nq * dim * c] : NULL; 284 cH = H ? &((PetscReal *)H)[nc * nq * dim * dim * c] : NULL; 285 } 286 if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Point %D not a cell [%D,%D), not implemented yet",cell,cStart,cEnd); 287 for (i = 0; i < nc * whol; i++) {work[i] = dafield->cornerCoeffs[i];} 288 for (j = 0; j < dim; j++) { 289 PetscReal e, d; 290 ijk[j] = (rem % cellsPer[j]); 291 rem /= cellsPer[j]; 292 293 e = 2. * (ijk[j] + first[j] + 0.5) * stepPer[j] - 1.; 294 d = stepPer[j]; 295 for (i = 0; i < half; i++) { 296 for (k = 0; k < nc; k++) { 297 cellCoeffs[ i * nc + k] = work[ 2 * i * nc + k] * d; 298 cellCoeffs[(i + half) * nc + k] = work[ 2 * i * nc + k] * e + work[(2 * i + 1) * nc + k]; 299 } 300 } 301 for (i = 0; i < whol * nc; i++) {work[i] = cellCoeffs[i];} 302 } 303 MultilinearEvaluate(dim,cellCoordRange,nc,cellCoeffs,dafield->work,nq,qs,datatype,cB,cD,cH); 304 } 305 if (!isStride) { 306 ierr = ISRestoreIndices(cellIS,&cells);CHKERRQ(ierr); 307 } 308 ierr = DMRestoreWorkArray(dm,(1 << dim) * nc,MPIU_SCALAR,&cellCoeffs);CHKERRQ(ierr); 309 #if defined(PETSC_USE_COMPLEX) 310 ierr = DMRestoreWorkArray(dm,nq * dim,MPIU_SCALAR,&qs);CHKERRQ(ierr); 311 #endif 312 PetscFunctionReturn(0); 313 } 314 315 static PetscErrorCode DMFieldEvaluateFV_DA(DMField field, IS cellIS, PetscDataType datatype, void *B, void *D, void *H) 316 { 317 PetscInt c, i, dim, cellsPer[3] = {0}, first[3] = {0}; 318 PetscReal stepPer[3] = {0.}; 319 DM dm; 320 DMDALocalInfo info; 321 PetscInt cStart, cEnd, numCells; 322 PetscInt nc; 323 PetscScalar *points; 324 DMField_DA *dafield; 325 PetscBool isStride; 326 const PetscInt *cells = NULL; 327 PetscInt sfirst = -1, stride = -1; 328 PetscErrorCode ierr; 329 330 PetscFunctionBegin; 331 dafield = (DMField_DA *) field->data; 332 dm = field->dm; 333 nc = field->numComponents; 334 ierr = DMDAGetLocalInfo(dm,&info);CHKERRQ(ierr); 335 dim = info.dim; 336 stepPer[0] = 1./ info.mx; 337 stepPer[1] = 1./ info.my; 338 stepPer[2] = 1./ info.mz; 339 first[0] = info.gxs; 340 first[1] = info.gys; 341 first[2] = info.gzs; 342 cellsPer[0] = info.gxm; 343 cellsPer[1] = info.gym; 344 cellsPer[2] = info.gzm; 345 ierr = DMDAGetHeightStratum(dm,0,&cStart,&cEnd);CHKERRQ(ierr); 346 ierr = ISGetLocalSize(cellIS,&numCells);CHKERRQ(ierr); 347 ierr = DMGetWorkArray(dm,dim * numCells,MPIU_SCALAR,&points);CHKERRQ(ierr); 348 ierr = PetscObjectTypeCompare((PetscObject)cellIS,ISSTRIDE,&isStride);CHKERRQ(ierr); 349 if (isStride) { 350 ierr = ISStrideGetInfo(cellIS,&sfirst,&stride);CHKERRQ(ierr); 351 } else { 352 ierr = ISGetIndices(cellIS,&cells);CHKERRQ(ierr); 353 } 354 for (c = 0; c < numCells; c++) { 355 PetscInt cell = isStride ? (sfirst + c * stride) : cells[c]; 356 PetscInt rem = cell; 357 PetscInt ijk[3] = {0}; 358 359 if (cell < cStart || cell >= cEnd) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Point %D not a cell [%D,%D), not implemented yet",cell,cStart,cEnd); 360 for (i = 0; i < dim; i++) { 361 ijk[i] = (rem % cellsPer[i]); 362 rem /= cellsPer[i]; 363 points[dim * c + i] = (ijk[i] + first[i] + 0.5) * stepPer[i]; 364 } 365 } 366 if (!isStride) { 367 ierr = ISRestoreIndices(cellIS,&cells);CHKERRQ(ierr); 368 } 369 MultilinearEvaluate(dim,dafield->coordRange,nc,dafield->cornerCoeffs,dafield->work,numCells,points,datatype,B,D,H); 370 ierr = DMRestoreWorkArray(dm,dim * numCells,MPIU_SCALAR,&points);CHKERRQ(ierr); 371 PetscFunctionReturn(0); 372 } 373 374 static PetscErrorCode DMFieldGetDegree_DA(DMField field, IS pointIS, PetscInt *minDegree, PetscInt *maxDegree) 375 { 376 DM dm; 377 PetscInt dim, h, imin; 378 PetscErrorCode ierr; 379 380 PetscFunctionBegin; 381 dm = field->dm; 382 ierr = ISGetMinMax(pointIS,&imin,NULL);CHKERRQ(ierr); 383 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 384 for (h = 0; h <= dim; h++) { 385 PetscInt hEnd; 386 387 ierr = DMDAGetHeightStratum(dm,h,NULL,&hEnd);CHKERRQ(ierr); 388 if (imin < hEnd) break; 389 } 390 dim -= h; 391 if (minDegree) *minDegree = 1; 392 if (maxDegree) *maxDegree = dim; 393 PetscFunctionReturn(0); 394 } 395 396 static PetscErrorCode DMFieldCreateDefaultQuadrature_DA(DMField field, IS cellIS, PetscQuadrature *quad) 397 { 398 PetscInt h, dim, imax, imin; 399 DM dm; 400 PetscErrorCode ierr; 401 402 PetscFunctionBegin; 403 dm = field->dm; 404 ierr = ISGetMinMax(cellIS,&imin,&imax);CHKERRQ(ierr); 405 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 406 *quad = NULL; 407 for (h = 0; h <= dim; h++) { 408 PetscInt hStart, hEnd; 409 410 ierr = DMDAGetHeightStratum(dm,h,&hStart,&hEnd);CHKERRQ(ierr); 411 if (imin >= hStart && imax < hEnd) break; 412 } 413 dim -= h; 414 if (dim > 0) { 415 ierr = PetscDTGaussTensorQuadrature(dim, 1, 1, -1.0, 1.0, quad);CHKERRQ(ierr); 416 } 417 418 PetscFunctionReturn(0); 419 } 420 421 static PetscErrorCode DMFieldInitialize_DA(DMField field) 422 { 423 DM dm; 424 Vec coords = NULL; 425 PetscInt dim, i, j, k; 426 DMField_DA *dafield = (DMField_DA *) field->data; 427 PetscErrorCode ierr; 428 429 PetscFunctionBegin; 430 field->ops->destroy = DMFieldDestroy_DA; 431 field->ops->evaluate = DMFieldEvaluate_DA; 432 field->ops->evaluateFE = DMFieldEvaluateFE_DA; 433 field->ops->evaluateFV = DMFieldEvaluateFV_DA; 434 field->ops->getDegree = DMFieldGetDegree_DA; 435 field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DA; 436 field->ops->view = DMFieldView_DA; 437 dm = field->dm; 438 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 439 if (dm->coordinates) coords = dm->coordinates; 440 else if (dm->coordinatesLocal) coords = dm->coordinatesLocal; 441 if (coords) { 442 PetscInt n; 443 const PetscScalar *array; 444 PetscReal mins[3][2] = {{PETSC_MAX_REAL,PETSC_MAX_REAL},{PETSC_MAX_REAL,PETSC_MAX_REAL},{PETSC_MAX_REAL,PETSC_MAX_REAL}}; 445 446 ierr = VecGetLocalSize(coords,&n);CHKERRQ(ierr); 447 n /= dim; 448 ierr = VecGetArrayRead(coords,&array);CHKERRQ(ierr); 449 for (i = 0, k = 0; i < n; i++) { 450 for (j = 0; j < dim; j++, k++) { 451 PetscReal val = PetscRealPart(array[k]); 452 453 mins[j][0] = PetscMin(mins[j][0],val); 454 mins[j][1] = PetscMin(mins[j][1],-val); 455 } 456 } 457 ierr = VecRestoreArrayRead(coords,&array);CHKERRQ(ierr); 458 ierr = MPIU_Allreduce((PetscReal *) mins,&(dafield->coordRange[0][0]),2*dim,MPIU_REAL,MPI_MIN,PetscObjectComm((PetscObject)dm));CHKERRMPI(ierr); 459 for (j = 0; j < dim; j++) { 460 dafield->coordRange[j][1] = -dafield->coordRange[j][1]; 461 } 462 } else { 463 for (j = 0; j < dim; j++) { 464 dafield->coordRange[j][0] = 0.; 465 dafield->coordRange[j][1] = 1.; 466 } 467 } 468 for (j = 0; j < dim; j++) { 469 PetscReal avg = 0.5 * (dafield->coordRange[j][1] + dafield->coordRange[j][0]); 470 PetscReal dif = 0.5 * (dafield->coordRange[j][1] - dafield->coordRange[j][0]); 471 472 dafield->coordRange[j][0] = avg; 473 dafield->coordRange[j][1] = dif; 474 } 475 PetscFunctionReturn(0); 476 } 477 478 PETSC_INTERN PetscErrorCode DMFieldCreate_DA(DMField field) 479 { 480 DMField_DA *dafield; 481 PetscErrorCode ierr; 482 483 PetscFunctionBegin; 484 ierr = PetscNewLog(field,&dafield);CHKERRQ(ierr); 485 field->data = dafield; 486 ierr = DMFieldInitialize_DA(field);CHKERRQ(ierr); 487 PetscFunctionReturn(0); 488 } 489 490 PetscErrorCode DMFieldCreateDA(DM dm, PetscInt nc, const PetscScalar *cornerValues,DMField *field) 491 { 492 DMField b; 493 DMField_DA *dafield; 494 PetscInt dim, nv, i, j, k; 495 PetscInt half; 496 PetscScalar *cv, *cf, *work; 497 PetscErrorCode ierr; 498 499 PetscFunctionBegin; 500 ierr = DMFieldCreate(dm,nc,DMFIELD_VERTEX,&b);CHKERRQ(ierr); 501 ierr = DMFieldSetType(b,DMFIELDDA);CHKERRQ(ierr); 502 dafield = (DMField_DA *) b->data; 503 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 504 nv = (1 << dim) * nc; 505 ierr = PetscMalloc3(nv,&cv,nv,&cf,nv,&work);CHKERRQ(ierr); 506 for (i = 0; i < nv; i++) cv[i] = cornerValues[i]; 507 for (i = 0; i < nv; i++) cf[i] = cv[i]; 508 dafield->cornerVals = cv; 509 dafield->cornerCoeffs = cf; 510 dafield->work = work; 511 half = (1 << (dim - 1)); 512 for (i = 0; i < dim; i++) { 513 PetscScalar *w; 514 515 w = work; 516 for (j = 0; j < half; j++) { 517 for (k = 0; k < nc; k++) { 518 w[j * nc + k] = 0.5 * (cf[(2 * j + 1) * nc + k] - cf[(2 * j) * nc + k]); 519 } 520 } 521 w = &work[j * nc]; 522 for (j = 0; j < half; j++) { 523 for (k = 0; k < nc; k++) { 524 w[j * nc + k] = 0.5 * (cf[(2 * j) * nc + k] + cf[(2 * j + 1) * nc + k]); 525 } 526 } 527 for (j = 0; j < nv; j++) {cf[j] = work[j];} 528 } 529 *field = b; 530 PetscFunctionReturn(0); 531 } 532