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