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