1 #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 #include <petscfe.h> 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "DMPlexGetScale" 7 PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale) 8 { 9 DM_Plex *mesh = (DM_Plex*) dm->data; 10 11 PetscFunctionBegin; 12 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 13 PetscValidPointer(scale, 3); 14 *scale = mesh->scale[unit]; 15 PetscFunctionReturn(0); 16 } 17 18 #undef __FUNCT__ 19 #define __FUNCT__ "DMPlexSetScale" 20 PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale) 21 { 22 DM_Plex *mesh = (DM_Plex*) dm->data; 23 24 PetscFunctionBegin; 25 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 26 mesh->scale[unit] = scale; 27 PetscFunctionReturn(0); 28 } 29 30 PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k) 31 { 32 switch (i) { 33 case 0: 34 switch (j) { 35 case 0: return 0; 36 case 1: 37 switch (k) { 38 case 0: return 0; 39 case 1: return 0; 40 case 2: return 1; 41 } 42 case 2: 43 switch (k) { 44 case 0: return 0; 45 case 1: return -1; 46 case 2: return 0; 47 } 48 } 49 case 1: 50 switch (j) { 51 case 0: 52 switch (k) { 53 case 0: return 0; 54 case 1: return 0; 55 case 2: return -1; 56 } 57 case 1: return 0; 58 case 2: 59 switch (k) { 60 case 0: return 1; 61 case 1: return 0; 62 case 2: return 0; 63 } 64 } 65 case 2: 66 switch (j) { 67 case 0: 68 switch (k) { 69 case 0: return 0; 70 case 1: return 1; 71 case 2: return 0; 72 } 73 case 1: 74 switch (k) { 75 case 0: return -1; 76 case 1: return 0; 77 case 2: return 0; 78 } 79 case 2: return 0; 80 } 81 } 82 return 0; 83 } 84 85 #undef __FUNCT__ 86 #define __FUNCT__ "DMPlexCreateRigidBody" 87 /*@C 88 DMPlexCreateRigidBody - create rigid body modes from coordinates 89 90 Collective on DM 91 92 Input Arguments: 93 + dm - the DM 94 . section - the local section associated with the rigid field, or NULL for the default section 95 - globalSection - the global section associated with the rigid field, or NULL for the default section 96 97 Output Argument: 98 . sp - the null space 99 100 Note: This is necessary to take account of Dirichlet conditions on the displacements 101 102 Level: advanced 103 104 .seealso: MatNullSpaceCreate() 105 @*/ 106 PetscErrorCode DMPlexCreateRigidBody(DM dm, PetscSection section, PetscSection globalSection, MatNullSpace *sp) 107 { 108 MPI_Comm comm; 109 Vec coordinates, localMode, mode[6]; 110 PetscSection coordSection; 111 PetscScalar *coords; 112 PetscInt dim, vStart, vEnd, v, n, m, d, i, j; 113 PetscErrorCode ierr; 114 115 PetscFunctionBegin; 116 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 117 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 118 if (dim == 1) { 119 ierr = MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);CHKERRQ(ierr); 120 PetscFunctionReturn(0); 121 } 122 if (!section) {ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr);} 123 if (!globalSection) {ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr);} 124 ierr = PetscSectionGetConstrainedStorageSize(globalSection, &n);CHKERRQ(ierr); 125 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 126 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 127 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 128 m = (dim*(dim+1))/2; 129 ierr = VecCreate(comm, &mode[0]);CHKERRQ(ierr); 130 ierr = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHKERRQ(ierr); 131 ierr = VecSetUp(mode[0]);CHKERRQ(ierr); 132 for (i = 1; i < m; ++i) {ierr = VecDuplicate(mode[0], &mode[i]);CHKERRQ(ierr);} 133 /* Assume P1 */ 134 ierr = DMGetLocalVector(dm, &localMode);CHKERRQ(ierr); 135 for (d = 0; d < dim; ++d) { 136 PetscScalar values[3] = {0.0, 0.0, 0.0}; 137 138 values[d] = 1.0; 139 ierr = VecSet(localMode, 0.0);CHKERRQ(ierr); 140 for (v = vStart; v < vEnd; ++v) { 141 ierr = DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);CHKERRQ(ierr); 142 } 143 ierr = DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 144 ierr = DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 145 } 146 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr); 147 for (d = dim; d < dim*(dim+1)/2; ++d) { 148 PetscInt i, j, k = dim > 2 ? d - dim : d; 149 150 ierr = VecSet(localMode, 0.0);CHKERRQ(ierr); 151 for (v = vStart; v < vEnd; ++v) { 152 PetscScalar values[3] = {0.0, 0.0, 0.0}; 153 PetscInt off; 154 155 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 156 for (i = 0; i < dim; ++i) { 157 for (j = 0; j < dim; ++j) { 158 values[j] += epsilon(i, j, k)*PetscRealPart(coords[off+i]); 159 } 160 } 161 ierr = DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);CHKERRQ(ierr); 162 } 163 ierr = DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 164 ierr = DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);CHKERRQ(ierr); 165 } 166 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr); 167 ierr = DMRestoreLocalVector(dm, &localMode);CHKERRQ(ierr); 168 for (i = 0; i < dim; ++i) {ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr);} 169 /* Orthonormalize system */ 170 for (i = dim; i < m; ++i) { 171 PetscScalar dots[6]; 172 173 ierr = VecMDot(mode[i], i, mode, dots);CHKERRQ(ierr); 174 for (j = 0; j < i; ++j) dots[j] *= -1.0; 175 ierr = VecMAXPY(mode[i], i, dots, mode);CHKERRQ(ierr); 176 ierr = VecNormalize(mode[i], NULL);CHKERRQ(ierr); 177 } 178 ierr = MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);CHKERRQ(ierr); 179 for (i = 0; i< m; ++i) {ierr = VecDestroy(&mode[i]);CHKERRQ(ierr);} 180 PetscFunctionReturn(0); 181 } 182 183 #undef __FUNCT__ 184 #define __FUNCT__ "DMPlexProjectFunctionLabelLocal" 185 PetscErrorCode DMPlexProjectFunctionLabelLocal(DM dm, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 186 { 187 PetscDualSpace *sp; 188 PetscSection section; 189 PetscScalar *values; 190 PetscReal *v0, *J, detJ; 191 PetscBool *fieldActive; 192 PetscInt numFields, numComp, dim, spDim, totDim = 0, numValues, cStart, cEnd, f, d, v, i, comp; 193 PetscErrorCode ierr; 194 195 PetscFunctionBegin; 196 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 197 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 198 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 199 ierr = PetscMalloc3(numFields,&sp,dim,&v0,dim*dim,&J);CHKERRQ(ierr); 200 for (f = 0; f < numFields; ++f) { 201 ierr = PetscFEGetDualSpace(fe[f], &sp[f]);CHKERRQ(ierr); 202 ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr); 203 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 204 totDim += spDim*numComp; 205 } 206 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 207 ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 208 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 209 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 210 ierr = DMGetWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 211 for (f = 0; f < numFields; ++f) fieldActive[f] = funcs[f] ? PETSC_TRUE : PETSC_FALSE; 212 for (i = 0; i < numIds; ++i) { 213 IS pointIS; 214 const PetscInt *points; 215 PetscInt n, p; 216 217 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 218 ierr = ISGetLocalSize(pointIS, &n);CHKERRQ(ierr); 219 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 220 for (p = 0; p < n; ++p) { 221 const PetscInt point = points[p]; 222 PetscCellGeometry geom; 223 224 if ((point < cStart) || (point >= cEnd)) continue; 225 ierr = DMPlexComputeCellGeometry(dm, point, v0, J, NULL, &detJ);CHKERRQ(ierr); 226 geom.v0 = v0; 227 geom.J = J; 228 geom.detJ = &detJ; 229 for (f = 0, v = 0; f < numFields; ++f) { 230 void * const ctx = ctxs ? ctxs[f] : NULL; 231 ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr); 232 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 233 for (d = 0; d < spDim; ++d) { 234 if (funcs[f]) { 235 ierr = PetscDualSpaceApply(sp[f], d, geom, numComp, funcs[f], ctx, &values[v]);CHKERRQ(ierr); 236 } else { 237 for (comp = 0; comp < numComp; ++comp) values[v+comp] = 0.0; 238 } 239 v += numComp; 240 } 241 } 242 ierr = DMPlexVecSetFieldClosure_Internal(dm, section, localX, fieldActive, point, values, mode);CHKERRQ(ierr); 243 } 244 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 245 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 246 } 247 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 248 ierr = DMRestoreWorkArray(dm, numFields, PETSC_BOOL, &fieldActive);CHKERRQ(ierr); 249 ierr = PetscFree3(sp,v0,J);CHKERRQ(ierr); 250 PetscFunctionReturn(0); 251 } 252 253 #undef __FUNCT__ 254 #define __FUNCT__ "DMPlexProjectFunctionLocal" 255 PetscErrorCode DMPlexProjectFunctionLocal(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 256 { 257 PetscDualSpace *sp; 258 PetscSection section; 259 PetscScalar *values; 260 PetscReal *v0, *J, detJ; 261 PetscInt numFields, numComp, dim, spDim, totDim = 0, numValues, cStart, cEnd, c, f, d, v, comp; 262 PetscErrorCode ierr; 263 264 PetscFunctionBegin; 265 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 266 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 267 ierr = PetscMalloc1(numFields, &sp);CHKERRQ(ierr); 268 for (f = 0; f < numFields; ++f) { 269 ierr = PetscFEGetDualSpace(fe[f], &sp[f]);CHKERRQ(ierr); 270 ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr); 271 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 272 totDim += spDim*numComp; 273 } 274 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 275 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 276 ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 277 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 278 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 279 ierr = PetscMalloc2(dim,&v0,dim*dim,&J);CHKERRQ(ierr); 280 for (c = cStart; c < cEnd; ++c) { 281 PetscCellGeometry geom; 282 283 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);CHKERRQ(ierr); 284 geom.v0 = v0; 285 geom.J = J; 286 geom.detJ = &detJ; 287 for (f = 0, v = 0; f < numFields; ++f) { 288 void * const ctx = ctxs ? ctxs[f] : NULL; 289 ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr); 290 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 291 for (d = 0; d < spDim; ++d) { 292 if (funcs[f]) { 293 ierr = PetscDualSpaceApply(sp[f], d, geom, numComp, funcs[f], ctx, &values[v]);CHKERRQ(ierr); 294 } else { 295 for (comp = 0; comp < numComp; ++comp) values[v+comp] = 0.0; 296 } 297 v += numComp; 298 } 299 } 300 ierr = DMPlexVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr); 301 } 302 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 303 ierr = PetscFree2(v0,J);CHKERRQ(ierr); 304 ierr = PetscFree(sp);CHKERRQ(ierr); 305 PetscFunctionReturn(0); 306 } 307 308 #undef __FUNCT__ 309 #define __FUNCT__ "DMPlexProjectFunction" 310 /*@C 311 DMPlexProjectFunction - This projects the given function into the function space provided. 312 313 Input Parameters: 314 + dm - The DM 315 . fe - The PetscFE associated with the field 316 . funcs - The coordinate functions to evaluate, one per field 317 . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null. 318 - mode - The insertion mode for values 319 320 Output Parameter: 321 . X - vector 322 323 Level: developer 324 325 .seealso: DMPlexComputeL2Diff() 326 @*/ 327 PetscErrorCode DMPlexProjectFunction(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X) 328 { 329 Vec localX; 330 PetscErrorCode ierr; 331 332 PetscFunctionBegin; 333 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 334 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 335 ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, ctxs, mode, localX);CHKERRQ(ierr); 336 ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr); 337 ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr); 338 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 339 PetscFunctionReturn(0); 340 } 341 342 #undef __FUNCT__ 343 #define __FUNCT__ "DMPlexInsertBoundaryValuesFEM" 344 PetscErrorCode DMPlexInsertBoundaryValuesFEM(DM dm, Vec localX) 345 { 346 void (**funcs)(const PetscReal x[], PetscScalar *u, void *ctx); 347 void **ctxs; 348 PetscFE *fe; 349 PetscInt numFields, f, numBd, b; 350 PetscErrorCode ierr; 351 352 PetscFunctionBegin; 353 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 354 PetscValidHeaderSpecific(localX, VEC_CLASSID, 2); 355 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 356 ierr = PetscMalloc3(numFields,&fe,numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 357 for (f = 0; f < numFields; ++f) {ierr = DMGetField(dm, f, (PetscObject *) &fe[f]);CHKERRQ(ierr);} 358 /* OPT: Could attempt to do multiple BCs at once */ 359 ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 360 for (b = 0; b < numBd; ++b) { 361 DMLabel label; 362 const PetscInt *ids; 363 const char *name; 364 PetscInt numids, field; 365 PetscBool isEssential; 366 void (*func)(); 367 void *ctx; 368 369 /* TODO: We need to set only the part indicated by the ids */ 370 ierr = DMPlexGetBoundary(dm, b, &isEssential, &name, &field, &func, &numids, &ids, &ctx);CHKERRQ(ierr); 371 ierr = DMPlexGetLabel(dm, name, &label);CHKERRQ(ierr); 372 for (f = 0; f < numFields; ++f) { 373 funcs[f] = field == f ? (void (*)(const PetscReal[], PetscScalar *, void *)) func : NULL; 374 ctxs[f] = field == f ? ctx : NULL; 375 } 376 ierr = DMPlexProjectFunctionLabelLocal(dm, label, numids, ids, fe, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 377 } 378 ierr = PetscFree3(fe,funcs,ctxs);CHKERRQ(ierr); 379 PetscFunctionReturn(0); 380 } 381 382 /* Assuming dim == 3 */ 383 typedef struct { 384 PetscScalar normal[3]; /* Area-scaled normals */ 385 PetscScalar centroid[3]; /* Location of centroid (quadrature point) */ 386 PetscScalar grad[2][3]; /* Face contribution to gradient in left and right cell */ 387 } FaceGeom; 388 389 #undef __FUNCT__ 390 #define __FUNCT__ "DMPlexInsertBoundaryValuesFVM" 391 PetscErrorCode DMPlexInsertBoundaryValuesFVM(DM dm, PetscReal time, Vec locX) 392 { 393 DM dmFace; 394 Vec faceGeometry; 395 DMLabel label; 396 const PetscScalar *facegeom; 397 PetscScalar *x; 398 PetscInt numBd, b, fStart, fEnd; 399 PetscErrorCode ierr; 400 401 PetscFunctionBegin; 402 /* TODO Pull this geometry calculation into the library */ 403 ierr = PetscObjectQuery((PetscObject) dm, "FaceGeometry", (PetscObject *) &faceGeometry);CHKERRQ(ierr); 404 ierr = DMPlexGetLabel(dm, "Face Sets", &label);CHKERRQ(ierr); 405 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 406 ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 407 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 408 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 409 ierr = VecGetArray(locX, &x);CHKERRQ(ierr); 410 for (b = 0; b < numBd; ++b) { 411 PetscErrorCode (*func)(PetscReal,const PetscScalar*,const PetscScalar*,const PetscScalar*,PetscScalar*,void*); 412 const PetscInt *ids; 413 PetscInt numids, i; 414 void *ctx; 415 416 ierr = DMPlexGetBoundary(dm, b, NULL, NULL, NULL, (void (**)()) &func, &numids, &ids, &ctx);CHKERRQ(ierr); 417 for (i = 0; i < numids; ++i) { 418 IS faceIS; 419 const PetscInt *faces; 420 PetscInt numFaces, f; 421 422 ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr); 423 if (!faceIS) continue; /* No points with that id on this process */ 424 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 425 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 426 for (f = 0; f < numFaces; ++f) { 427 const PetscInt face = faces[f], *cells; 428 const PetscScalar *xI; 429 PetscScalar *xG; 430 const FaceGeom *fg; 431 432 if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */ 433 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 434 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 435 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr); 436 ierr = DMPlexPointLocalRef(dm, cells[1], x, &xG);CHKERRQ(ierr); 437 ierr = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);CHKERRQ(ierr); 438 } 439 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 440 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 441 } 442 } 443 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 444 ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr); 445 PetscFunctionReturn(0); 446 } 447 448 #undef __FUNCT__ 449 #define __FUNCT__ "DMPlexComputeL2Diff" 450 /*@C 451 DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h. 452 453 Input Parameters: 454 + dm - The DM 455 . fe - The PetscFE object for each field 456 . funcs - The functions to evaluate for each field component 457 . ctxs - Optional array of contexts to pass to each function, or NULL. 458 - X - The coefficient vector u_h 459 460 Output Parameter: 461 . diff - The diff ||u - u_h||_2 462 463 Level: developer 464 465 .seealso: DMPlexProjectFunction(), DMPlexComputeL2GradientDiff() 466 @*/ 467 PetscErrorCode DMPlexComputeL2Diff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 468 { 469 const PetscInt debug = 0; 470 PetscSection section; 471 PetscQuadrature quad; 472 Vec localX; 473 PetscScalar *funcVal; 474 PetscReal *coords, *v0, *J, *invJ, detJ; 475 PetscReal localDiff = 0.0; 476 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp; 477 PetscErrorCode ierr; 478 479 PetscFunctionBegin; 480 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 481 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 482 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 483 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 484 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 485 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 486 for (field = 0; field < numFields; ++field) { 487 PetscInt Nc; 488 489 ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr); 490 numComponents += Nc; 491 } 492 ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 493 ierr = PetscMalloc5(numComponents,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 494 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 495 ierr = PetscFEGetQuadrature(fe[0], &quad);CHKERRQ(ierr); 496 for (c = cStart; c < cEnd; ++c) { 497 PetscScalar *x = NULL; 498 PetscReal elemDiff = 0.0; 499 500 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr); 501 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 502 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 503 504 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 505 void * const ctx = ctxs ? ctxs[field] : NULL; 506 const PetscReal *quadPoints, *quadWeights; 507 PetscReal *basis; 508 PetscInt numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f; 509 510 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 511 ierr = PetscFEGetDimension(fe[field], &numBasisFuncs);CHKERRQ(ierr); 512 ierr = PetscFEGetNumComponents(fe[field], &numBasisComps);CHKERRQ(ierr); 513 ierr = PetscFEGetDefaultTabulation(fe[field], &basis, NULL, NULL);CHKERRQ(ierr); 514 if (debug) { 515 char title[1024]; 516 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 517 ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr); 518 } 519 for (q = 0; q < numQuadPoints; ++q) { 520 for (d = 0; d < dim; d++) { 521 coords[d] = v0[d]; 522 for (e = 0; e < dim; e++) { 523 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 524 } 525 } 526 (*funcs[field])(coords, funcVal, ctx); 527 for (fc = 0; fc < numBasisComps; ++fc) { 528 PetscScalar interpolant = 0.0; 529 530 for (f = 0; f < numBasisFuncs; ++f) { 531 const PetscInt fidx = f*numBasisComps+fc; 532 interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx]; 533 } 534 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);} 535 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 536 } 537 } 538 comp += numBasisComps; 539 fieldOffset += numBasisFuncs*numBasisComps; 540 } 541 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 542 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 543 localDiff += elemDiff; 544 } 545 ierr = PetscFree5(funcVal,coords,v0,J,invJ);CHKERRQ(ierr); 546 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 547 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 548 *diff = PetscSqrtReal(*diff); 549 PetscFunctionReturn(0); 550 } 551 552 #undef __FUNCT__ 553 #define __FUNCT__ "DMPlexComputeL2GradientDiff" 554 /*@C 555 DMPlexComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h. 556 557 Input Parameters: 558 + dm - The DM 559 . fe - The PetscFE object for each field 560 . funcs - The gradient functions to evaluate for each field component 561 . ctxs - Optional array of contexts to pass to each function, or NULL. 562 . X - The coefficient vector u_h 563 - n - The vector to project along 564 565 Output Parameter: 566 . diff - The diff ||(grad u - grad u_h) . n||_2 567 568 Level: developer 569 570 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff() 571 @*/ 572 PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff) 573 { 574 const PetscInt debug = 0; 575 PetscSection section; 576 PetscQuadrature quad; 577 Vec localX; 578 PetscScalar *funcVal, *interpolantVec; 579 PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ; 580 PetscReal localDiff = 0.0; 581 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp; 582 PetscErrorCode ierr; 583 584 PetscFunctionBegin; 585 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 586 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 587 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 588 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 589 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 590 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 591 for (field = 0; field < numFields; ++field) { 592 PetscInt Nc; 593 594 ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr); 595 numComponents += Nc; 596 } 597 /* ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 598 ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr); 599 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 600 ierr = PetscFEGetQuadrature(fe[0], &quad);CHKERRQ(ierr); 601 for (c = cStart; c < cEnd; ++c) { 602 PetscScalar *x = NULL; 603 PetscReal elemDiff = 0.0; 604 605 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr); 606 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 607 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 608 609 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 610 void * const ctx = ctxs ? ctxs[field] : NULL; 611 const PetscReal *quadPoints, *quadWeights; 612 PetscReal *basisDer; 613 PetscInt numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g; 614 615 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 616 ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr); 617 ierr = PetscFEGetNumComponents(fe[field], &Ncomp);CHKERRQ(ierr); 618 ierr = PetscFEGetDefaultTabulation(fe[field], NULL, &basisDer, NULL);CHKERRQ(ierr); 619 if (debug) { 620 char title[1024]; 621 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 622 ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr); 623 } 624 for (q = 0; q < numQuadPoints; ++q) { 625 for (d = 0; d < dim; d++) { 626 coords[d] = v0[d]; 627 for (e = 0; e < dim; e++) { 628 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 629 } 630 } 631 (*funcs[field])(coords, n, funcVal, ctx); 632 for (fc = 0; fc < Ncomp; ++fc) { 633 PetscScalar interpolant = 0.0; 634 635 for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0; 636 for (f = 0; f < Nb; ++f) { 637 const PetscInt fidx = f*Ncomp+fc; 638 639 for (d = 0; d < dim; ++d) { 640 realSpaceDer[d] = 0.0; 641 for (g = 0; g < dim; ++g) { 642 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g]; 643 } 644 interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d]; 645 } 646 } 647 for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d]; 648 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d fieldDer %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);} 649 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 650 } 651 } 652 comp += Ncomp; 653 fieldOffset += Nb*Ncomp; 654 } 655 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 656 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 657 localDiff += elemDiff; 658 } 659 ierr = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr); 660 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 661 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 662 *diff = PetscSqrtReal(*diff); 663 PetscFunctionReturn(0); 664 } 665 666 #undef __FUNCT__ 667 #define __FUNCT__ "DMPlexComputeL2FieldDiff" 668 PetscErrorCode DMPlexComputeL2FieldDiff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[]) 669 { 670 const PetscInt debug = 0; 671 PetscSection section; 672 PetscQuadrature quad; 673 Vec localX; 674 PetscScalar *funcVal; 675 PetscReal *coords, *v0, *J, *invJ, detJ; 676 PetscReal *localDiff; 677 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp; 678 PetscErrorCode ierr; 679 680 PetscFunctionBegin; 681 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 682 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 683 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 684 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 685 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 686 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 687 for (field = 0; field < numFields; ++field) { 688 PetscInt Nc; 689 690 ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr); 691 numComponents += Nc; 692 } 693 ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 694 ierr = PetscCalloc6(numFields,&localDiff,numComponents,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 695 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 696 ierr = PetscFEGetQuadrature(fe[0], &quad);CHKERRQ(ierr); 697 for (c = cStart; c < cEnd; ++c) { 698 PetscScalar *x = NULL; 699 700 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr); 701 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 702 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 703 704 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 705 void * const ctx = ctxs ? ctxs[field] : NULL; 706 const PetscReal *quadPoints, *quadWeights; 707 PetscReal *basis, elemDiff = 0.0; 708 PetscInt numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f; 709 710 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 711 ierr = PetscFEGetDimension(fe[field], &numBasisFuncs);CHKERRQ(ierr); 712 ierr = PetscFEGetNumComponents(fe[field], &numBasisComps);CHKERRQ(ierr); 713 ierr = PetscFEGetDefaultTabulation(fe[field], &basis, NULL, NULL);CHKERRQ(ierr); 714 if (debug) { 715 char title[1024]; 716 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 717 ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr); 718 } 719 for (q = 0; q < numQuadPoints; ++q) { 720 for (d = 0; d < dim; d++) { 721 coords[d] = v0[d]; 722 for (e = 0; e < dim; e++) { 723 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 724 } 725 } 726 (*funcs[field])(coords, funcVal, ctx); 727 for (fc = 0; fc < numBasisComps; ++fc) { 728 PetscScalar interpolant = 0.0; 729 730 for (f = 0; f < numBasisFuncs; ++f) { 731 const PetscInt fidx = f*numBasisComps+fc; 732 interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx]; 733 } 734 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ);CHKERRQ(ierr);} 735 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 736 } 737 } 738 comp += numBasisComps; 739 fieldOffset += numBasisFuncs*numBasisComps; 740 localDiff[field] += elemDiff; 741 } 742 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 743 } 744 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 745 ierr = MPI_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 746 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 747 ierr = PetscFree6(localDiff,funcVal,coords,v0,J,invJ);CHKERRQ(ierr); 748 PetscFunctionReturn(0); 749 } 750 751 #undef __FUNCT__ 752 #define __FUNCT__ "DMPlexComputeIntegralFEM" 753 /*@ 754 DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user 755 756 Input Parameters: 757 + dm - The mesh 758 . X - Local input vector 759 . obj - The functions to integrate for each field 760 - user - The user context 761 762 Output Parameter: 763 . integral - Local integral for each field 764 765 Note: 766 The first member of the user context must be an FEMContext. 767 768 Level: developer 769 770 .seealso: DMPlexComputeResidualFEM() 771 @*/ 772 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, void (**obj)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscReal *integral, void *user) 773 { 774 DM_Plex *mesh = (DM_Plex *) dm->data; 775 PetscFEM *fem = (PetscFEM *) user; 776 PetscFE *fe = fem->fe; 777 PetscFE *feAux = fem->feAux; 778 DM dmAux; 779 Vec localX, A; 780 PetscQuadrature q; 781 PetscCellGeometry geom; 782 PetscSection section, sectionAux; 783 PetscReal *v0, *J, *invJ, *detJ; 784 PetscScalar *u, *a = NULL; 785 PetscInt dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, c; 786 PetscInt cellDof = 0, numComponents = 0; 787 PetscInt cellDofAux = 0, numComponentsAux = 0; 788 PetscErrorCode ierr; 789 790 PetscFunctionBegin; 791 /*ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr);*/ 792 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 793 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 794 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 795 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 796 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 797 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 798 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 799 numCells = cEnd - cStart; 800 for (f = 0; f < Nf; ++f) { 801 PetscInt Nb, Nc; 802 803 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 804 ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr); 805 cellDof += Nb*Nc; 806 numComponents += Nc; 807 integral[f] = 0.0; 808 } 809 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 810 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 811 if (dmAux) { 812 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 813 ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr); 814 } 815 for (f = 0; f < NfAux; ++f) { 816 PetscInt Nb, Nc; 817 818 ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr); 819 ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr); 820 cellDofAux += Nb*Nc; 821 numComponentsAux += Nc; 822 } 823 ierr = DMPlexInsertBoundaryValuesFEM(dm, localX);CHKERRQ(ierr); 824 ierr = PetscMalloc5(numCells*cellDof,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ);CHKERRQ(ierr); 825 if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);} 826 for (c = cStart; c < cEnd; ++c) { 827 PetscScalar *x = NULL; 828 PetscInt i; 829 830 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 831 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 832 ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 833 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 834 ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 835 if (dmAux) { 836 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 837 for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x[i]; 838 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 839 } 840 } 841 for (f = 0; f < Nf; ++f) { 842 PetscInt numQuadPoints, Nb; 843 /* Conforming batches */ 844 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 845 /* Remainder */ 846 PetscInt Nr, offset; 847 848 ierr = PetscFEGetQuadrature(fe[f], &q);CHKERRQ(ierr); 849 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 850 ierr = PetscFEGetTileSizes(fe[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 851 ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 852 blockSize = Nb*numQuadPoints; 853 batchSize = numBlocks * blockSize; 854 ierr = PetscFESetTileSizes(fe[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 855 numChunks = numCells / (numBatches*batchSize); 856 Ne = numChunks*numBatches*batchSize; 857 Nr = numCells % (numBatches*batchSize); 858 offset = numCells - Nr; 859 geom.v0 = v0; 860 geom.J = J; 861 geom.invJ = invJ; 862 geom.detJ = detJ; 863 ierr = PetscFEIntegrate(fe[f], Ne, Nf, fe, f, geom, u, NfAux, feAux, a, obj[f], integral);CHKERRQ(ierr); 864 geom.v0 = &v0[offset*dim]; 865 geom.J = &J[offset*dim*dim]; 866 geom.invJ = &invJ[offset*dim*dim]; 867 geom.detJ = &detJ[offset]; 868 ierr = PetscFEIntegrate(fe[f], Nr, Nf, fe, f, geom, &u[offset*cellDof], NfAux, feAux, &a[offset*cellDofAux], obj[f], integral);CHKERRQ(ierr); 869 } 870 ierr = PetscFree5(u,v0,J,invJ,detJ);CHKERRQ(ierr); 871 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 872 if (mesh->printFEM) { 873 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr); 874 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", integral[f]);CHKERRQ(ierr);} 875 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 876 } 877 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 878 /* TODO: Allreduce for integral */ 879 /*ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr);*/ 880 PetscFunctionReturn(0); 881 } 882 883 #undef __FUNCT__ 884 #define __FUNCT__ "DMPlexComputeResidualFEM" 885 /*@ 886 DMPlexComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user 887 888 Input Parameters: 889 + dm - The mesh 890 . X - Local input vector 891 - user - The user context 892 893 Output Parameter: 894 . F - Local output vector 895 896 Note: 897 The first member of the user context must be an FEMContext. 898 899 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 900 like a GPU, or vectorize on a multicore machine. 901 902 Level: developer 903 904 .seealso: DMPlexComputeJacobianActionFEM() 905 @*/ 906 PetscErrorCode DMPlexComputeResidualFEM(DM dm, Vec X, Vec F, void *user) 907 { 908 DM_Plex *mesh = (DM_Plex *) dm->data; 909 PetscFEM *fem = (PetscFEM *) user; 910 PetscFE *fe = fem->fe; 911 PetscFE *feAux = fem->feAux; 912 PetscFE *feBd = fem->feBd; 913 const char *name = "Residual"; 914 DM dmAux; 915 Vec A; 916 PetscQuadrature q; 917 PetscCellGeometry geom; 918 PetscSection section, sectionAux; 919 PetscReal *v0, *J, *invJ, *detJ; 920 PetscScalar *elemVec, *u, *a = NULL; 921 PetscInt dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, c; 922 PetscInt cellDof = 0, numComponents = 0; 923 PetscInt cellDofAux = 0, numComponentsAux = 0; 924 PetscErrorCode ierr; 925 926 PetscFunctionBegin; 927 ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr); 928 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 929 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 930 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 931 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 932 numCells = cEnd - cStart; 933 for (f = 0; f < Nf; ++f) { 934 PetscInt Nb, Nc; 935 936 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 937 ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr); 938 cellDof += Nb*Nc; 939 numComponents += Nc; 940 } 941 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 942 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 943 if (dmAux) { 944 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 945 ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr); 946 } 947 for (f = 0; f < NfAux; ++f) { 948 PetscInt Nb, Nc; 949 950 ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr); 951 ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr); 952 cellDofAux += Nb*Nc; 953 numComponentsAux += Nc; 954 } 955 ierr = DMPlexInsertBoundaryValuesFEM(dm, X);CHKERRQ(ierr); 956 ierr = VecSet(F, 0.0);CHKERRQ(ierr); 957 ierr = PetscMalloc6(numCells*cellDof,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof,&elemVec);CHKERRQ(ierr); 958 if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);} 959 for (c = cStart; c < cEnd; ++c) { 960 PetscScalar *x = NULL; 961 PetscInt i; 962 963 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 964 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 965 ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 966 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 967 ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 968 if (dmAux) { 969 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 970 for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x[i]; 971 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 972 } 973 } 974 for (f = 0; f < Nf; ++f) { 975 void (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0Funcs[f]; 976 void (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1Funcs[f]; 977 PetscInt numQuadPoints, Nb; 978 /* Conforming batches */ 979 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 980 /* Remainder */ 981 PetscInt Nr, offset; 982 983 ierr = PetscFEGetQuadrature(fe[f], &q);CHKERRQ(ierr); 984 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 985 ierr = PetscFEGetTileSizes(fe[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 986 ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 987 blockSize = Nb*numQuadPoints; 988 batchSize = numBlocks * blockSize; 989 ierr = PetscFESetTileSizes(fe[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 990 numChunks = numCells / (numBatches*batchSize); 991 Ne = numChunks*numBatches*batchSize; 992 Nr = numCells % (numBatches*batchSize); 993 offset = numCells - Nr; 994 geom.v0 = v0; 995 geom.J = J; 996 geom.invJ = invJ; 997 geom.detJ = detJ; 998 ierr = PetscFEIntegrateResidual(fe[f], Ne, Nf, fe, f, geom, u, NfAux, feAux, a, f0, f1, elemVec);CHKERRQ(ierr); 999 geom.v0 = &v0[offset*dim]; 1000 geom.J = &J[offset*dim*dim]; 1001 geom.invJ = &invJ[offset*dim*dim]; 1002 geom.detJ = &detJ[offset]; 1003 ierr = PetscFEIntegrateResidual(fe[f], Nr, Nf, fe, f, geom, &u[offset*cellDof], NfAux, feAux, &a[offset*cellDofAux], f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr); 1004 } 1005 for (c = cStart; c < cEnd; ++c) { 1006 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);} 1007 ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr); 1008 } 1009 ierr = PetscFree6(u,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr); 1010 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1011 if (feBd) { 1012 DMLabel depth; 1013 PetscInt numBd, bd; 1014 1015 for (f = 0, cellDof = 0, numComponents = 0; f < Nf; ++f) { 1016 PetscInt Nb, Nc; 1017 1018 ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr); 1019 ierr = PetscFEGetNumComponents(feBd[f], &Nc);CHKERRQ(ierr); 1020 cellDof += Nb*Nc; 1021 numComponents += Nc; 1022 } 1023 ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr); 1024 ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 1025 for (bd = 0; bd < numBd; ++bd) { 1026 const char *bdLabel; 1027 DMLabel label; 1028 IS pointIS; 1029 const PetscInt *points; 1030 const PetscInt *values; 1031 PetscReal *n; 1032 PetscInt field, numValues, numPoints, p, dep, numFaces; 1033 PetscBool isEssential; 1034 1035 ierr = DMPlexGetBoundary(dm, bd, &isEssential, &bdLabel, &field, NULL, &numValues, &values, NULL);CHKERRQ(ierr); 1036 if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this"); 1037 ierr = DMPlexGetLabel(dm, bdLabel, &label);CHKERRQ(ierr); 1038 ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr); 1039 ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr); 1040 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 1041 for (p = 0, numFaces = 0; p < numPoints; ++p) { 1042 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 1043 if (dep == dim-1) ++numFaces; 1044 } 1045 ierr = PetscMalloc7(numFaces*cellDof,&u,numFaces*dim,&v0,numFaces*dim,&n,numFaces*dim*dim,&J,numFaces*dim*dim,&invJ,numFaces,&detJ,numFaces*cellDof,&elemVec);CHKERRQ(ierr); 1046 for (p = 0, f = 0; p < numPoints; ++p) { 1047 const PetscInt point = points[p]; 1048 PetscScalar *x = NULL; 1049 PetscInt i; 1050 1051 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 1052 if (dep != dim-1) continue; 1053 ierr = DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr); 1054 ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]); 1055 if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point); 1056 ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 1057 for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i]; 1058 ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 1059 ++f; 1060 } 1061 for (f = 0; f < Nf; ++f) { 1062 void (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdFuncs[f]; 1063 void (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdFuncs[f]; 1064 PetscInt numQuadPoints, Nb; 1065 /* Conforming batches */ 1066 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 1067 /* Remainder */ 1068 PetscInt Nr, offset; 1069 1070 ierr = PetscFEGetQuadrature(feBd[f], &q);CHKERRQ(ierr); 1071 ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr); 1072 ierr = PetscFEGetTileSizes(feBd[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1073 ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 1074 blockSize = Nb*numQuadPoints; 1075 batchSize = numBlocks * blockSize; 1076 ierr = PetscFESetTileSizes(feBd[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1077 numChunks = numFaces / (numBatches*batchSize); 1078 Ne = numChunks*numBatches*batchSize; 1079 Nr = numFaces % (numBatches*batchSize); 1080 offset = numFaces - Nr; 1081 geom.v0 = v0; 1082 geom.n = n; 1083 geom.J = J; 1084 geom.invJ = invJ; 1085 geom.detJ = detJ; 1086 ierr = PetscFEIntegrateBdResidual(feBd[f], Ne, Nf, feBd, f, geom, u, 0, NULL, NULL, f0, f1, elemVec);CHKERRQ(ierr); 1087 geom.v0 = &v0[offset*dim]; 1088 geom.n = &n[offset*dim]; 1089 geom.J = &J[offset*dim*dim]; 1090 geom.invJ = &invJ[offset*dim*dim]; 1091 geom.detJ = &detJ[offset]; 1092 ierr = PetscFEIntegrateBdResidual(feBd[f], Nr, Nf, feBd, f, geom, &u[offset*cellDof], 0, NULL, NULL, f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr); 1093 } 1094 for (p = 0, f = 0; p < numPoints; ++p) { 1095 const PetscInt point = points[p]; 1096 1097 ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr); 1098 if (dep != dim-1) continue; 1099 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "BdResidual", cellDof, &elemVec[f*cellDof]);CHKERRQ(ierr);} 1100 ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*cellDof], ADD_VALUES);CHKERRQ(ierr); 1101 ++f; 1102 } 1103 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 1104 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 1105 ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr); 1106 } 1107 } 1108 if (mesh->printFEM) {ierr = DMPrintLocalVec(dm, name, mesh->printTol, F);CHKERRQ(ierr);} 1109 ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr); 1110 PetscFunctionReturn(0); 1111 } 1112 1113 #undef __FUNCT__ 1114 #define __FUNCT__ "DMPlexComputeIFunctionFEM" 1115 /*@ 1116 DMPlexComputeIFunctionFEM - Form the local implicit function F from the local input X, X_t using pointwise functions specified by the user 1117 1118 Input Parameters: 1119 + dm - The mesh 1120 . time - The current time 1121 . X - Local input vector 1122 . X_t - Time derivative of the local input vector 1123 - user - The user context 1124 1125 Output Parameter: 1126 . F - Local output vector 1127 1128 Note: 1129 The first member of the user context must be an FEMContext. 1130 1131 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 1132 like a GPU, or vectorize on a multicore machine. 1133 1134 Level: developer 1135 1136 .seealso: DMPlexComputeResidualFEM() 1137 @*/ 1138 PetscErrorCode DMPlexComputeIFunctionFEM(DM dm, PetscReal time, Vec X, Vec X_t, Vec F, void *user) 1139 { 1140 DM_Plex *mesh = (DM_Plex *) dm->data; 1141 PetscFEM *fem = (PetscFEM *) user; 1142 PetscFE *fe = fem->fe; 1143 PetscFE *feAux = fem->feAux; 1144 PetscFE *feBd = fem->feBd; 1145 const char *name = "Residual"; 1146 DM dmAux; 1147 Vec A; 1148 PetscQuadrature q; 1149 PetscCellGeometry geom; 1150 PetscSection section, sectionAux; 1151 PetscReal *v0, *J, *invJ, *detJ; 1152 PetscScalar *elemVec, *u, *u_t, *a = NULL; 1153 PetscInt dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, c; 1154 PetscInt cellDof = 0, numComponents = 0; 1155 PetscInt cellDofAux = 0, numComponentsAux = 0; 1156 PetscErrorCode ierr; 1157 1158 PetscFunctionBegin; 1159 ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr); 1160 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1161 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1162 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1163 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1164 numCells = cEnd - cStart; 1165 for (f = 0; f < Nf; ++f) { 1166 PetscInt Nb, Nc; 1167 1168 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 1169 ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr); 1170 cellDof += Nb*Nc; 1171 numComponents += Nc; 1172 } 1173 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1174 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1175 if (dmAux) { 1176 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1177 ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr); 1178 } 1179 for (f = 0; f < NfAux; ++f) { 1180 PetscInt Nb, Nc; 1181 1182 ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr); 1183 ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr); 1184 cellDofAux += Nb*Nc; 1185 numComponentsAux += Nc; 1186 } 1187 ierr = DMPlexInsertBoundaryValuesFEM(dm, X);CHKERRQ(ierr); 1188 ierr = VecSet(F, 0.0);CHKERRQ(ierr); 1189 ierr = PetscMalloc7(numCells*cellDof,&u,numCells*cellDof,&u_t,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof,&elemVec);CHKERRQ(ierr); 1190 if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);} 1191 for (c = cStart; c < cEnd; ++c) { 1192 PetscScalar *x = NULL, *x_t = NULL; 1193 PetscInt i; 1194 1195 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 1196 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 1197 ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 1198 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 1199 ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 1200 ierr = DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr); 1201 for (i = 0; i < cellDof; ++i) u_t[c*cellDof+i] = x_t[i]; 1202 ierr = DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr); 1203 if (dmAux) { 1204 PetscScalar *x_a = NULL; 1205 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x_a);CHKERRQ(ierr); 1206 for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x_a[i]; 1207 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x_a);CHKERRQ(ierr); 1208 } 1209 } 1210 for (f = 0; f < Nf; ++f) { 1211 void (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0IFuncs[f]; 1212 void (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1IFuncs[f]; 1213 PetscInt numQuadPoints, Nb; 1214 /* Conforming batches */ 1215 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 1216 /* Remainder */ 1217 PetscInt Nr, offset; 1218 1219 ierr = PetscFEGetQuadrature(fe[f], &q);CHKERRQ(ierr); 1220 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 1221 ierr = PetscFEGetTileSizes(fe[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1222 ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 1223 blockSize = Nb*numQuadPoints; 1224 batchSize = numBlocks * blockSize; 1225 ierr = PetscFESetTileSizes(fe[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1226 numChunks = numCells / (numBatches*batchSize); 1227 Ne = numChunks*numBatches*batchSize; 1228 Nr = numCells % (numBatches*batchSize); 1229 offset = numCells - Nr; 1230 geom.v0 = v0; 1231 geom.J = J; 1232 geom.invJ = invJ; 1233 geom.detJ = detJ; 1234 ierr = PetscFEIntegrateIFunction(fe[f], Ne, Nf, fe, f, geom, u, u_t, NfAux, feAux, a, f0, f1, elemVec);CHKERRQ(ierr); 1235 geom.v0 = &v0[offset*dim]; 1236 geom.J = &J[offset*dim*dim]; 1237 geom.invJ = &invJ[offset*dim*dim]; 1238 geom.detJ = &detJ[offset]; 1239 ierr = PetscFEIntegrateIFunction(fe[f], Nr, Nf, fe, f, geom, &u[offset*cellDof], &u_t[offset*cellDof], NfAux, feAux, &a[offset*cellDofAux], f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr); 1240 } 1241 for (c = cStart; c < cEnd; ++c) { 1242 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);} 1243 ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr); 1244 } 1245 ierr = PetscFree7(u,u_t,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr); 1246 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1247 if (feBd) { 1248 DMLabel label, depth; 1249 IS pointIS; 1250 const PetscInt *points; 1251 PetscInt dep, numPoints, p, numFaces; 1252 PetscReal *n; 1253 1254 ierr = DMPlexGetLabel(dm, "boundary", &label);CHKERRQ(ierr); 1255 ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr); 1256 ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr); 1257 ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr); 1258 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 1259 for (f = 0, cellDof = 0, numComponents = 0; f < Nf; ++f) { 1260 PetscInt Nb, Nc; 1261 1262 ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr); 1263 ierr = PetscFEGetNumComponents(feBd[f], &Nc);CHKERRQ(ierr); 1264 cellDof += Nb*Nc; 1265 numComponents += Nc; 1266 } 1267 for (p = 0, numFaces = 0; p < numPoints; ++p) { 1268 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 1269 if (dep == dim-1) ++numFaces; 1270 } 1271 ierr = PetscMalloc7(numFaces*cellDof,&u,numFaces*dim,&v0,numFaces*dim,&n,numFaces*dim*dim,&J,numFaces*dim*dim,&invJ,numFaces,&detJ,numFaces*cellDof,&elemVec);CHKERRQ(ierr); 1272 ierr = PetscMalloc1(numFaces*cellDof,&u_t);CHKERRQ(ierr); 1273 for (p = 0, f = 0; p < numPoints; ++p) { 1274 const PetscInt point = points[p]; 1275 PetscScalar *x = NULL; 1276 PetscInt i; 1277 1278 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 1279 if (dep != dim-1) continue; 1280 ierr = DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr); 1281 ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]); 1282 if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point); 1283 ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 1284 for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i]; 1285 ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 1286 ierr = DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr); 1287 for (i = 0; i < cellDof; ++i) u_t[f*cellDof+i] = x[i]; 1288 ierr = DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr); 1289 ++f; 1290 } 1291 for (f = 0; f < Nf; ++f) { 1292 void (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdIFuncs[f]; 1293 void (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdIFuncs[f]; 1294 PetscInt numQuadPoints, Nb; 1295 /* Conforming batches */ 1296 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 1297 /* Remainder */ 1298 PetscInt Nr, offset; 1299 1300 ierr = PetscFEGetQuadrature(feBd[f], &q);CHKERRQ(ierr); 1301 ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr); 1302 ierr = PetscFEGetTileSizes(feBd[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1303 ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 1304 blockSize = Nb*numQuadPoints; 1305 batchSize = numBlocks * blockSize; 1306 ierr = PetscFESetTileSizes(feBd[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1307 numChunks = numFaces / (numBatches*batchSize); 1308 Ne = numChunks*numBatches*batchSize; 1309 Nr = numFaces % (numBatches*batchSize); 1310 offset = numFaces - Nr; 1311 geom.v0 = v0; 1312 geom.n = n; 1313 geom.J = J; 1314 geom.invJ = invJ; 1315 geom.detJ = detJ; 1316 ierr = PetscFEIntegrateBdIFunction(feBd[f], Ne, Nf, feBd, f, geom, u, u_t, 0, NULL, NULL, f0, f1, elemVec);CHKERRQ(ierr); 1317 geom.v0 = &v0[offset*dim]; 1318 geom.n = &n[offset*dim]; 1319 geom.J = &J[offset*dim*dim]; 1320 geom.invJ = &invJ[offset*dim*dim]; 1321 geom.detJ = &detJ[offset]; 1322 ierr = PetscFEIntegrateBdIFunction(feBd[f], Nr, Nf, feBd, f, geom, &u[offset*cellDof], &u_t[offset*cellDof], 0, NULL, NULL, f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr); 1323 } 1324 for (p = 0, f = 0; p < numPoints; ++p) { 1325 const PetscInt point = points[p]; 1326 1327 ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr); 1328 if (dep != dim-1) continue; 1329 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "BdResidual", cellDof, &elemVec[f*cellDof]);CHKERRQ(ierr);} 1330 ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*cellDof], ADD_VALUES);CHKERRQ(ierr); 1331 ++f; 1332 } 1333 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 1334 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 1335 ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr); 1336 ierr = PetscFree(u_t);CHKERRQ(ierr); 1337 } 1338 if (mesh->printFEM) {ierr = DMPrintLocalVec(dm, name, mesh->printTol, F);CHKERRQ(ierr);} 1339 ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr); 1340 PetscFunctionReturn(0); 1341 } 1342 1343 #undef __FUNCT__ 1344 #define __FUNCT__ "DMPlexComputeJacobianActionFEM" 1345 /*@C 1346 DMPlexComputeJacobianActionFEM - Form the local action of Jacobian J(u) on the local input X using pointwise functions specified by the user 1347 1348 Input Parameters: 1349 + dm - The mesh 1350 . J - The Jacobian shell matrix 1351 . X - Local input vector 1352 - user - The user context 1353 1354 Output Parameter: 1355 . F - Local output vector 1356 1357 Note: 1358 The first member of the user context must be an FEMContext. 1359 1360 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 1361 like a GPU, or vectorize on a multicore machine. 1362 1363 Level: developer 1364 1365 .seealso: DMPlexComputeResidualFEM() 1366 @*/ 1367 PetscErrorCode DMPlexComputeJacobianActionFEM(DM dm, Mat Jac, Vec X, Vec F, void *user) 1368 { 1369 DM_Plex *mesh = (DM_Plex *) dm->data; 1370 PetscFEM *fem = (PetscFEM *) user; 1371 PetscFE *fe = fem->fe; 1372 PetscQuadrature quad; 1373 PetscCellGeometry geom; 1374 PetscSection section; 1375 JacActionCtx *jctx; 1376 PetscReal *v0, *J, *invJ, *detJ; 1377 PetscScalar *elemVec, *u, *a; 1378 PetscInt dim, numFields, field, numCells, cStart, cEnd, c; 1379 PetscInt cellDof = 0; 1380 PetscErrorCode ierr; 1381 1382 PetscFunctionBegin; 1383 /* ierr = PetscLogEventBegin(DMPLEX_JacobianActionFEM,dm,0,0,0);CHKERRQ(ierr); */ 1384 ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr); 1385 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1386 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1387 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1388 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1389 numCells = cEnd - cStart; 1390 for (field = 0; field < numFields; ++field) { 1391 PetscInt Nb, Nc; 1392 1393 ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr); 1394 ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr); 1395 cellDof += Nb*Nc; 1396 } 1397 ierr = VecSet(F, 0.0);CHKERRQ(ierr); 1398 ierr = PetscMalloc7(numCells*cellDof,&u,numCells*cellDof,&a,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof,&elemVec);CHKERRQ(ierr); 1399 for (c = cStart; c < cEnd; ++c) { 1400 PetscScalar *x = NULL; 1401 PetscInt i; 1402 1403 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 1404 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 1405 ierr = DMPlexVecGetClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr); 1406 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 1407 ierr = DMPlexVecRestoreClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr); 1408 ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 1409 for (i = 0; i < cellDof; ++i) a[c*cellDof+i] = x[i]; 1410 ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 1411 } 1412 for (field = 0; field < numFields; ++field) { 1413 PetscInt numQuadPoints, Nb; 1414 /* Conforming batches */ 1415 PetscInt numBlocks = 1; 1416 PetscInt numBatches = 1; 1417 PetscInt numChunks, Ne, blockSize, batchSize; 1418 /* Remainder */ 1419 PetscInt Nr, offset; 1420 1421 ierr = PetscFEGetQuadrature(fe[field], &quad);CHKERRQ(ierr); 1422 ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr); 1423 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 1424 blockSize = Nb*numQuadPoints; 1425 batchSize = numBlocks * blockSize; 1426 numChunks = numCells / (numBatches*batchSize); 1427 Ne = numChunks*numBatches*batchSize; 1428 Nr = numCells % (numBatches*batchSize); 1429 offset = numCells - Nr; 1430 geom.v0 = v0; 1431 geom.J = J; 1432 geom.invJ = invJ; 1433 geom.detJ = detJ; 1434 ierr = PetscFEIntegrateJacobianAction(fe[field], Ne, numFields, fe, field, geom, u, a, fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, elemVec);CHKERRQ(ierr); 1435 geom.v0 = &v0[offset*dim]; 1436 geom.J = &J[offset*dim*dim]; 1437 geom.invJ = &invJ[offset*dim*dim]; 1438 geom.detJ = &detJ[offset]; 1439 ierr = PetscFEIntegrateJacobianAction(fe[field], Nr, numFields, fe, field, geom, &u[offset*cellDof], &a[offset*cellDof], 1440 fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, &elemVec[offset*cellDof]);CHKERRQ(ierr); 1441 } 1442 for (c = cStart; c < cEnd; ++c) { 1443 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Jacobian Action", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);} 1444 ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr); 1445 } 1446 ierr = PetscFree7(u,a,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr); 1447 if (mesh->printFEM) { 1448 PetscMPIInt rank, numProcs; 1449 PetscInt p; 1450 1451 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 1452 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr); 1453 ierr = PetscPrintf(PetscObjectComm((PetscObject)dm), "Jacobian Action:\n");CHKERRQ(ierr); 1454 for (p = 0; p < numProcs; ++p) { 1455 if (p == rank) {ierr = VecView(F, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);} 1456 ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr); 1457 } 1458 } 1459 /* ierr = PetscLogEventEnd(DMPLEX_JacobianActionFEM,dm,0,0,0);CHKERRQ(ierr); */ 1460 PetscFunctionReturn(0); 1461 } 1462 1463 #undef __FUNCT__ 1464 #define __FUNCT__ "DMPlexComputeJacobianFEM" 1465 /*@ 1466 DMPlexComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user. 1467 1468 Input Parameters: 1469 + dm - The mesh 1470 . X - Local input vector 1471 - user - The user context 1472 1473 Output Parameter: 1474 . Jac - Jacobian matrix 1475 1476 Note: 1477 The first member of the user context must be an FEMContext. 1478 1479 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 1480 like a GPU, or vectorize on a multicore machine. 1481 1482 Level: developer 1483 1484 .seealso: FormFunctionLocal() 1485 @*/ 1486 PetscErrorCode DMPlexComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user) 1487 { 1488 DM_Plex *mesh = (DM_Plex *) dm->data; 1489 PetscFEM *fem = (PetscFEM *) user; 1490 PetscFE *fe = fem->fe; 1491 PetscFE *feAux = fem->feAux; 1492 PetscFE *feBd = fem->feBd; 1493 const char *name = "Jacobian"; 1494 DM dmAux; 1495 Vec A; 1496 PetscQuadrature quad; 1497 PetscCellGeometry geom; 1498 PetscSection section, globalSection, sectionAux; 1499 PetscReal *v0, *J, *invJ, *detJ; 1500 PetscScalar *elemMat, *u, *a; 1501 PetscInt dim, Nf, NfAux = 0, f, fieldI, fieldJ, numCells, cStart, cEnd, c; 1502 PetscInt cellDof = 0, numComponents = 0; 1503 PetscInt cellDofAux = 0, numComponentsAux = 0; 1504 PetscBool isShell; 1505 PetscErrorCode ierr; 1506 1507 PetscFunctionBegin; 1508 ierr = PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr); 1509 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1510 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1511 ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr); 1512 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1513 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1514 numCells = cEnd - cStart; 1515 for (f = 0; f < Nf; ++f) { 1516 PetscInt Nb, Nc; 1517 1518 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 1519 ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr); 1520 cellDof += Nb*Nc; 1521 numComponents += Nc; 1522 } 1523 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1524 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1525 if (dmAux) { 1526 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1527 ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr); 1528 } 1529 for (f = 0; f < NfAux; ++f) { 1530 PetscInt Nb, Nc; 1531 1532 ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr); 1533 ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr); 1534 cellDofAux += Nb*Nc; 1535 numComponentsAux += Nc; 1536 } 1537 ierr = DMPlexInsertBoundaryValuesFEM(dm, X);CHKERRQ(ierr); 1538 ierr = MatZeroEntries(JacP);CHKERRQ(ierr); 1539 ierr = PetscMalloc6(numCells*cellDof,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof*cellDof,&elemMat);CHKERRQ(ierr); 1540 if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);} 1541 for (c = cStart; c < cEnd; ++c) { 1542 PetscScalar *x = NULL; 1543 PetscInt i; 1544 1545 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 1546 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 1547 ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 1548 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 1549 ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 1550 if (dmAux) { 1551 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1552 for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x[i]; 1553 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1554 } 1555 } 1556 ierr = PetscMemzero(elemMat, numCells*cellDof*cellDof * sizeof(PetscScalar));CHKERRQ(ierr); 1557 for (fieldI = 0; fieldI < Nf; ++fieldI) { 1558 PetscInt numQuadPoints, Nb; 1559 /* Conforming batches */ 1560 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 1561 /* Remainder */ 1562 PetscInt Nr, offset; 1563 1564 ierr = PetscFEGetQuadrature(fe[fieldI], &quad);CHKERRQ(ierr); 1565 ierr = PetscFEGetDimension(fe[fieldI], &Nb);CHKERRQ(ierr); 1566 ierr = PetscFEGetTileSizes(fe[fieldI], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1567 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 1568 blockSize = Nb*numQuadPoints; 1569 batchSize = numBlocks * blockSize; 1570 ierr = PetscFESetTileSizes(fe[fieldI], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1571 numChunks = numCells / (numBatches*batchSize); 1572 Ne = numChunks*numBatches*batchSize; 1573 Nr = numCells % (numBatches*batchSize); 1574 offset = numCells - Nr; 1575 for (fieldJ = 0; fieldJ < Nf; ++fieldJ) { 1576 void (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g0Funcs[fieldI*Nf+fieldJ]; 1577 void (*g1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g1Funcs[fieldI*Nf+fieldJ]; 1578 void (*g2)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g2Funcs[fieldI*Nf+fieldJ]; 1579 void (*g3)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g3Funcs[fieldI*Nf+fieldJ]; 1580 1581 geom.v0 = v0; 1582 geom.J = J; 1583 geom.invJ = invJ; 1584 geom.detJ = detJ; 1585 ierr = PetscFEIntegrateJacobian(fe[fieldI], Ne, Nf, fe, fieldI, fieldJ, geom, u, NfAux, feAux, a, g0, g1, g2, g3, elemMat);CHKERRQ(ierr); 1586 geom.v0 = &v0[offset*dim]; 1587 geom.J = &J[offset*dim*dim]; 1588 geom.invJ = &invJ[offset*dim*dim]; 1589 geom.detJ = &detJ[offset]; 1590 ierr = PetscFEIntegrateJacobian(fe[fieldI], Nr, Nf, fe, fieldI, fieldJ, geom, &u[offset*cellDof], NfAux, feAux, &a[offset*cellDofAux], g0, g1, g2, g3, &elemMat[offset*cellDof*cellDof]);CHKERRQ(ierr); 1591 } 1592 } 1593 for (c = cStart; c < cEnd; ++c) { 1594 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, name, cellDof, cellDof, &elemMat[c*cellDof*cellDof]);CHKERRQ(ierr);} 1595 ierr = DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMat[c*cellDof*cellDof], ADD_VALUES);CHKERRQ(ierr); 1596 } 1597 ierr = PetscFree6(u,v0,J,invJ,detJ,elemMat);CHKERRQ(ierr); 1598 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1599 if (feBd) { 1600 DMLabel depth; 1601 PetscInt numBd, bd; 1602 1603 for (f = 0, cellDof = 0, numComponents = 0; f < Nf; ++f) { 1604 PetscInt Nb, Nc; 1605 1606 ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr); 1607 ierr = PetscFEGetNumComponents(feBd[f], &Nc);CHKERRQ(ierr); 1608 cellDof += Nb*Nc; 1609 numComponents += Nc; 1610 } 1611 ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr); 1612 ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 1613 for (bd = 0; bd < numBd; ++bd) { 1614 const char *bdLabel; 1615 DMLabel label; 1616 IS pointIS; 1617 const PetscInt *points; 1618 const PetscInt *values; 1619 PetscReal *n; 1620 PetscInt field, numValues, numPoints, p, dep, numFaces; 1621 PetscBool isEssential; 1622 1623 ierr = DMPlexGetBoundary(dm, bd, &isEssential, &bdLabel, &field, NULL, &numValues, &values, NULL);CHKERRQ(ierr); 1624 if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this"); 1625 ierr = DMPlexGetLabel(dm, bdLabel, &label);CHKERRQ(ierr); 1626 ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr); 1627 ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr); 1628 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 1629 for (p = 0, numFaces = 0; p < numPoints; ++p) { 1630 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 1631 if (dep == dim-1) ++numFaces; 1632 } 1633 ierr = PetscMalloc7(numFaces*cellDof,&u,numFaces*dim,&v0,numFaces*dim,&n,numFaces*dim*dim,&J,numFaces*dim*dim,&invJ,numFaces,&detJ,numFaces*cellDof*cellDof,&elemMat);CHKERRQ(ierr); 1634 for (p = 0, f = 0; p < numPoints; ++p) { 1635 const PetscInt point = points[p]; 1636 PetscScalar *x = NULL; 1637 PetscInt i; 1638 1639 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 1640 if (dep != dim-1) continue; 1641 ierr = DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr); 1642 ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]); 1643 if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point); 1644 ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 1645 for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i]; 1646 ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 1647 ++f; 1648 } 1649 ierr = PetscMemzero(elemMat, numFaces*cellDof*cellDof * sizeof(PetscScalar));CHKERRQ(ierr); 1650 for (fieldI = 0; fieldI < Nf; ++fieldI) { 1651 PetscInt numQuadPoints, Nb; 1652 /* Conforming batches */ 1653 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 1654 /* Remainder */ 1655 PetscInt Nr, offset; 1656 1657 ierr = PetscFEGetQuadrature(feBd[fieldI], &quad);CHKERRQ(ierr); 1658 ierr = PetscFEGetDimension(feBd[fieldI], &Nb);CHKERRQ(ierr); 1659 ierr = PetscFEGetTileSizes(feBd[fieldI], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1660 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 1661 blockSize = Nb*numQuadPoints; 1662 batchSize = numBlocks * blockSize; 1663 ierr = PetscFESetTileSizes(feBd[fieldI], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1664 numChunks = numFaces / (numBatches*batchSize); 1665 Ne = numChunks*numBatches*batchSize; 1666 Nr = numFaces % (numBatches*batchSize); 1667 offset = numFaces - Nr; 1668 for (fieldJ = 0; fieldJ < Nf; ++fieldJ) { 1669 void (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->g0BdFuncs[fieldI*Nf+fieldJ]; 1670 void (*g1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->g1BdFuncs[fieldI*Nf+fieldJ]; 1671 void (*g2)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->g2BdFuncs[fieldI*Nf+fieldJ]; 1672 void (*g3)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->g3BdFuncs[fieldI*Nf+fieldJ]; 1673 1674 geom.v0 = v0; 1675 geom.n = n; 1676 geom.J = J; 1677 geom.invJ = invJ; 1678 geom.detJ = detJ; 1679 ierr = PetscFEIntegrateBdJacobian(feBd[fieldI], Ne, Nf, feBd, fieldI, fieldJ, geom, u, 0, NULL, NULL, g0, g1, g2, g3, elemMat);CHKERRQ(ierr); 1680 geom.v0 = &v0[offset*dim]; 1681 geom.n = &n[offset*dim]; 1682 geom.J = &J[offset*dim*dim]; 1683 geom.invJ = &invJ[offset*dim*dim]; 1684 geom.detJ = &detJ[offset]; 1685 ierr = PetscFEIntegrateBdJacobian(feBd[fieldI], Nr, Nf, feBd, fieldI, fieldJ, geom, &u[offset*cellDof], 0, NULL, NULL, g0, g1, g2, g3, &elemMat[offset*cellDof*cellDof]);CHKERRQ(ierr); 1686 } 1687 } 1688 for (p = 0, f = 0; p < numPoints; ++p) { 1689 const PetscInt point = points[p]; 1690 1691 ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr); 1692 if (dep != dim-1) continue; 1693 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(point, "BdJacobian", cellDof, cellDof, &elemMat[f*cellDof*cellDof]);CHKERRQ(ierr);} 1694 ierr = DMPlexMatSetClosure(dm, section, globalSection, JacP, point, &elemMat[f*cellDof*cellDof], ADD_VALUES);CHKERRQ(ierr); 1695 ++f; 1696 } 1697 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 1698 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 1699 ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemMat);CHKERRQ(ierr); 1700 } 1701 } 1702 ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1703 ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1704 if (mesh->printFEM) { 1705 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1706 ierr = MatChop(JacP, 1.0e-10);CHKERRQ(ierr); 1707 ierr = MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1708 } 1709 ierr = PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr); 1710 ierr = PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);CHKERRQ(ierr); 1711 if (isShell) { 1712 JacActionCtx *jctx; 1713 1714 ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr); 1715 ierr = VecCopy(X, jctx->u);CHKERRQ(ierr); 1716 } 1717 PetscFunctionReturn(0); 1718 } 1719 1720 #undef __FUNCT__ 1721 #define __FUNCT__ "DMPlexComputeInterpolatorFEM" 1722 /*@ 1723 DMPlexComputeInterpolatorFEM - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1724 1725 Input Parameters: 1726 + dmf - The fine mesh 1727 . dmc - The coarse mesh 1728 - user - The user context 1729 1730 Output Parameter: 1731 . In - The interpolation matrix 1732 1733 Note: 1734 The first member of the user context must be an FEMContext. 1735 1736 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 1737 like a GPU, or vectorize on a multicore machine. 1738 1739 Level: developer 1740 1741 .seealso: DMPlexComputeJacobianFEM() 1742 @*/ 1743 PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user) 1744 { 1745 DM_Plex *mesh = (DM_Plex *) dmc->data; 1746 PetscFEM *fem = (PetscFEM *) user; 1747 PetscFE *fe = fem->fe; 1748 const char *name = "Interpolator"; 1749 PetscFE *feRef; 1750 PetscSection fsection, fglobalSection; 1751 PetscSection csection, cglobalSection; 1752 PetscScalar *elemMat; 1753 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, c; 1754 PetscInt rCellDof = 0, cCellDof = 0; 1755 PetscErrorCode ierr; 1756 1757 PetscFunctionBegin; 1758 #if 0 1759 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1760 #endif 1761 ierr = DMPlexGetDimension(dmf, &dim);CHKERRQ(ierr); 1762 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1763 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1764 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1765 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1766 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1767 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1768 ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr); 1769 for (f = 0; f < Nf; ++f) { 1770 PetscInt rNb, cNb, Nc; 1771 1772 ierr = PetscFERefine(fe[f], &feRef[f]);CHKERRQ(ierr); 1773 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1774 ierr = PetscFEGetDimension(fe[f], &cNb);CHKERRQ(ierr); 1775 ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr); 1776 rCellDof += rNb*Nc; 1777 cCellDof += cNb*Nc; 1778 } 1779 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1780 ierr = PetscMalloc1(rCellDof*cCellDof,&elemMat);CHKERRQ(ierr); 1781 ierr = PetscMemzero(elemMat, rCellDof*cCellDof * sizeof(PetscScalar));CHKERRQ(ierr); 1782 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1783 PetscDualSpace Qref; 1784 PetscQuadrature f; 1785 const PetscReal *qpoints, *qweights; 1786 PetscReal *points; 1787 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1788 1789 /* Compose points from all dual basis functionals */ 1790 ierr = PetscFEGetNumComponents(fe[fieldI], &Nc);CHKERRQ(ierr); 1791 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1792 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1793 for (i = 0; i < fpdim; ++i) { 1794 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1795 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1796 npoints += Np; 1797 } 1798 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1799 for (i = 0, k = 0; i < fpdim; ++i) { 1800 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1801 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1802 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1803 } 1804 1805 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1806 PetscReal *B; 1807 PetscInt NcJ, cpdim, j; 1808 1809 /* Evaluate basis at points */ 1810 ierr = PetscFEGetNumComponents(fe[fieldJ], &NcJ);CHKERRQ(ierr); 1811 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); 1812 ierr = PetscFEGetDimension(fe[fieldJ], &cpdim);CHKERRQ(ierr); 1813 /* For now, fields only interpolate themselves */ 1814 if (fieldI == fieldJ) { 1815 ierr = PetscFEGetTabulation(fe[fieldJ], npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1816 for (i = 0, k = 0; i < fpdim; ++i) { 1817 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1818 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1819 for (p = 0; p < Np; ++p, ++k) { 1820 for (j = 0; j < cpdim; ++j) { 1821 for (c = 0; c < Nc; ++c) elemMat[(offsetI + i*Nc + c)*cCellDof + offsetJ + j*NcJ + c] += B[k*cpdim*NcJ+j*Nc+c]*qweights[p]; 1822 } 1823 } 1824 } 1825 ierr = PetscFERestoreTabulation(fe[fieldJ], npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1826 } 1827 offsetJ += cpdim*NcJ; 1828 } 1829 offsetI += fpdim*Nc; 1830 ierr = PetscFree(points);CHKERRQ(ierr); 1831 } 1832 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rCellDof, cCellDof, elemMat);CHKERRQ(ierr);} 1833 for (c = cStart; c < cEnd; ++c) { 1834 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1835 } 1836 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1837 ierr = PetscFree(feRef);CHKERRQ(ierr); 1838 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1839 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1840 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1841 if (mesh->printFEM) { 1842 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1843 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1844 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1845 } 1846 #if 0 1847 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1848 #endif 1849 PetscFunctionReturn(0); 1850 } 1851 1852 #undef __FUNCT__ 1853 #define __FUNCT__ "DMPlexAddBoundary" 1854 /* The ids can be overridden by the command line option -bc_<boundary name> */ 1855 PetscErrorCode DMPlexAddBoundary(DM dm, PetscBool isEssential, const char name[], PetscInt field, void (*bcFunc)(), PetscInt numids, const PetscInt *ids, void *ctx) 1856 { 1857 DM_Plex *mesh = (DM_Plex *) dm->data; 1858 DMBoundary b; 1859 PetscErrorCode ierr; 1860 1861 PetscFunctionBegin; 1862 ierr = PetscNew(&b);CHKERRQ(ierr); 1863 ierr = PetscStrallocpy(name, (char **) &b->name);CHKERRQ(ierr); 1864 ierr = PetscMalloc1(numids, &b->ids);CHKERRQ(ierr); 1865 ierr = PetscMemcpy(b->ids, ids, numids*sizeof(PetscInt));CHKERRQ(ierr); 1866 b->essential = isEssential; 1867 b->field = field; 1868 b->func = bcFunc; 1869 b->numids = numids; 1870 b->ctx = ctx; 1871 b->next = mesh->boundary; 1872 mesh->boundary = b; 1873 PetscFunctionReturn(0); 1874 } 1875 1876 #undef __FUNCT__ 1877 #define __FUNCT__ "DMPlexGetNumBoundary" 1878 PetscErrorCode DMPlexGetNumBoundary(DM dm, PetscInt *numBd) 1879 { 1880 DM_Plex *mesh = (DM_Plex *) dm->data; 1881 DMBoundary b = mesh->boundary; 1882 1883 PetscFunctionBegin; 1884 *numBd = 0; 1885 while (b) {++(*numBd); b = b->next;} 1886 PetscFunctionReturn(0); 1887 } 1888 1889 #undef __FUNCT__ 1890 #define __FUNCT__ "DMPlexGetBoundary" 1891 PetscErrorCode DMPlexGetBoundary(DM dm, PetscInt bd, PetscBool *isEssential, const char **name, PetscInt *field, void (**func)(), PetscInt *numids, const PetscInt **ids, void **ctx) 1892 { 1893 DM_Plex *mesh = (DM_Plex *) dm->data; 1894 DMBoundary b = mesh->boundary; 1895 PetscInt n = 0; 1896 1897 PetscFunctionBegin; 1898 while (b) { 1899 if (n == bd) break; 1900 b = b->next; 1901 ++n; 1902 } 1903 if (n != bd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n); 1904 if (isEssential) { 1905 PetscValidPointer(isEssential, 3); 1906 *isEssential = b->essential; 1907 } 1908 if (name) { 1909 PetscValidPointer(name, 4); 1910 *name = b->name; 1911 } 1912 if (field) { 1913 PetscValidPointer(field, 5); 1914 *field = b->field; 1915 } 1916 if (func) { 1917 PetscValidPointer(func, 6); 1918 *func = b->func; 1919 } 1920 if (numids) { 1921 PetscValidPointer(numids, 7); 1922 *numids = b->numids; 1923 } 1924 if (ids) { 1925 PetscValidPointer(ids, 8); 1926 *ids = b->ids; 1927 } 1928 if (ctx) { 1929 PetscValidPointer(ctx, 9); 1930 *ctx = b->ctx; 1931 } 1932 PetscFunctionReturn(0); 1933 } 1934