1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <petscsf.h> 3 4 #include <petsc/private/petscfeimpl.h> 5 #include <petsc/private/petscfvimpl.h> 6 7 /*@ 8 DMPlexGetScale - Get the scale for the specified fundamental unit 9 10 Not collective 11 12 Input Arguments: 13 + dm - the DM 14 - unit - The SI unit 15 16 Output Argument: 17 . scale - The value used to scale all quantities with this unit 18 19 Level: advanced 20 21 .seealso: DMPlexSetScale(), PetscUnit 22 @*/ 23 PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale) 24 { 25 DM_Plex *mesh = (DM_Plex*) dm->data; 26 27 PetscFunctionBegin; 28 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 29 PetscValidPointer(scale, 3); 30 *scale = mesh->scale[unit]; 31 PetscFunctionReturn(0); 32 } 33 34 /*@ 35 DMPlexSetScale - Set the scale for the specified fundamental unit 36 37 Not collective 38 39 Input Arguments: 40 + dm - the DM 41 . unit - The SI unit 42 - scale - The value used to scale all quantities with this unit 43 44 Level: advanced 45 46 .seealso: DMPlexGetScale(), PetscUnit 47 @*/ 48 PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale) 49 { 50 DM_Plex *mesh = (DM_Plex*) dm->data; 51 52 PetscFunctionBegin; 53 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 54 mesh->scale[unit] = scale; 55 PetscFunctionReturn(0); 56 } 57 58 PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k) 59 { 60 switch (i) { 61 case 0: 62 switch (j) { 63 case 0: return 0; 64 case 1: 65 switch (k) { 66 case 0: return 0; 67 case 1: return 0; 68 case 2: return 1; 69 } 70 case 2: 71 switch (k) { 72 case 0: return 0; 73 case 1: return -1; 74 case 2: return 0; 75 } 76 } 77 case 1: 78 switch (j) { 79 case 0: 80 switch (k) { 81 case 0: return 0; 82 case 1: return 0; 83 case 2: return -1; 84 } 85 case 1: return 0; 86 case 2: 87 switch (k) { 88 case 0: return 1; 89 case 1: return 0; 90 case 2: return 0; 91 } 92 } 93 case 2: 94 switch (j) { 95 case 0: 96 switch (k) { 97 case 0: return 0; 98 case 1: return 1; 99 case 2: return 0; 100 } 101 case 1: 102 switch (k) { 103 case 0: return -1; 104 case 1: return 0; 105 case 2: return 0; 106 } 107 case 2: return 0; 108 } 109 } 110 return 0; 111 } 112 113 static PetscErrorCode DMPlexProjectRigidBody_Private(PetscInt dim, PetscReal t, const PetscReal X[], PetscInt Nf, PetscScalar *mode, void *ctx) 114 { 115 PetscInt *ctxInt = (PetscInt *) ctx; 116 PetscInt dim2 = ctxInt[0]; 117 PetscInt d = ctxInt[1]; 118 PetscInt i, j, k = dim > 2 ? d - dim : d; 119 120 PetscFunctionBegin; 121 if (dim != dim2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Input dimension %d does not match context dimension %d", dim, dim2); 122 for (i = 0; i < dim; i++) mode[i] = 0.; 123 if (d < dim) { 124 mode[d] = 1.; 125 } else { 126 for (i = 0; i < dim; i++) { 127 for (j = 0; j < dim; j++) { 128 mode[j] += epsilon(i, j, k)*X[i]; 129 } 130 } 131 } 132 PetscFunctionReturn(0); 133 } 134 135 /*@C 136 DMPlexCreateRigidBody - for the default global section, create rigid body modes from coordinates 137 138 Collective on DM 139 140 Input Arguments: 141 . dm - the DM 142 143 Output Argument: 144 . sp - the null space 145 146 Note: This is necessary to take account of Dirichlet conditions on the displacements 147 148 Level: advanced 149 150 .seealso: MatNullSpaceCreate() 151 @*/ 152 PetscErrorCode DMPlexCreateRigidBody(DM dm, MatNullSpace *sp) 153 { 154 MPI_Comm comm; 155 Vec mode[6]; 156 PetscSection section, globalSection; 157 PetscInt dim, dimEmbed, n, m, d, i, j; 158 PetscErrorCode ierr; 159 160 PetscFunctionBegin; 161 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 162 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 163 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr); 164 if (dim == 1) { 165 ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr); 166 PetscFunctionReturn(0); 167 } 168 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 169 ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr); 170 ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr); 171 m = (dim*(dim+1))/2; 172 ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr); 173 ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr); 174 ierr = VecSetUp(mode[0]);CHKERRQ(ierr); 175 for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);} 176 for (d = 0; d < m; d++) { 177 PetscInt ctx[2]; 178 PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal *, PetscInt, PetscScalar *, void *) = DMPlexProjectRigidBody_Private; 179 void *voidctx = (void *) (&ctx[0]); 180 181 ctx[0] = dimEmbed; 182 ctx[1] = d; 183 ierr = DMProjectFunction(dm, 0.0, &func, &voidctx, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 184 } 185 for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);} 186 /* Orthonormalize system */ 187 for (i = dim; i < m; ++i) { 188 PetscScalar dots[6]; 189 190 ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr); 191 for (j = 0; j < i; ++j) dots[j] *= -1.0; 192 ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr); 193 ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr); 194 } 195 ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr); 196 for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);} 197 PetscFunctionReturn(0); 198 } 199 200 /*@ 201 DMPlexSetMaxProjectionHeight - In DMPlexProjectXXXLocal() functions, the projected values of a basis function's dofs 202 are computed by associating the basis function with one of the mesh points in its transitively-closed support, and 203 evaluating the dual space basis of that point. A basis function is associated with the point in its 204 transitively-closed support whose mesh height is highest (w.r.t. DAG height), but not greater than the maximum 205 projection height, which is set with this function. By default, the maximum projection height is zero, which means 206 that only mesh cells are used to project basis functions. A height of one, for example, evaluates a cell-interior 207 basis functions using its cells dual space basis, but all other basis functions with the dual space basis of a face. 208 209 Input Parameters: 210 + dm - the DMPlex object 211 - height - the maximum projection height >= 0 212 213 Level: advanced 214 215 .seealso: DMPlexGetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 216 @*/ 217 PetscErrorCode DMPlexSetMaxProjectionHeight(DM dm, PetscInt height) 218 { 219 DM_Plex *plex = (DM_Plex *) dm->data; 220 221 PetscFunctionBegin; 222 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 223 plex->maxProjectionHeight = height; 224 PetscFunctionReturn(0); 225 } 226 227 /*@ 228 DMPlexGetMaxProjectionHeight - Get the maximum height (w.r.t. DAG) of mesh points used to evaluate dual bases in 229 DMPlexProjectXXXLocal() functions. 230 231 Input Parameters: 232 . dm - the DMPlex object 233 234 Output Parameters: 235 . height - the maximum projection height 236 237 Level: intermediate 238 239 .seealso: DMPlexSetMaxProjectionHeight(), DMProjectFunctionLocal(), DMProjectFunctionLabelLocal() 240 @*/ 241 PetscErrorCode DMPlexGetMaxProjectionHeight(DM dm, PetscInt *height) 242 { 243 DM_Plex *plex = (DM_Plex *) dm->data; 244 245 PetscFunctionBegin; 246 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 247 *height = plex->maxProjectionHeight; 248 PetscFunctionReturn(0); 249 } 250 251 /*@C 252 DMPlexInsertBoundaryValuesEssential - Insert boundary values into a local vector 253 254 Input Parameters: 255 + dm - The DM, with a PetscDS that matches the problem being constrained 256 . time - The time 257 . field - The field to constrain 258 . label - The DMLabel defining constrained points 259 . numids - The number of DMLabel ids for constrained points 260 . ids - An array of ids for constrained points 261 . func - A pointwise function giving boundary values 262 - ctx - An optional user context for bcFunc 263 264 Output Parameter: 265 . locX - A local vector to receives the boundary values 266 267 Level: developer 268 269 .seealso: DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary() 270 @*/ 271 PetscErrorCode DMPlexInsertBoundaryValuesEssential(DM dm, PetscReal time, PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *), void *ctx, Vec locX) 272 { 273 PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal x[], PetscInt, PetscScalar *u, void *ctx); 274 void **ctxs; 275 PetscInt numFields; 276 PetscErrorCode ierr; 277 278 PetscFunctionBegin; 279 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 280 ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 281 funcs[field] = func; 282 ctxs[field] = ctx; 283 ierr = DMProjectFunctionLabelLocal(dm, time, label, numids, ids, funcs, ctxs, INSERT_BC_VALUES, locX);CHKERRQ(ierr); 284 ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr); 285 PetscFunctionReturn(0); 286 } 287 288 /*@C 289 DMPlexInsertBoundaryValuesEssentialField - Insert boundary values into a local vector 290 291 Input Parameters: 292 + dm - The DM, with a PetscDS that matches the problem being constrained 293 . time - The time 294 . locU - A local vector with the input solution values 295 . field - The field to constrain 296 . label - The DMLabel defining constrained points 297 . numids - The number of DMLabel ids for constrained points 298 . ids - An array of ids for constrained points 299 . func - A pointwise function giving boundary values 300 - ctx - An optional user context for bcFunc 301 302 Output Parameter: 303 . locX - A local vector to receives the boundary values 304 305 Level: developer 306 307 .seealso: DMPlexInsertBoundaryValuesEssential(), DMAddBoundary() 308 @*/ 309 PetscErrorCode DMPlexInsertBoundaryValuesEssentialField(DM dm, PetscReal time, Vec locU, PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], 310 void (*func)(PetscInt, PetscInt, PetscInt, 311 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 312 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 313 PetscReal, const PetscReal[], 314 PetscScalar[]), 315 void *ctx, Vec locX) 316 { 317 void (**funcs)(PetscInt, PetscInt, PetscInt, 318 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 319 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 320 PetscReal, const PetscReal[], PetscScalar[]); 321 void **ctxs; 322 PetscInt numFields; 323 PetscErrorCode ierr; 324 325 PetscFunctionBegin; 326 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 327 ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 328 funcs[field] = func; 329 ctxs[field] = ctx; 330 ierr = DMProjectFieldLabelLocal(dm, time, label, numids, ids, locU, funcs, INSERT_BC_VALUES, locX);CHKERRQ(ierr); 331 ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr); 332 PetscFunctionReturn(0); 333 } 334 335 /*@C 336 DMPlexInsertBoundaryValuesRiemann - Insert boundary values into a local vector 337 338 Input Parameters: 339 + dm - The DM, with a PetscDS that matches the problem being constrained 340 . time - The time 341 . faceGeometry - A vector with the FVM face geometry information 342 . cellGeometry - A vector with the FVM cell geometry information 343 . Grad - A vector with the FVM cell gradient information 344 . field - The field to constrain 345 . label - The DMLabel defining constrained points 346 . numids - The number of DMLabel ids for constrained points 347 . ids - An array of ids for constrained points 348 . func - A pointwise function giving boundary values 349 - ctx - An optional user context for bcFunc 350 351 Output Parameter: 352 . locX - A local vector to receives the boundary values 353 354 Note: This implementation currently ignores the numcomps/comps argument from DMAddBoundary() 355 356 Level: developer 357 358 .seealso: DMPlexInsertBoundaryValuesEssential(), DMPlexInsertBoundaryValuesEssentialField(), DMAddBoundary() 359 @*/ 360 PetscErrorCode DMPlexInsertBoundaryValuesRiemann(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad, PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], 361 PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX) 362 { 363 PetscDS prob; 364 PetscSF sf; 365 DM dmFace, dmCell, dmGrad; 366 const PetscScalar *facegeom, *cellgeom = NULL, *grad; 367 const PetscInt *leaves; 368 PetscScalar *x, *fx; 369 PetscInt dim, nleaves, loc, fStart, fEnd, pdim, i; 370 PetscErrorCode ierr, ierru = 0; 371 372 PetscFunctionBegin; 373 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr); 374 ierr = PetscSFGetGraph(sf, NULL, &nleaves, &leaves, NULL);CHKERRQ(ierr); 375 nleaves = PetscMax(0, nleaves); 376 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 377 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 378 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 379 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 380 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 381 if (cellGeometry) { 382 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 383 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 384 } 385 if (Grad) { 386 PetscFV fv; 387 388 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fv);CHKERRQ(ierr); 389 ierr = VecGetDM(Grad, &dmGrad);CHKERRQ(ierr); 390 ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr); 391 ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr); 392 ierr = DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 393 } 394 ierr = VecGetArray(locX, &x);CHKERRQ(ierr); 395 for (i = 0; i < numids; ++i) { 396 IS faceIS; 397 const PetscInt *faces; 398 PetscInt numFaces, f; 399 400 ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr); 401 if (!faceIS) continue; /* No points with that id on this process */ 402 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 403 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 404 for (f = 0; f < numFaces; ++f) { 405 const PetscInt face = faces[f], *cells; 406 PetscFVFaceGeom *fg; 407 408 if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */ 409 ierr = PetscFindInt(face, nleaves, (PetscInt *) leaves, &loc);CHKERRQ(ierr); 410 if (loc >= 0) continue; 411 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 412 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 413 if (Grad) { 414 PetscFVCellGeom *cg; 415 PetscScalar *cx, *cgrad; 416 PetscScalar *xG; 417 PetscReal dx[3]; 418 PetscInt d; 419 420 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr); 421 ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr); 422 ierr = DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr); 423 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 424 DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx); 425 for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx); 426 ierru = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx); 427 if (ierru) { 428 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 429 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 430 goto cleanup; 431 } 432 } else { 433 PetscScalar *xI; 434 PetscScalar *xG; 435 436 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr); 437 ierr = DMPlexPointLocalFieldRef(dm, cells[1], field, x, &xG);CHKERRQ(ierr); 438 ierru = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx); 439 if (ierru) { 440 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 441 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 442 goto cleanup; 443 } 444 } 445 } 446 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 447 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 448 } 449 cleanup: 450 ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr); 451 if (Grad) { 452 ierr = DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 453 ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr); 454 } 455 if (cellGeometry) {ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr);} 456 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 457 CHKERRQ(ierru); 458 PetscFunctionReturn(0); 459 } 460 461 /*@ 462 DMPlexInsertBoundaryValues - Puts coefficients which represent boundary values into the local solution vector 463 464 Input Parameters: 465 + dm - The DM 466 . insertEssential - Should I insert essential (e.g. Dirichlet) or inessential (e.g. Neumann) boundary conditions 467 . time - The time 468 . faceGeomFVM - Face geometry data for FV discretizations 469 . cellGeomFVM - Cell geometry data for FV discretizations 470 - gradFVM - Gradient reconstruction data for FV discretizations 471 472 Output Parameters: 473 . locX - Solution updated with boundary values 474 475 Level: developer 476 477 .seealso: DMProjectFunctionLabelLocal() 478 @*/ 479 PetscErrorCode DMPlexInsertBoundaryValues(DM dm, PetscBool insertEssential, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM) 480 { 481 PetscInt numBd, b; 482 PetscErrorCode ierr; 483 484 PetscFunctionBegin; 485 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 486 PetscValidHeaderSpecific(locX, VEC_CLASSID, 2); 487 if (faceGeomFVM) {PetscValidHeaderSpecific(faceGeomFVM, VEC_CLASSID, 4);} 488 if (cellGeomFVM) {PetscValidHeaderSpecific(cellGeomFVM, VEC_CLASSID, 5);} 489 if (gradFVM) {PetscValidHeaderSpecific(gradFVM, VEC_CLASSID, 6);} 490 ierr = PetscDSGetNumBoundary(dm->prob, &numBd);CHKERRQ(ierr); 491 for (b = 0; b < numBd; ++b) { 492 DMBoundaryConditionType type; 493 const char *labelname; 494 DMLabel label; 495 PetscInt field; 496 PetscObject obj; 497 PetscClassId id; 498 void (*func)(); 499 PetscInt numids; 500 const PetscInt *ids; 501 void *ctx; 502 503 ierr = DMGetBoundary(dm, b, &type, NULL, &labelname, &field, NULL, NULL, &func, &numids, &ids, &ctx);CHKERRQ(ierr); 504 if (insertEssential != (type & DM_BC_ESSENTIAL)) continue; 505 ierr = DMGetLabel(dm, labelname, &label);CHKERRQ(ierr); 506 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 507 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 508 if (id == PETSCFE_CLASSID) { 509 switch (type) { 510 /* for FEM, there is no insertion to be done for non-essential boundary conditions */ 511 case DM_BC_ESSENTIAL: 512 ierr = DMPlexLabelAddCells(dm,label);CHKERRQ(ierr); 513 ierr = DMPlexInsertBoundaryValuesEssential(dm, time, field, label, numids, ids, (PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar *, void *)) func, ctx, locX);CHKERRQ(ierr); 514 ierr = DMPlexLabelClearCells(dm,label);CHKERRQ(ierr); 515 break; 516 case DM_BC_ESSENTIAL_FIELD: 517 ierr = DMPlexLabelAddCells(dm,label);CHKERRQ(ierr); 518 ierr = DMPlexInsertBoundaryValuesEssentialField(dm, time, locX, field, label, numids, ids, 519 (void (*)(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 520 const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], 521 PetscReal, const PetscReal[], PetscScalar[])) func, ctx, locX);CHKERRQ(ierr); 522 ierr = DMPlexLabelClearCells(dm,label);CHKERRQ(ierr); 523 break; 524 default: break; 525 } 526 } else if (id == PETSCFV_CLASSID) { 527 if (!faceGeomFVM) continue; 528 ierr = DMPlexInsertBoundaryValuesRiemann(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, field, label, numids, ids, 529 (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);CHKERRQ(ierr); 530 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 531 } 532 PetscFunctionReturn(0); 533 } 534 535 PetscErrorCode DMComputeL2Diff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 536 { 537 const PetscInt debug = 0; 538 PetscSection section; 539 PetscQuadrature quad; 540 Vec localX; 541 PetscScalar *funcVal, *interpolant; 542 PetscReal *coords, *detJ, *J; 543 PetscReal localDiff = 0.0; 544 const PetscReal *quadWeights; 545 PetscInt dim, coordDim, numFields, numComponents = 0, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset; 546 PetscErrorCode ierr; 547 548 PetscFunctionBegin; 549 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 550 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 551 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 552 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 553 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 554 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 555 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 556 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 557 for (field = 0; field < numFields; ++field) { 558 PetscObject obj; 559 PetscClassId id; 560 PetscInt Nc; 561 562 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 563 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 564 if (id == PETSCFE_CLASSID) { 565 PetscFE fe = (PetscFE) obj; 566 567 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 568 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 569 } else if (id == PETSCFV_CLASSID) { 570 PetscFV fv = (PetscFV) obj; 571 572 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 573 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 574 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 575 numComponents += Nc; 576 } 577 ierr = PetscQuadratureGetData(quad, NULL, &Nq, NULL, &quadWeights);CHKERRQ(ierr); 578 ierr = PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);CHKERRQ(ierr); 579 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 580 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 581 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 582 for (c = cStart; c < cEnd; ++c) { 583 PetscScalar *x = NULL; 584 PetscReal elemDiff = 0.0; 585 586 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);CHKERRQ(ierr); 587 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 588 589 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 590 PetscObject obj; 591 PetscClassId id; 592 void * const ctx = ctxs ? ctxs[field] : NULL; 593 PetscInt Nb, Nc, q, fc; 594 595 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 596 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 597 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 598 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 599 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 600 if (debug) { 601 char title[1024]; 602 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 603 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 604 } 605 for (q = 0; q < Nq; ++q) { 606 if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, point %D", detJ[q], c, q); 607 ierr = (*funcs[field])(coordDim, time, &coords[coordDim * q], Nc, funcVal, ctx); 608 if (ierr) { 609 PetscErrorCode ierr2; 610 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 611 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 612 ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2); 613 CHKERRQ(ierr); 614 } 615 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 616 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 617 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 618 for (fc = 0; fc < Nc; ++fc) { 619 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ[q]);CHKERRQ(ierr);} 620 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ[q]; 621 } 622 } 623 fieldOffset += Nb*Nc; 624 } 625 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 626 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 627 localDiff += elemDiff; 628 } 629 ierr = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr); 630 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 631 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 632 *diff = PetscSqrtReal(*diff); 633 PetscFunctionReturn(0); 634 } 635 636 PetscErrorCode DMComputeL2GradientDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff) 637 { 638 const PetscInt debug = 0; 639 PetscSection section; 640 PetscQuadrature quad; 641 Vec localX; 642 PetscScalar *funcVal, *interpolantVec; 643 PetscReal *coords, *realSpaceDer, *J, *invJ, *detJ; 644 PetscReal localDiff = 0.0; 645 PetscInt dim, coordDim, Nq = 0, numFields, numComponents = 0, cStart, cEnd, cEndInterior, c, field, fieldOffset, comp; 646 PetscErrorCode ierr; 647 648 PetscFunctionBegin; 649 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 650 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 651 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 652 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 653 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 654 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 655 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 656 for (field = 0; field < numFields; ++field) { 657 PetscFE fe; 658 PetscInt Nc; 659 660 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 661 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 662 ierr = PetscQuadratureGetData(quad, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 663 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 664 numComponents += Nc; 665 } 666 /* ierr = DMProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 667 ierr = PetscMalloc7(numComponents,&funcVal,coordDim*Nq,&coords,coordDim*Nq,&realSpaceDer,coordDim*coordDim*Nq,&J,coordDim*coordDim*Nq,&invJ,coordDim,&interpolantVec,Nq,&detJ);CHKERRQ(ierr); 668 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 669 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 670 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 671 for (c = cStart; c < cEnd; ++c) { 672 PetscScalar *x = NULL; 673 PetscReal elemDiff = 0.0; 674 675 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, invJ, detJ);CHKERRQ(ierr); 676 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 677 678 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 679 PetscFE fe; 680 void * const ctx = ctxs ? ctxs[field] : NULL; 681 const PetscReal *quadPoints, *quadWeights; 682 PetscReal *basisDer; 683 PetscInt numQuadPoints, Nb, Ncomp, q, d, fc, f, g; 684 685 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 686 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 687 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 688 ierr = PetscFEGetNumComponents(fe, &Ncomp);CHKERRQ(ierr); 689 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 690 if (debug) { 691 char title[1024]; 692 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 693 ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr); 694 } 695 for (q = 0; q < numQuadPoints; ++q) { 696 if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", detJ[q], c, q); 697 ierr = (*funcs[field])(coordDim, time, &coords[q*coordDim], n, numFields, funcVal, ctx); 698 if (ierr) { 699 PetscErrorCode ierr2; 700 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 701 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 702 ierr2 = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolantVec,detJ);CHKERRQ(ierr2); 703 CHKERRQ(ierr); 704 } 705 for (fc = 0; fc < Ncomp; ++fc) { 706 PetscScalar interpolant = 0.0; 707 708 for (d = 0; d < coordDim; ++d) interpolantVec[d] = 0.0; 709 for (f = 0; f < Nb; ++f) { 710 const PetscInt fidx = f*Ncomp+fc; 711 712 for (d = 0; d < coordDim; ++d) { 713 realSpaceDer[d] = 0.0; 714 for (g = 0; g < dim; ++g) { 715 realSpaceDer[d] += invJ[coordDim * coordDim * q + g*coordDim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g]; 716 } 717 interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d]; 718 } 719 } 720 for (d = 0; d < coordDim; ++d) interpolant += interpolantVec[d]*n[d]; 721 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ[q]);CHKERRQ(ierr);} 722 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ[q]; 723 } 724 } 725 comp += Ncomp; 726 fieldOffset += Nb*Ncomp; 727 } 728 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 729 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 730 localDiff += elemDiff; 731 } 732 ierr = PetscFree7(funcVal,coords,realSpaceDer,J,invJ,interpolantVec,detJ);CHKERRQ(ierr); 733 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 734 ierr = MPIU_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 735 *diff = PetscSqrtReal(*diff); 736 PetscFunctionReturn(0); 737 } 738 739 PetscErrorCode DMComputeL2FieldDiff_Plex(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 740 { 741 const PetscInt debug = 0; 742 PetscSection section; 743 PetscQuadrature quad; 744 Vec localX; 745 PetscScalar *funcVal, *interpolant; 746 PetscReal *coords, *detJ, *J; 747 PetscReal *localDiff; 748 const PetscReal *quadPoints, *quadWeights; 749 PetscInt dim, coordDim, numFields, numComponents = 0, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset; 750 PetscErrorCode ierr; 751 752 PetscFunctionBegin; 753 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 754 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 755 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 756 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 757 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 758 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 759 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 760 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 761 for (field = 0; field < numFields; ++field) { 762 PetscObject obj; 763 PetscClassId id; 764 PetscInt Nc; 765 766 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 767 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 768 if (id == PETSCFE_CLASSID) { 769 PetscFE fe = (PetscFE) obj; 770 771 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 772 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 773 } else if (id == PETSCFV_CLASSID) { 774 PetscFV fv = (PetscFV) obj; 775 776 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 777 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 778 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 779 numComponents += Nc; 780 } 781 ierr = PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 782 ierr = PetscCalloc6(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);CHKERRQ(ierr); 783 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 784 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 785 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 786 for (c = cStart; c < cEnd; ++c) { 787 PetscScalar *x = NULL; 788 789 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);CHKERRQ(ierr); 790 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 791 792 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 793 PetscObject obj; 794 PetscClassId id; 795 void * const ctx = ctxs ? ctxs[field] : NULL; 796 PetscInt Nb, Nc, q, fc; 797 798 PetscReal elemDiff = 0.0; 799 800 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 801 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 802 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 803 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 804 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 805 if (debug) { 806 char title[1024]; 807 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 808 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 809 } 810 for (q = 0; q < Nq; ++q) { 811 if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature point %D", detJ, c, q); 812 ierr = (*funcs[field])(coordDim, time, &coords[coordDim*q], numFields, funcVal, ctx); 813 if (ierr) { 814 PetscErrorCode ierr2; 815 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 816 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 817 ierr2 = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr2); 818 CHKERRQ(ierr); 819 } 820 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 821 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 822 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 823 for (fc = 0; fc < Nc; ++fc) { 824 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d point %g %g %g diff %g\n", c, field, coordDim > 0 ? coords[0] : 0., coordDim > 1 ? coords[1] : 0., coordDim > 2 ? coords[2] : 0., PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ[q]);CHKERRQ(ierr);} 825 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ[q]; 826 } 827 } 828 fieldOffset += Nb*Nc; 829 localDiff[field] += elemDiff; 830 } 831 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 832 } 833 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 834 ierr = MPIU_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 835 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 836 ierr = PetscFree6(localDiff,funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr); 837 PetscFunctionReturn(0); 838 } 839 840 /*@C 841 DMPlexComputeL2DiffVec - This function computes the cellwise L_2 difference between a function u and an FEM interpolant solution u_h, and stores it in a Vec. 842 843 Input Parameters: 844 + dm - The DM 845 . time - The time 846 . funcs - The functions to evaluate for each field component: NULL means that component does not contribute to error calculation 847 . ctxs - Optional array of contexts to pass to each function, or NULL. 848 - X - The coefficient vector u_h 849 850 Output Parameter: 851 . D - A Vec which holds the difference ||u - u_h||_2 for each cell 852 853 Level: developer 854 855 .seealso: DMProjectFunction(), DMComputeL2Diff(), DMPlexComputeL2FieldDiff(), DMComputeL2GradientDiff() 856 @*/ 857 PetscErrorCode DMPlexComputeL2DiffVec(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, Vec D) 858 { 859 PetscSection section; 860 PetscQuadrature quad; 861 Vec localX; 862 PetscScalar *funcVal, *interpolant; 863 PetscReal *coords, *detJ, *J; 864 const PetscReal *quadPoints, *quadWeights; 865 PetscInt dim, coordDim, numFields, numComponents = 0, Nq, cStart, cEnd, cEndInterior, c, field, fieldOffset; 866 PetscErrorCode ierr; 867 868 PetscFunctionBegin; 869 ierr = VecSet(D, 0.0);CHKERRQ(ierr); 870 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 871 ierr = DMGetCoordinateDim(dm, &coordDim);CHKERRQ(ierr); 872 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 873 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 874 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 875 ierr = DMProjectFunctionLocal(dm, time, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 876 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 877 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 878 for (field = 0; field < numFields; ++field) { 879 PetscObject obj; 880 PetscClassId id; 881 PetscInt Nc; 882 883 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 884 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 885 if (id == PETSCFE_CLASSID) { 886 PetscFE fe = (PetscFE) obj; 887 888 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 889 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 890 } else if (id == PETSCFV_CLASSID) { 891 PetscFV fv = (PetscFV) obj; 892 893 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 894 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 895 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 896 numComponents += Nc; 897 } 898 ierr = PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr); 899 ierr = PetscMalloc5(numComponents,&funcVal,numComponents,&interpolant,coordDim*Nq,&coords,Nq,&detJ,coordDim*coordDim*Nq,&J);CHKERRQ(ierr); 900 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 901 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 902 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 903 for (c = cStart; c < cEnd; ++c) { 904 PetscScalar *x = NULL; 905 PetscScalar elemDiff = 0.0; 906 907 ierr = DMPlexComputeCellGeometryFEM(dm, c, quad, coords, J, NULL, detJ);CHKERRQ(ierr); 908 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 909 910 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 911 PetscObject obj; 912 PetscClassId id; 913 void * const ctx = ctxs ? ctxs[field] : NULL; 914 PetscInt Nb, Nc, q, fc; 915 916 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 917 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 918 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 919 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 920 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 921 if (funcs[field]) { 922 for (q = 0; q < Nq; ++q) { 923 if (detJ[q] <= 0.0) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %D, quadrature points %D", detJ[q], c, q); 924 ierr = (*funcs[field])(dim, time, &coords[q*coordDim], Nc, funcVal, ctx); 925 if (ierr) { 926 PetscErrorCode ierr2; 927 ierr2 = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr2); 928 ierr2 = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr); 929 ierr2 = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr2); 930 CHKERRQ(ierr); 931 } 932 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 933 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 934 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 935 for (fc = 0; fc < Nc; ++fc) { 936 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ[q]; 937 } 938 } 939 } 940 fieldOffset += Nb*Nc; 941 } 942 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 943 ierr = VecSetValue(D, c - cStart, elemDiff, INSERT_VALUES);CHKERRQ(ierr); 944 } 945 ierr = PetscFree5(funcVal,interpolant,coords,detJ,J);CHKERRQ(ierr); 946 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 947 ierr = VecSqrtAbs(D);CHKERRQ(ierr); 948 PetscFunctionReturn(0); 949 } 950 951 /*@ 952 DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user 953 954 Input Parameters: 955 + dm - The mesh 956 . X - Local input vector 957 - user - The user context 958 959 Output Parameter: 960 . integral - Local integral for each field 961 962 Level: developer 963 964 .seealso: DMPlexComputeResidualFEM() 965 @*/ 966 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user) 967 { 968 DM_Plex *mesh = (DM_Plex *) dm->data; 969 DM dmAux, dmGrad; 970 Vec localX, A, cellGeometryFVM = NULL, faceGeometryFVM = NULL, locGrad = NULL; 971 PetscDS prob, probAux = NULL; 972 PetscSection section, sectionAux; 973 PetscFV fvm = NULL; 974 PetscFECellGeom *cgeomFEM; 975 PetscFVFaceGeom *fgeomFVM; 976 PetscFVCellGeom *cgeomFVM; 977 PetscScalar *u, *a = NULL; 978 const PetscScalar *lgrad; 979 PetscReal *lintegral; 980 PetscInt *uOff, *uOff_x, *aOff = NULL; 981 PetscInt dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, cEndInterior, c; 982 PetscInt totDim, totDimAux; 983 PetscBool useFVM = PETSC_FALSE; 984 PetscErrorCode ierr; 985 986 PetscFunctionBegin; 987 ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 988 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 989 ierr = DMPlexInsertBoundaryValues(dm, PETSC_TRUE, localX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 990 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 991 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 992 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 993 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 994 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 995 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 996 ierr = PetscDSGetComponentOffsets(prob, &uOff);CHKERRQ(ierr); 997 ierr = PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);CHKERRQ(ierr); 998 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 999 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1000 ierr = DMPlexGetHybridBounds(dm, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1001 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1002 numCells = cEnd - cStart; 1003 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1004 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1005 if (dmAux) { 1006 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1007 ierr = PetscDSGetNumFields(probAux, &NfAux);CHKERRQ(ierr); 1008 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1009 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1010 ierr = PetscDSGetComponentOffsets(probAux, &aOff);CHKERRQ(ierr); 1011 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, NULL);CHKERRQ(ierr); 1012 } 1013 for (f = 0; f < Nf; ++f) { 1014 PetscObject obj; 1015 PetscClassId id; 1016 1017 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1018 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1019 if (id == PETSCFV_CLASSID) {useFVM = PETSC_TRUE; fvm = (PetscFV) obj;} 1020 } 1021 if (useFVM) { 1022 Vec grad; 1023 PetscInt fStart, fEnd; 1024 PetscBool compGrad; 1025 1026 ierr = PetscFVGetComputeGradients(fvm, &compGrad);CHKERRQ(ierr); 1027 ierr = PetscFVSetComputeGradients(fvm, PETSC_TRUE);CHKERRQ(ierr); 1028 ierr = DMPlexComputeGeometryFVM(dm, &cellGeometryFVM, &faceGeometryFVM);CHKERRQ(ierr); 1029 ierr = DMPlexComputeGradientFVM(dm, fvm, faceGeometryFVM, cellGeometryFVM, &dmGrad);CHKERRQ(ierr); 1030 ierr = PetscFVSetComputeGradients(fvm, compGrad);CHKERRQ(ierr); 1031 ierr = VecGetArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);CHKERRQ(ierr); 1032 ierr = VecGetArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1033 /* Reconstruct and limit cell gradients */ 1034 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 1035 ierr = DMGetGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1036 ierr = DMPlexReconstructGradients_Internal(dm, fvm, fStart, fEnd, faceGeometryFVM, cellGeometryFVM, localX, grad);CHKERRQ(ierr); 1037 /* Communicate gradient values */ 1038 ierr = DMGetLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1039 ierr = DMGlobalToLocalBegin(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1040 ierr = DMGlobalToLocalEnd(dmGrad, grad, INSERT_VALUES, locGrad);CHKERRQ(ierr); 1041 ierr = DMRestoreGlobalVector(dmGrad, &grad);CHKERRQ(ierr); 1042 /* Handle non-essential (e.g. outflow) boundary values */ 1043 ierr = DMPlexInsertBoundaryValues(dm, PETSC_FALSE, localX, 0.0, faceGeometryFVM, cellGeometryFVM, locGrad);CHKERRQ(ierr); 1044 ierr = VecGetArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1045 } 1046 ierr = PetscMalloc3(Nf,&lintegral,numCells*totDim,&u,numCells,&cgeomFEM);CHKERRQ(ierr); 1047 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1048 for (f = 0; f < Nf; ++f) {lintegral[f] = 0.0;} 1049 for (c = cStart; c < cEnd; ++c) { 1050 PetscScalar *x = NULL; 1051 PetscInt i; 1052 1053 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeomFEM[c].v0, cgeomFEM[c].J, cgeomFEM[c].invJ, &cgeomFEM[c].detJ);CHKERRQ(ierr); 1054 if (cgeomFEM[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeomFEM[c].detJ, c); 1055 ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1056 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1057 ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1058 if (dmAux) { 1059 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1060 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1061 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1062 } 1063 } 1064 for (f = 0; f < Nf; ++f) { 1065 PetscObject obj; 1066 PetscClassId id; 1067 PetscInt numChunks, numBatches, batchSize, numBlocks, blockSize, Ne, Nr, offset; 1068 1069 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1070 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1071 if (id == PETSCFE_CLASSID) { 1072 PetscFE fe = (PetscFE) obj; 1073 PetscQuadrature q; 1074 PetscInt Nq, Nb; 1075 1076 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1077 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1078 ierr = PetscQuadratureGetData(q, NULL, &Nq, NULL, NULL);CHKERRQ(ierr); 1079 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1080 blockSize = Nb*Nq; 1081 batchSize = numBlocks * blockSize; 1082 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1083 numChunks = numCells / (numBatches*batchSize); 1084 Ne = numChunks*numBatches*batchSize; 1085 Nr = numCells % (numBatches*batchSize); 1086 offset = numCells - Nr; 1087 ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeomFEM, u, probAux, a, lintegral);CHKERRQ(ierr); 1088 ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeomFEM[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], lintegral);CHKERRQ(ierr); 1089 } else if (id == PETSCFV_CLASSID) { 1090 /* PetscFV fv = (PetscFV) obj; */ 1091 PetscInt foff; 1092 PetscPointFunc obj_func; 1093 PetscScalar lint; 1094 1095 ierr = PetscDSGetObjective(prob, f, &obj_func);CHKERRQ(ierr); 1096 ierr = PetscDSGetFieldOffset(prob, f, &foff);CHKERRQ(ierr); 1097 if (obj_func) { 1098 for (c = 0; c < numCells; ++c) { 1099 PetscScalar *u_x; 1100 1101 ierr = DMPlexPointLocalRead(dmGrad, c, lgrad, &u_x);CHKERRQ(ierr); 1102 obj_func(dim, Nf, NfAux, uOff, uOff_x, &u[totDim*c+foff], NULL, u_x, aOff, NULL, a, NULL, NULL, 0.0, cgeomFVM[c].centroid, &lint); 1103 lintegral[f] += PetscRealPart(lint)*cgeomFVM[c].volume; 1104 } 1105 } 1106 } else SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 1107 } 1108 if (useFVM) { 1109 ierr = VecRestoreArrayRead(locGrad, &lgrad);CHKERRQ(ierr); 1110 ierr = VecRestoreArrayRead(faceGeometryFVM, (const PetscScalar **) &fgeomFVM);CHKERRQ(ierr); 1111 ierr = VecRestoreArrayRead(cellGeometryFVM, (const PetscScalar **) &cgeomFVM);CHKERRQ(ierr); 1112 ierr = DMRestoreLocalVector(dmGrad, &locGrad);CHKERRQ(ierr); 1113 ierr = VecDestroy(&faceGeometryFVM);CHKERRQ(ierr); 1114 ierr = VecDestroy(&cellGeometryFVM);CHKERRQ(ierr); 1115 ierr = DMDestroy(&dmGrad);CHKERRQ(ierr); 1116 } 1117 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1118 if (mesh->printFEM) { 1119 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr); 1120 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", lintegral[f]);CHKERRQ(ierr);} 1121 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1122 } 1123 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1124 ierr = MPIU_Allreduce(lintegral, integral, Nf, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr); 1125 ierr = PetscFree3(lintegral,u,cgeomFEM);CHKERRQ(ierr); 1126 ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr); 1127 PetscFunctionReturn(0); 1128 } 1129 1130 /*@ 1131 DMPlexComputeInterpolatorNested - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1132 1133 Input Parameters: 1134 + dmf - The fine mesh 1135 . dmc - The coarse mesh 1136 - user - The user context 1137 1138 Output Parameter: 1139 . In - The interpolation matrix 1140 1141 Level: developer 1142 1143 .seealso: DMPlexComputeInterpolatorGeneral(), DMPlexComputeJacobianFEM() 1144 @*/ 1145 PetscErrorCode DMPlexComputeInterpolatorNested(DM dmc, DM dmf, Mat In, void *user) 1146 { 1147 DM_Plex *mesh = (DM_Plex *) dmc->data; 1148 const char *name = "Interpolator"; 1149 PetscDS prob; 1150 PetscFE *feRef; 1151 PetscFV *fvRef; 1152 PetscSection fsection, fglobalSection; 1153 PetscSection csection, cglobalSection; 1154 PetscScalar *elemMat; 1155 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, cEndInterior, c; 1156 PetscInt cTotDim, rTotDim = 0; 1157 PetscErrorCode ierr; 1158 1159 PetscFunctionBegin; 1160 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1161 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1162 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1163 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1164 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1165 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1166 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1167 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1168 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1169 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1170 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1171 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1172 for (f = 0; f < Nf; ++f) { 1173 PetscObject obj; 1174 PetscClassId id; 1175 PetscInt rNb = 0, Nc = 0; 1176 1177 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1178 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1179 if (id == PETSCFE_CLASSID) { 1180 PetscFE fe = (PetscFE) obj; 1181 1182 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1183 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1184 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1185 } else if (id == PETSCFV_CLASSID) { 1186 PetscFV fv = (PetscFV) obj; 1187 PetscDualSpace Q; 1188 1189 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1190 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1191 ierr = PetscDualSpaceGetDimension(Q, &rNb);CHKERRQ(ierr); 1192 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1193 } 1194 rTotDim += rNb*Nc; 1195 } 1196 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1197 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1198 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1199 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1200 PetscDualSpace Qref; 1201 PetscQuadrature f; 1202 const PetscReal *qpoints, *qweights; 1203 PetscReal *points; 1204 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1205 1206 /* Compose points from all dual basis functionals */ 1207 if (feRef[fieldI]) { 1208 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1209 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1210 } else { 1211 ierr = PetscFVGetDualSpace(fvRef[fieldI], &Qref);CHKERRQ(ierr); 1212 ierr = PetscFVGetNumComponents(fvRef[fieldI], &Nc);CHKERRQ(ierr); 1213 } 1214 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1215 for (i = 0; i < fpdim; ++i) { 1216 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1217 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1218 npoints += Np; 1219 } 1220 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1221 for (i = 0, k = 0; i < fpdim; ++i) { 1222 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1223 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1224 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1225 } 1226 1227 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1228 PetscObject obj; 1229 PetscClassId id; 1230 PetscReal *B; 1231 PetscInt NcJ = 0, cpdim = 0, j; 1232 1233 ierr = PetscDSGetDiscretization(prob, fieldJ, &obj);CHKERRQ(ierr); 1234 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1235 if (id == PETSCFE_CLASSID) { 1236 PetscFE fe = (PetscFE) obj; 1237 1238 /* Evaluate basis at points */ 1239 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1240 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1241 /* For now, fields only interpolate themselves */ 1242 if (fieldI == fieldJ) { 1243 if (Nc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %d does not match coarse field %d", Nc, NcJ); 1244 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1245 for (i = 0, k = 0; i < fpdim; ++i) { 1246 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1247 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1248 for (p = 0; p < Np; ++p, ++k) { 1249 for (j = 0; j < cpdim; ++j) { 1250 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p]; 1251 } 1252 } 1253 } 1254 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1255 } 1256 } else if (id == PETSCFV_CLASSID) { 1257 PetscFV fv = (PetscFV) obj; 1258 1259 /* Evaluate constant function at points */ 1260 ierr = PetscFVGetNumComponents(fv, &NcJ);CHKERRQ(ierr); 1261 cpdim = 1; 1262 /* For now, fields only interpolate themselves */ 1263 if (fieldI == fieldJ) { 1264 if (Nc != NcJ) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %d does not match coarse field %d", Nc, NcJ); 1265 for (i = 0, k = 0; i < fpdim; ++i) { 1266 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1267 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1268 for (p = 0; p < Np; ++p, ++k) { 1269 for (j = 0; j < cpdim; ++j) { 1270 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cTotDim + offsetJ + j*NcJ + c] += 1.0*qweights[p]; 1271 } 1272 } 1273 } 1274 } 1275 } 1276 offsetJ += cpdim*NcJ; 1277 } 1278 offsetI += fpdim*Nc; 1279 ierr = PetscFree(points);CHKERRQ(ierr); 1280 } 1281 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1282 /* Preallocate matrix */ 1283 { 1284 Mat preallocator; 1285 PetscScalar *vals; 1286 PetscInt *cellCIndices, *cellFIndices; 1287 PetscInt locRows, locCols, cell; 1288 1289 ierr = MatGetLocalSize(In, &locRows, &locCols);CHKERRQ(ierr); 1290 ierr = MatCreate(PetscObjectComm((PetscObject) In), &preallocator);CHKERRQ(ierr); 1291 ierr = MatSetType(preallocator, MATPREALLOCATOR);CHKERRQ(ierr); 1292 ierr = MatSetSizes(preallocator, locRows, locCols, PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr); 1293 ierr = MatSetUp(preallocator);CHKERRQ(ierr); 1294 ierr = PetscCalloc3(rTotDim*cTotDim, &vals,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1295 for (cell = cStart; cell < cEnd; ++cell) { 1296 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1297 ierr = MatSetValues(preallocator, rTotDim, cellFIndices, cTotDim, cellCIndices, vals, INSERT_VALUES);CHKERRQ(ierr); 1298 } 1299 ierr = PetscFree3(vals,cellCIndices,cellFIndices);CHKERRQ(ierr); 1300 ierr = MatAssemblyBegin(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1301 ierr = MatAssemblyEnd(preallocator, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1302 ierr = MatPreallocatorPreallocate(preallocator, PETSC_TRUE, In);CHKERRQ(ierr); 1303 ierr = MatDestroy(&preallocator);CHKERRQ(ierr); 1304 } 1305 /* Fill matrix */ 1306 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1307 for (c = cStart; c < cEnd; ++c) { 1308 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1309 } 1310 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1311 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1312 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1313 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1314 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1315 if (mesh->printFEM) { 1316 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1317 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1318 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1319 } 1320 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1321 PetscFunctionReturn(0); 1322 } 1323 1324 /*@ 1325 DMPlexComputeInterpolatorGeneral - Form the local portion of the interpolation matrix I from the coarse DM to a non-nested fine DM. 1326 1327 Input Parameters: 1328 + dmf - The fine mesh 1329 . dmc - The coarse mesh 1330 - user - The user context 1331 1332 Output Parameter: 1333 . In - The interpolation matrix 1334 1335 Level: developer 1336 1337 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1338 @*/ 1339 PetscErrorCode DMPlexComputeInterpolatorGeneral(DM dmc, DM dmf, Mat In, void *user) 1340 { 1341 DM_Plex *mesh = (DM_Plex *) dmf->data; 1342 const char *name = "Interpolator"; 1343 PetscDS prob; 1344 PetscSection fsection, csection, globalFSection, globalCSection; 1345 PetscHashJK ht; 1346 PetscLayout rLayout; 1347 PetscInt *dnz, *onz; 1348 PetscInt locRows, rStart, rEnd; 1349 PetscReal *x, *v0, *J, *invJ, detJ; 1350 PetscReal *v0c, *Jc, *invJc, detJc; 1351 PetscScalar *elemMat; 1352 PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, ccell; 1353 PetscErrorCode ierr; 1354 1355 PetscFunctionBegin; 1356 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1357 ierr = DMGetCoordinateDim(dmc, &dim);CHKERRQ(ierr); 1358 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1359 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 1360 ierr = PetscDSGetNumFields(prob, &Nf);CHKERRQ(ierr); 1361 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 1362 ierr = PetscMalloc3(dim,&v0c,dim*dim,&Jc,dim*dim,&invJc);CHKERRQ(ierr); 1363 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1364 ierr = DMGetDefaultGlobalSection(dmf, &globalFSection);CHKERRQ(ierr); 1365 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1366 ierr = DMGetDefaultGlobalSection(dmc, &globalCSection);CHKERRQ(ierr); 1367 ierr = DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);CHKERRQ(ierr); 1368 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 1369 ierr = PetscMalloc1(totDim*totDim, &elemMat);CHKERRQ(ierr); 1370 1371 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1372 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1373 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1374 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1375 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1376 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1377 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1378 ierr = PetscCalloc2(locRows,&dnz,locRows,&onz);CHKERRQ(ierr); 1379 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1380 for (field = 0; field < Nf; ++field) { 1381 PetscObject obj; 1382 PetscClassId id; 1383 PetscDualSpace Q = NULL; 1384 PetscQuadrature f; 1385 const PetscReal *qpoints; 1386 PetscInt Nc, Np, fpdim, i, d; 1387 1388 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1389 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1390 if (id == PETSCFE_CLASSID) { 1391 PetscFE fe = (PetscFE) obj; 1392 1393 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1394 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1395 } else if (id == PETSCFV_CLASSID) { 1396 PetscFV fv = (PetscFV) obj; 1397 1398 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1399 Nc = 1; 1400 } 1401 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1402 /* For each fine grid cell */ 1403 for (cell = cStart; cell < cEnd; ++cell) { 1404 PetscInt *findices, *cindices; 1405 PetscInt numFIndices, numCIndices; 1406 1407 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1408 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1409 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1410 for (i = 0; i < fpdim; ++i) { 1411 Vec pointVec; 1412 PetscScalar *pV; 1413 PetscSF coarseCellSF = NULL; 1414 const PetscSFNode *coarseCells; 1415 PetscInt numCoarseCells, q, r, c; 1416 1417 /* Get points from the dual basis functional quadrature */ 1418 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1419 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1420 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1421 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1422 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1423 for (q = 0; q < Np; ++q) { 1424 /* Transform point to real space */ 1425 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1426 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1427 } 1428 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1429 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1430 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1431 ierr = PetscSFViewFromOptions(coarseCellSF, NULL, "-interp_sf_view");CHKERRQ(ierr); 1432 /* Update preallocation info */ 1433 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1434 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1435 for (r = 0; r < Nc; ++r) { 1436 PetscHashJKKey key; 1437 PetscHashJKIter missing, iter; 1438 1439 key.j = findices[i*Nc+r]; 1440 if (key.j < 0) continue; 1441 /* Get indices for coarse elements */ 1442 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1443 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1444 for (c = 0; c < numCIndices; ++c) { 1445 key.k = cindices[c]; 1446 if (key.k < 0) continue; 1447 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1448 if (missing) { 1449 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1450 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1451 else ++onz[key.j-rStart]; 1452 } 1453 } 1454 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1455 } 1456 } 1457 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1458 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1459 } 1460 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1461 } 1462 } 1463 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1464 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1465 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1466 ierr = PetscFree2(dnz,onz);CHKERRQ(ierr); 1467 for (field = 0; field < Nf; ++field) { 1468 PetscObject obj; 1469 PetscClassId id; 1470 PetscDualSpace Q = NULL; 1471 PetscQuadrature f; 1472 const PetscReal *qpoints, *qweights; 1473 PetscInt Nc, Np, fpdim, i, d; 1474 1475 ierr = PetscDSGetDiscretization(prob, field, &obj);CHKERRQ(ierr); 1476 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1477 if (id == PETSCFE_CLASSID) { 1478 PetscFE fe = (PetscFE) obj; 1479 1480 ierr = PetscFEGetDualSpace(fe, &Q);CHKERRQ(ierr); 1481 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1482 } else if (id == PETSCFV_CLASSID) { 1483 PetscFV fv = (PetscFV) obj; 1484 1485 ierr = PetscFVGetDualSpace(fv, &Q);CHKERRQ(ierr); 1486 Nc = 1; 1487 } 1488 ierr = PetscDualSpaceGetDimension(Q, &fpdim);CHKERRQ(ierr); 1489 /* For each fine grid cell */ 1490 for (cell = cStart; cell < cEnd; ++cell) { 1491 PetscInt *findices, *cindices; 1492 PetscInt numFIndices, numCIndices; 1493 1494 ierr = DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1495 ierr = DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 1496 if (numFIndices != fpdim*Nc) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of fine indices %d != %d dual basis vecs", numFIndices, fpdim*Nc); 1497 for (i = 0; i < fpdim; ++i) { 1498 Vec pointVec; 1499 PetscScalar *pV; 1500 PetscSF coarseCellSF = NULL; 1501 const PetscSFNode *coarseCells; 1502 PetscInt numCoarseCells, cpdim, q, c, j; 1503 1504 /* Get points from the dual basis functional quadrature */ 1505 ierr = PetscDualSpaceGetFunctional(Q, i, &f);CHKERRQ(ierr); 1506 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, &qweights);CHKERRQ(ierr); 1507 ierr = VecCreateSeq(PETSC_COMM_SELF, Np*dim, &pointVec);CHKERRQ(ierr); 1508 ierr = VecSetBlockSize(pointVec, dim);CHKERRQ(ierr); 1509 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1510 for (q = 0; q < Np; ++q) { 1511 /* Transform point to real space */ 1512 CoordinatesRefToReal(dim, dim, v0, J, &qpoints[q*dim], x); 1513 for (d = 0; d < dim; ++d) pV[q*dim+d] = x[d]; 1514 } 1515 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1516 /* Get set of coarse cells that overlap points (would like to group points by coarse cell) */ 1517 ierr = DMLocatePoints(dmc, pointVec, DM_POINTLOCATION_NEAREST, &coarseCellSF);CHKERRQ(ierr); 1518 /* Update preallocation info */ 1519 ierr = PetscSFGetGraph(coarseCellSF, NULL, &numCoarseCells, NULL, &coarseCells);CHKERRQ(ierr); 1520 if (numCoarseCells != Np) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Not all closure points located"); 1521 ierr = VecGetArray(pointVec, &pV);CHKERRQ(ierr); 1522 for (ccell = 0; ccell < numCoarseCells; ++ccell) { 1523 PetscReal pVReal[3]; 1524 1525 ierr = DMPlexGetClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1526 /* Transform points from real space to coarse reference space */ 1527 ierr = DMPlexComputeCellGeometryFEM(dmc, coarseCells[ccell].index, NULL, v0c, Jc, invJc, &detJc);CHKERRQ(ierr); 1528 for (d = 0; d < dim; ++d) pVReal[d] = PetscRealPart(pV[ccell*dim+d]); 1529 CoordinatesRealToRef(dim, dim, v0c, invJc, pVReal, x); 1530 1531 if (id == PETSCFE_CLASSID) { 1532 PetscFE fe = (PetscFE) obj; 1533 PetscReal *B; 1534 1535 /* Evaluate coarse basis on contained point */ 1536 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1537 ierr = PetscFEGetTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr); 1538 ierr = PetscMemzero(elemMat, cpdim*Nc*Nc * sizeof(PetscScalar));CHKERRQ(ierr); 1539 /* Get elemMat entries by multiplying by weight */ 1540 for (j = 0; j < cpdim; ++j) { 1541 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = B[j*Nc + c]*qweights[ccell]; 1542 } 1543 ierr = PetscFERestoreTabulation(fe, 1, x, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1544 } else { 1545 cpdim = 1; 1546 for (j = 0; j < cpdim; ++j) { 1547 for (c = 0; c < Nc; ++c) elemMat[(c*cpdim + j)*Nc + c] = 1.0*qweights[ccell]; 1548 } 1549 } 1550 /* Update interpolator */ 1551 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(cell, name, Nc, numCIndices, elemMat);CHKERRQ(ierr);} 1552 if (numCIndices != cpdim*Nc) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of element matrix columns %D != %D*%D", numCIndices, cpdim, Nc); 1553 ierr = MatSetValues(In, Nc, &findices[i*Nc], numCIndices, cindices, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1554 ierr = DMPlexRestoreClosureIndices(dmc, csection, globalCSection, coarseCells[ccell].index, &numCIndices, &cindices, NULL);CHKERRQ(ierr); 1555 } 1556 ierr = VecRestoreArray(pointVec, &pV);CHKERRQ(ierr); 1557 ierr = PetscSFDestroy(&coarseCellSF);CHKERRQ(ierr); 1558 ierr = VecDestroy(&pointVec);CHKERRQ(ierr); 1559 } 1560 ierr = DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);CHKERRQ(ierr); 1561 } 1562 } 1563 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 1564 ierr = PetscFree3(v0c,Jc,invJc);CHKERRQ(ierr); 1565 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1566 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1567 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1568 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1569 PetscFunctionReturn(0); 1570 } 1571 1572 /*@ 1573 DMPlexComputeInjectorFEM - Compute a mapping from coarse unknowns to fine unknowns 1574 1575 Input Parameters: 1576 + dmc - The coarse mesh 1577 - dmf - The fine mesh 1578 - user - The user context 1579 1580 Output Parameter: 1581 . sc - The mapping 1582 1583 Level: developer 1584 1585 .seealso: DMPlexComputeInterpolatorNested(), DMPlexComputeJacobianFEM() 1586 @*/ 1587 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 1588 { 1589 PetscDS prob; 1590 PetscFE *feRef; 1591 PetscFV *fvRef; 1592 Vec fv, cv; 1593 IS fis, cis; 1594 PetscSection fsection, fglobalSection, csection, cglobalSection; 1595 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 1596 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, cEndInterior, c, dim, d, startC, endC, offsetC, offsetF, m; 1597 PetscErrorCode ierr; 1598 1599 PetscFunctionBegin; 1600 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1601 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1602 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1603 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1604 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1605 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1606 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1607 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1608 ierr = DMPlexGetHybridBounds(dmc, &cEndInterior, NULL, NULL, NULL);CHKERRQ(ierr); 1609 cEnd = cEndInterior < 0 ? cEnd : cEndInterior; 1610 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1611 ierr = PetscCalloc2(Nf,&feRef,Nf,&fvRef);CHKERRQ(ierr); 1612 for (f = 0; f < Nf; ++f) { 1613 PetscObject obj; 1614 PetscClassId id; 1615 PetscInt fNb = 0, Nc = 0; 1616 1617 ierr = PetscDSGetDiscretization(prob, f, &obj);CHKERRQ(ierr); 1618 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 1619 if (id == PETSCFE_CLASSID) { 1620 PetscFE fe = (PetscFE) obj; 1621 1622 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1623 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 1624 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1625 } else if (id == PETSCFV_CLASSID) { 1626 PetscFV fv = (PetscFV) obj; 1627 PetscDualSpace Q; 1628 1629 ierr = PetscFVRefine(fv, &fvRef[f]);CHKERRQ(ierr); 1630 ierr = PetscFVGetDualSpace(fvRef[f], &Q);CHKERRQ(ierr); 1631 ierr = PetscDualSpaceGetDimension(Q, &fNb);CHKERRQ(ierr); 1632 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 1633 } 1634 fTotDim += fNb*Nc; 1635 } 1636 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1637 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 1638 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 1639 PetscFE feC; 1640 PetscFV fvC; 1641 PetscDualSpace QF, QC; 1642 PetscInt NcF, NcC, fpdim, cpdim; 1643 1644 if (feRef[field]) { 1645 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 1646 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 1647 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 1648 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 1649 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1650 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 1651 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1652 } else { 1653 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &fvC);CHKERRQ(ierr); 1654 ierr = PetscFVGetNumComponents(fvC, &NcC);CHKERRQ(ierr); 1655 ierr = PetscFVGetNumComponents(fvRef[field], &NcF);CHKERRQ(ierr); 1656 ierr = PetscFVGetDualSpace(fvRef[field], &QF);CHKERRQ(ierr); 1657 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1658 ierr = PetscFVGetDualSpace(fvC, &QC);CHKERRQ(ierr); 1659 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1660 } 1661 if (NcF != NcC) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components in fine space field %d does not match coarse field %d", NcF, NcC); 1662 for (c = 0; c < cpdim; ++c) { 1663 PetscQuadrature cfunc; 1664 const PetscReal *cqpoints; 1665 PetscInt NpC; 1666 PetscBool found = PETSC_FALSE; 1667 1668 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 1669 ierr = PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 1670 if (NpC != 1 && feRef[field]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 1671 for (f = 0; f < fpdim; ++f) { 1672 PetscQuadrature ffunc; 1673 const PetscReal *fqpoints; 1674 PetscReal sum = 0.0; 1675 PetscInt NpF, comp; 1676 1677 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 1678 ierr = PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 1679 if (NpC != NpF) continue; 1680 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 1681 if (sum > 1.0e-9) continue; 1682 for (comp = 0; comp < NcC; ++comp) { 1683 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 1684 } 1685 found = PETSC_TRUE; 1686 break; 1687 } 1688 if (!found) { 1689 /* TODO We really want the average here, but some asshole put VecScatter in the interface */ 1690 if (fvRef[field]) { 1691 PetscInt comp; 1692 for (comp = 0; comp < NcC; ++comp) { 1693 cmap[(offsetC+c)*NcC+comp] = (offsetF+0)*NcF+comp; 1694 } 1695 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Could not locate matching functional for injection"); 1696 } 1697 } 1698 offsetC += cpdim*NcC; 1699 offsetF += fpdim*NcF; 1700 } 1701 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);ierr = PetscFVDestroy(&fvRef[f]);CHKERRQ(ierr);} 1702 ierr = PetscFree2(feRef,fvRef);CHKERRQ(ierr); 1703 1704 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 1705 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 1706 ierr = VecGetOwnershipRange(cv, &startC, &endC);CHKERRQ(ierr); 1707 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 1708 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 1709 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 1710 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 1711 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 1712 for (c = cStart; c < cEnd; ++c) { 1713 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 1714 for (d = 0; d < cTotDim; ++d) { 1715 if ((cellCIndices[d] < startC) || (cellCIndices[d] >= endC)) continue; 1716 if ((findices[cellCIndices[d]-startC] >= 0) && (findices[cellCIndices[d]-startC] != cellFIndices[cmap[d]])) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Coarse dof %d maps to both %d and %d", cindices[cellCIndices[d]-startC], findices[cellCIndices[d]-startC], cellFIndices[cmap[d]]); 1717 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 1718 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 1719 } 1720 } 1721 ierr = PetscFree(cmap);CHKERRQ(ierr); 1722 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 1723 1724 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 1725 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 1726 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 1727 ierr = ISDestroy(&cis);CHKERRQ(ierr); 1728 ierr = ISDestroy(&fis);CHKERRQ(ierr); 1729 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 1730 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 1731 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1732 PetscFunctionReturn(0); 1733 } 1734