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 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 DMFieldGetFEInvariance_DA(DMField field, IS pointIS, PetscBool *isConstant, PetscBool *isAffine, PetscBool *isQuadratic) 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 (isConstant) *isConstant = (dim < 1) ? PETSC_TRUE : PETSC_FALSE; 392 if (isAffine) *isAffine = (dim < 2) ? PETSC_TRUE : PETSC_FALSE; 393 if (isQuadratic) *isQuadratic = (dim < 3) ? PETSC_TRUE : PETSC_FALSE; 394 PetscFunctionReturn(0); 395 } 396 397 static PetscErrorCode DMFieldCreateDefaultQuadrature_DA(DMField field, IS cellIS, PetscQuadrature *quad) 398 { 399 PetscInt h, dim, imax, imin; 400 DM dm; 401 PetscErrorCode ierr; 402 403 PetscFunctionBegin; 404 dm = field->dm; 405 ierr = ISGetMinMax(cellIS,&imin,&imax);CHKERRQ(ierr); 406 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 407 *quad = NULL; 408 for (h = 0; h <= dim; h++) { 409 PetscInt hStart, hEnd; 410 411 ierr = DMDAGetHeightStratum(dm,h,&hStart,&hEnd);CHKERRQ(ierr); 412 if (imin >= hStart && imax < hEnd) break; 413 } 414 dim -= h; 415 if (dim > 0) { 416 ierr = PetscDTGaussTensorQuadrature(dim, 1, 1, -1.0, 1.0, quad);CHKERRQ(ierr); 417 } 418 419 PetscFunctionReturn(0); 420 } 421 422 static PetscErrorCode DMFieldInitialize_DA(DMField field) 423 { 424 DM dm; 425 Vec coords = NULL; 426 PetscInt dim, i, j, k; 427 DMField_DA *dafield = (DMField_DA *) field->data; 428 PetscErrorCode ierr; 429 430 PetscFunctionBegin; 431 field->ops->destroy = DMFieldDestroy_DA; 432 field->ops->evaluate = DMFieldEvaluate_DA; 433 field->ops->evaluateFE = DMFieldEvaluateFE_DA; 434 field->ops->evaluateFV = DMFieldEvaluateFV_DA; 435 field->ops->getFEInvariance = DMFieldGetFEInvariance_DA; 436 field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DA; 437 field->ops->view = DMFieldView_DA; 438 dm = field->dm; 439 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 440 if (dm->coordinates) coords = dm->coordinates; 441 else if (dm->coordinatesLocal) coords = dm->coordinatesLocal; 442 if (coords) { 443 PetscInt n; 444 const PetscScalar *array; 445 PetscReal mins[3][2] = {{PETSC_MAX_REAL,PETSC_MAX_REAL},{PETSC_MAX_REAL,PETSC_MAX_REAL},{PETSC_MAX_REAL,PETSC_MAX_REAL}}; 446 447 ierr = VecGetLocalSize(coords,&n);CHKERRQ(ierr); 448 n /= dim; 449 ierr = VecGetArrayRead(coords,&array);CHKERRQ(ierr); 450 for (i = 0, k = 0; i < n; i++) { 451 for (j = 0; j < dim; j++, k++) { 452 PetscReal val = PetscRealPart(array[k]); 453 454 mins[j][0] = PetscMin(mins[j][0],val); 455 mins[j][1] = PetscMin(mins[j][1],-val); 456 } 457 } 458 ierr = VecRestoreArrayRead(coords,&array);CHKERRQ(ierr); 459 ierr = MPIU_Allreduce((PetscReal *) mins,&(dafield->coordRange[0][0]),2*dim,MPIU_REAL,MPI_MIN,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 460 for (j = 0; j < dim; j++) { 461 dafield->coordRange[j][1] = -dafield->coordRange[j][1]; 462 } 463 } else { 464 for (j = 0; j < dim; j++) { 465 dafield->coordRange[j][0] = 0.; 466 dafield->coordRange[j][1] = 1.; 467 } 468 } 469 for (j = 0; j < dim; j++) { 470 PetscReal avg = 0.5 * (dafield->coordRange[j][1] + dafield->coordRange[j][0]); 471 PetscReal dif = 0.5 * (dafield->coordRange[j][1] - dafield->coordRange[j][0]); 472 473 dafield->coordRange[j][0] = avg; 474 dafield->coordRange[j][1] = dif; 475 } 476 PetscFunctionReturn(0); 477 } 478 479 PETSC_INTERN PetscErrorCode DMFieldCreate_DA(DMField field) 480 { 481 DMField_DA *dafield; 482 PetscErrorCode ierr; 483 484 PetscFunctionBegin; 485 ierr = PetscNewLog(field,&dafield);CHKERRQ(ierr); 486 field->data = dafield; 487 ierr = DMFieldInitialize_DA(field);CHKERRQ(ierr); 488 PetscFunctionReturn(0); 489 } 490 491 PetscErrorCode DMFieldCreateDA(DM dm, PetscInt nc, const PetscScalar *cornerValues,DMField *field) 492 { 493 DMField b; 494 DMField_DA *dafield; 495 PetscInt dim, nv, i, j, k; 496 PetscInt half; 497 PetscScalar *cv, *cf, *work; 498 PetscErrorCode ierr; 499 500 PetscFunctionBegin; 501 ierr = DMFieldCreate(dm,nc,DMFIELD_VERTEX,&b);CHKERRQ(ierr); 502 ierr = DMFieldSetType(b,DMFIELDDA);CHKERRQ(ierr); 503 dafield = (DMField_DA *) b->data; 504 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 505 nv = (1 << dim) * nc; 506 ierr = PetscMalloc3(nv,&cv,nv,&cf,nv,&work);CHKERRQ(ierr); 507 for (i = 0; i < nv; i++) cv[i] = cornerValues[i]; 508 for (i = 0; i < nv; i++) cf[i] = cv[i]; 509 dafield->cornerVals = cv; 510 dafield->cornerCoeffs = cf; 511 dafield->work = work; 512 half = (1 << (dim - 1)); 513 for (i = 0; i < dim; i++) { 514 PetscScalar *w; 515 516 w = work; 517 for (j = 0; j < half; j++) { 518 for (k = 0; k < nc; k++) { 519 w[j * nc + k] = 0.5 * (cf[(2 * j + 1) * nc + k] - cf[(2 * j) * nc + k]); 520 } 521 } 522 w = &work[j * nc]; 523 for (j = 0; j < half; j++) { 524 for (k = 0; k < nc; k++) { 525 w[j * nc + k] = 0.5 * (cf[(2 * j) * nc + k] + cf[(2 * j + 1) * nc + k]); 526 } 527 } 528 for (j = 0; j < nv; j++) {cf[j] = work[j];} 529 } 530 *field = b; 531 PetscFunctionReturn(0); 532 } 533