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