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 PetscInt numFields, numComp, dim, spDim, totDim = 0, numValues, cStart, cEnd, f, d, v, i, comp; 192 PetscErrorCode ierr; 193 194 PetscFunctionBegin; 195 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 196 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 197 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 198 ierr = PetscMalloc3(numFields,&sp,dim,&v0,dim*dim,&J);CHKERRQ(ierr); 199 for (f = 0; f < numFields; ++f) { 200 ierr = PetscFEGetDualSpace(fe[f], &sp[f]);CHKERRQ(ierr); 201 ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr); 202 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 203 totDim += spDim*numComp; 204 } 205 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 206 ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 207 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 208 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 209 for (i = 0; i < numIds; ++i) { 210 IS pointIS; 211 const PetscInt *points; 212 PetscInt n, p; 213 214 ierr = DMLabelGetStratumIS(label, ids[i], &pointIS);CHKERRQ(ierr); 215 ierr = ISGetLocalSize(pointIS, &n);CHKERRQ(ierr); 216 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 217 for (p = 0; p < n; ++p) { 218 const PetscInt point = points[p]; 219 PetscCellGeometry geom; 220 221 if ((point < cStart) || (point >= cEnd)) continue; 222 ierr = DMPlexComputeCellGeometry(dm, point, v0, J, NULL, &detJ);CHKERRQ(ierr); 223 geom.v0 = v0; 224 geom.J = J; 225 geom.detJ = &detJ; 226 for (f = 0, v = 0; f < numFields; ++f) { 227 void * const ctx = ctxs ? ctxs[f] : NULL; 228 ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr); 229 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 230 for (d = 0; d < spDim; ++d) { 231 if (funcs[f]) { 232 ierr = PetscDualSpaceApply(sp[f], d, geom, numComp, funcs[f], ctx, &values[v]);CHKERRQ(ierr); 233 } else { 234 for (comp = 0; comp < numComp; ++comp) values[v+comp] = 0.0; 235 } 236 v += numComp; 237 } 238 } 239 ierr = DMPlexVecSetClosure(dm, section, localX, point, values, mode);CHKERRQ(ierr); 240 } 241 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 242 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 243 } 244 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 245 ierr = PetscFree3(sp,v0,J);CHKERRQ(ierr); 246 PetscFunctionReturn(0); 247 } 248 249 #undef __FUNCT__ 250 #define __FUNCT__ "DMPlexProjectFunctionLocal" 251 PetscErrorCode DMPlexProjectFunctionLocal(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 252 { 253 PetscDualSpace *sp; 254 PetscSection section; 255 PetscScalar *values; 256 PetscReal *v0, *J, detJ; 257 PetscInt numFields, numComp, dim, spDim, totDim = 0, numValues, cStart, cEnd, c, f, d, v, comp; 258 PetscErrorCode ierr; 259 260 PetscFunctionBegin; 261 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 262 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 263 ierr = PetscMalloc1(numFields, &sp);CHKERRQ(ierr); 264 for (f = 0; f < numFields; ++f) { 265 ierr = PetscFEGetDualSpace(fe[f], &sp[f]);CHKERRQ(ierr); 266 ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr); 267 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 268 totDim += spDim*numComp; 269 } 270 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 271 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 272 ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 273 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 274 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 275 ierr = PetscMalloc2(dim,&v0,dim*dim,&J);CHKERRQ(ierr); 276 for (c = cStart; c < cEnd; ++c) { 277 PetscCellGeometry geom; 278 279 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);CHKERRQ(ierr); 280 geom.v0 = v0; 281 geom.J = J; 282 geom.detJ = &detJ; 283 for (f = 0, v = 0; f < numFields; ++f) { 284 void * const ctx = ctxs ? ctxs[f] : NULL; 285 ierr = PetscFEGetNumComponents(fe[f], &numComp);CHKERRQ(ierr); 286 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 287 for (d = 0; d < spDim; ++d) { 288 if (funcs[f]) { 289 ierr = PetscDualSpaceApply(sp[f], d, geom, numComp, funcs[f], ctx, &values[v]);CHKERRQ(ierr); 290 } else { 291 for (comp = 0; comp < numComp; ++comp) values[v+comp] = 0.0; 292 } 293 v += numComp; 294 } 295 } 296 ierr = DMPlexVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr); 297 } 298 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 299 ierr = PetscFree2(v0,J);CHKERRQ(ierr); 300 ierr = PetscFree(sp);CHKERRQ(ierr); 301 PetscFunctionReturn(0); 302 } 303 304 #undef __FUNCT__ 305 #define __FUNCT__ "DMPlexProjectFunction" 306 /*@C 307 DMPlexProjectFunction - This projects the given function into the function space provided. 308 309 Input Parameters: 310 + dm - The DM 311 . fe - The PetscFE associated with the field 312 . funcs - The coordinate functions to evaluate, one per field 313 . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null. 314 - mode - The insertion mode for values 315 316 Output Parameter: 317 . X - vector 318 319 Level: developer 320 321 .seealso: DMPlexComputeL2Diff() 322 @*/ 323 PetscErrorCode DMPlexProjectFunction(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X) 324 { 325 Vec localX; 326 PetscErrorCode ierr; 327 328 PetscFunctionBegin; 329 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 330 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 331 ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, ctxs, mode, localX);CHKERRQ(ierr); 332 ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr); 333 ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr); 334 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 335 PetscFunctionReturn(0); 336 } 337 338 #undef __FUNCT__ 339 #define __FUNCT__ "DMPlexInsertBoundaryValuesFEM" 340 PetscErrorCode DMPlexInsertBoundaryValuesFEM(DM dm, Vec localX) 341 { 342 void (**funcs)(const PetscReal x[], PetscScalar *u, void *ctx); 343 void **ctxs; 344 PetscFE *fe; 345 PetscInt numFields, f, numBd, b; 346 PetscErrorCode ierr; 347 348 PetscFunctionBegin; 349 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 350 PetscValidHeaderSpecific(localX, VEC_CLASSID, 2); 351 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 352 ierr = PetscMalloc3(numFields,&fe,numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 353 for (f = 0; f < numFields; ++f) {ierr = DMGetField(dm, f, (PetscObject *) &fe[f]);CHKERRQ(ierr);} 354 /* OPT: Could attempt to do multiple BCs at once */ 355 ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 356 for (b = 0; b < numBd; ++b) { 357 DMLabel label; 358 const PetscInt *ids; 359 const char *name; 360 PetscInt numids, field; 361 PetscBool isEssential; 362 void (*func)(); 363 void *ctx; 364 365 /* TODO: We need to set only the part indicated by the ids */ 366 ierr = DMPlexGetBoundary(dm, b, &isEssential, &name, &field, &func, &numids, &ids, &ctx);CHKERRQ(ierr); 367 ierr = DMPlexGetLabel(dm, name, &label);CHKERRQ(ierr); 368 for (f = 0; f < numFields; ++f) { 369 funcs[f] = field == f ? (void (*)(const PetscReal[], PetscScalar *, void *)) func : NULL; 370 ctxs[f] = field == f ? ctx : NULL; 371 } 372 ierr = DMPlexProjectFunctionLabelLocal(dm, label, numids, ids, fe, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 373 } 374 ierr = PetscFree3(fe,funcs,ctxs);CHKERRQ(ierr); 375 PetscFunctionReturn(0); 376 } 377 378 /* Assuming dim == 3 */ 379 typedef struct { 380 PetscScalar normal[3]; /* Area-scaled normals */ 381 PetscScalar centroid[3]; /* Location of centroid (quadrature point) */ 382 PetscScalar grad[2][3]; /* Face contribution to gradient in left and right cell */ 383 } FaceGeom; 384 385 #undef __FUNCT__ 386 #define __FUNCT__ "DMPlexInsertBoundaryValuesFVM" 387 PetscErrorCode DMPlexInsertBoundaryValuesFVM(DM dm, PetscReal time, Vec locX) 388 { 389 DM dmFace; 390 Vec faceGeometry; 391 DMLabel label; 392 const PetscScalar *facegeom; 393 PetscScalar *x; 394 PetscInt numBd, b, fStart, fEnd; 395 PetscErrorCode ierr; 396 397 PetscFunctionBegin; 398 /* TODO Pull this geometry calculation into the library */ 399 ierr = PetscObjectQuery((PetscObject) dm, "FaceGeometry", (PetscObject *) &faceGeometry);CHKERRQ(ierr); 400 ierr = DMPlexGetLabel(dm, "Face Sets", &label);CHKERRQ(ierr); 401 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 402 ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 403 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 404 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 405 ierr = VecGetArray(locX, &x);CHKERRQ(ierr); 406 for (b = 0; b < numBd; ++b) { 407 PetscErrorCode (*func)(PetscReal,const PetscScalar*,const PetscScalar*,const PetscScalar*,PetscScalar*,void*); 408 const PetscInt *ids; 409 PetscInt numids, i; 410 void *ctx; 411 412 ierr = DMPlexGetBoundary(dm, b, NULL, NULL, NULL, (void (**)()) &func, &numids, &ids, &ctx);CHKERRQ(ierr); 413 for (i = 0; i < numids; ++i) { 414 IS faceIS; 415 const PetscInt *faces; 416 PetscInt numFaces, f; 417 418 ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr); 419 if (!faceIS) continue; /* No points with that id on this process */ 420 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 421 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 422 for (f = 0; f < numFaces; ++f) { 423 const PetscInt face = faces[f], *cells; 424 const PetscScalar *xI; 425 PetscScalar *xG; 426 const FaceGeom *fg; 427 428 if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */ 429 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 430 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 431 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr); 432 ierr = DMPlexPointLocalRef(dm, cells[1], x, &xG);CHKERRQ(ierr); 433 ierr = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);CHKERRQ(ierr); 434 } 435 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 436 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 437 } 438 } 439 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 440 ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr); 441 PetscFunctionReturn(0); 442 } 443 444 #undef __FUNCT__ 445 #define __FUNCT__ "DMPlexComputeL2Diff" 446 /*@C 447 DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h. 448 449 Input Parameters: 450 + dm - The DM 451 . fe - The PetscFE object for each field 452 . funcs - The functions to evaluate for each field component 453 . ctxs - Optional array of contexts to pass to each function, or NULL. 454 - X - The coefficient vector u_h 455 456 Output Parameter: 457 . diff - The diff ||u - u_h||_2 458 459 Level: developer 460 461 .seealso: DMPlexProjectFunction(), DMPlexComputeL2GradientDiff() 462 @*/ 463 PetscErrorCode DMPlexComputeL2Diff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 464 { 465 const PetscInt debug = 0; 466 PetscSection section; 467 PetscQuadrature quad; 468 Vec localX; 469 PetscScalar *funcVal; 470 PetscReal *coords, *v0, *J, *invJ, detJ; 471 PetscReal localDiff = 0.0; 472 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp; 473 PetscErrorCode ierr; 474 475 PetscFunctionBegin; 476 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 477 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 478 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 479 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 480 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 481 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 482 for (field = 0; field < numFields; ++field) { 483 PetscInt Nc; 484 485 ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr); 486 numComponents += Nc; 487 } 488 ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 489 ierr = PetscMalloc5(numComponents,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 490 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 491 ierr = PetscFEGetQuadrature(fe[0], &quad);CHKERRQ(ierr); 492 for (c = cStart; c < cEnd; ++c) { 493 PetscScalar *x = NULL; 494 PetscReal elemDiff = 0.0; 495 496 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr); 497 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 498 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 499 500 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 501 void * const ctx = ctxs ? ctxs[field] : NULL; 502 const PetscReal *quadPoints, *quadWeights; 503 PetscReal *basis; 504 PetscInt numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f; 505 506 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 507 ierr = PetscFEGetDimension(fe[field], &numBasisFuncs);CHKERRQ(ierr); 508 ierr = PetscFEGetNumComponents(fe[field], &numBasisComps);CHKERRQ(ierr); 509 ierr = PetscFEGetDefaultTabulation(fe[field], &basis, NULL, NULL);CHKERRQ(ierr); 510 if (debug) { 511 char title[1024]; 512 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 513 ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr); 514 } 515 for (q = 0; q < numQuadPoints; ++q) { 516 for (d = 0; d < dim; d++) { 517 coords[d] = v0[d]; 518 for (e = 0; e < dim; e++) { 519 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 520 } 521 } 522 (*funcs[field])(coords, funcVal, ctx); 523 for (fc = 0; fc < numBasisComps; ++fc) { 524 PetscScalar interpolant = 0.0; 525 526 for (f = 0; f < numBasisFuncs; ++f) { 527 const PetscInt fidx = f*numBasisComps+fc; 528 interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx]; 529 } 530 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);} 531 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 532 } 533 } 534 comp += numBasisComps; 535 fieldOffset += numBasisFuncs*numBasisComps; 536 } 537 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 538 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 539 localDiff += elemDiff; 540 } 541 ierr = PetscFree5(funcVal,coords,v0,J,invJ);CHKERRQ(ierr); 542 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 543 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 544 *diff = PetscSqrtReal(*diff); 545 PetscFunctionReturn(0); 546 } 547 548 #undef __FUNCT__ 549 #define __FUNCT__ "DMPlexComputeL2GradientDiff" 550 /*@C 551 DMPlexComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h. 552 553 Input Parameters: 554 + dm - The DM 555 . fe - The PetscFE object for each field 556 . funcs - The gradient functions to evaluate for each field component 557 . ctxs - Optional array of contexts to pass to each function, or NULL. 558 . X - The coefficient vector u_h 559 - n - The vector to project along 560 561 Output Parameter: 562 . diff - The diff ||(grad u - grad u_h) . n||_2 563 564 Level: developer 565 566 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff() 567 @*/ 568 PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff) 569 { 570 const PetscInt debug = 0; 571 PetscSection section; 572 PetscQuadrature quad; 573 Vec localX; 574 PetscScalar *funcVal, *interpolantVec; 575 PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ; 576 PetscReal localDiff = 0.0; 577 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp; 578 PetscErrorCode ierr; 579 580 PetscFunctionBegin; 581 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 582 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 583 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 584 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 585 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 586 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 587 for (field = 0; field < numFields; ++field) { 588 PetscInt Nc; 589 590 ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr); 591 numComponents += Nc; 592 } 593 /* ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 594 ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr); 595 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 596 ierr = PetscFEGetQuadrature(fe[0], &quad);CHKERRQ(ierr); 597 for (c = cStart; c < cEnd; ++c) { 598 PetscScalar *x = NULL; 599 PetscReal elemDiff = 0.0; 600 601 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr); 602 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 603 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 604 605 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 606 void * const ctx = ctxs ? ctxs[field] : NULL; 607 const PetscReal *quadPoints, *quadWeights; 608 PetscReal *basisDer; 609 PetscInt numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g; 610 611 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 612 ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr); 613 ierr = PetscFEGetNumComponents(fe[field], &Ncomp);CHKERRQ(ierr); 614 ierr = PetscFEGetDefaultTabulation(fe[field], NULL, &basisDer, NULL);CHKERRQ(ierr); 615 if (debug) { 616 char title[1024]; 617 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 618 ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr); 619 } 620 for (q = 0; q < numQuadPoints; ++q) { 621 for (d = 0; d < dim; d++) { 622 coords[d] = v0[d]; 623 for (e = 0; e < dim; e++) { 624 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 625 } 626 } 627 (*funcs[field])(coords, n, funcVal, ctx); 628 for (fc = 0; fc < Ncomp; ++fc) { 629 PetscScalar interpolant = 0.0; 630 631 for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0; 632 for (f = 0; f < Nb; ++f) { 633 const PetscInt fidx = f*Ncomp+fc; 634 635 for (d = 0; d < dim; ++d) { 636 realSpaceDer[d] = 0.0; 637 for (g = 0; g < dim; ++g) { 638 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g]; 639 } 640 interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d]; 641 } 642 } 643 for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d]; 644 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);} 645 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 646 } 647 } 648 comp += Ncomp; 649 fieldOffset += Nb*Ncomp; 650 } 651 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 652 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 653 localDiff += elemDiff; 654 } 655 ierr = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr); 656 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 657 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 658 *diff = PetscSqrtReal(*diff); 659 PetscFunctionReturn(0); 660 } 661 662 #undef __FUNCT__ 663 #define __FUNCT__ "DMPlexComputeResidualFEM" 664 /*@ 665 DMPlexComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user 666 667 Input Parameters: 668 + dm - The mesh 669 . X - Local input vector 670 - user - The user context 671 672 Output Parameter: 673 . F - Local output vector 674 675 Note: 676 The first member of the user context must be an FEMContext. 677 678 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 679 like a GPU, or vectorize on a multicore machine. 680 681 Level: developer 682 683 .seealso: DMPlexComputeJacobianActionFEM() 684 @*/ 685 PetscErrorCode DMPlexComputeResidualFEM(DM dm, Vec X, Vec F, void *user) 686 { 687 DM_Plex *mesh = (DM_Plex *) dm->data; 688 PetscFEM *fem = (PetscFEM *) user; 689 PetscFE *fe = fem->fe; 690 PetscFE *feAux = fem->feAux; 691 PetscFE *feBd = fem->feBd; 692 const char *name = "Residual"; 693 DM dmAux; 694 Vec A; 695 PetscQuadrature q; 696 PetscCellGeometry geom; 697 PetscSection section, sectionAux; 698 PetscReal *v0, *J, *invJ, *detJ; 699 PetscScalar *elemVec, *u, *a = NULL; 700 PetscInt dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, c; 701 PetscInt cellDof = 0, numComponents = 0; 702 PetscInt cellDofAux = 0, numComponentsAux = 0; 703 PetscErrorCode ierr; 704 705 PetscFunctionBegin; 706 ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr); 707 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 708 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 709 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 710 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 711 numCells = cEnd - cStart; 712 for (f = 0; f < Nf; ++f) { 713 PetscInt Nb, Nc; 714 715 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 716 ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr); 717 cellDof += Nb*Nc; 718 numComponents += Nc; 719 } 720 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 721 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 722 if (dmAux) { 723 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 724 ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr); 725 } 726 for (f = 0; f < NfAux; ++f) { 727 PetscInt Nb, Nc; 728 729 ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr); 730 ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr); 731 cellDofAux += Nb*Nc; 732 numComponentsAux += Nc; 733 } 734 ierr = DMPlexInsertBoundaryValuesFEM(dm, X);CHKERRQ(ierr); 735 ierr = VecSet(F, 0.0);CHKERRQ(ierr); 736 ierr = PetscMalloc6(numCells*cellDof,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof,&elemVec);CHKERRQ(ierr); 737 if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);} 738 for (c = cStart; c < cEnd; ++c) { 739 PetscScalar *x = NULL; 740 PetscInt i; 741 742 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 743 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 744 ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 745 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 746 ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 747 if (dmAux) { 748 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 749 for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x[i]; 750 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 751 } 752 } 753 for (f = 0; f < Nf; ++f) { 754 void (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0Funcs[f]; 755 void (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1Funcs[f]; 756 PetscInt numQuadPoints, Nb; 757 /* Conforming batches */ 758 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 759 /* Remainder */ 760 PetscInt Nr, offset; 761 762 ierr = PetscFEGetQuadrature(fe[f], &q);CHKERRQ(ierr); 763 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 764 ierr = PetscFEGetTileSizes(fe[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 765 ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 766 blockSize = Nb*numQuadPoints; 767 batchSize = numBlocks * blockSize; 768 ierr = PetscFESetTileSizes(fe[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 769 numChunks = numCells / (numBatches*batchSize); 770 Ne = numChunks*numBatches*batchSize; 771 Nr = numCells % (numBatches*batchSize); 772 offset = numCells - Nr; 773 geom.v0 = v0; 774 geom.J = J; 775 geom.invJ = invJ; 776 geom.detJ = detJ; 777 ierr = PetscFEIntegrateResidual(fe[f], Ne, Nf, fe, f, geom, u, NfAux, feAux, a, f0, f1, elemVec);CHKERRQ(ierr); 778 geom.v0 = &v0[offset*dim]; 779 geom.J = &J[offset*dim*dim]; 780 geom.invJ = &invJ[offset*dim*dim]; 781 geom.detJ = &detJ[offset]; 782 ierr = PetscFEIntegrateResidual(fe[f], Nr, Nf, fe, f, geom, &u[offset*cellDof], NfAux, feAux, &a[offset*cellDofAux], f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr); 783 } 784 for (c = cStart; c < cEnd; ++c) { 785 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);} 786 ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr); 787 } 788 ierr = PetscFree6(u,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr); 789 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 790 if (feBd) { 791 DMLabel label, depth; 792 IS pointIS; 793 const PetscInt *points; 794 PetscInt dep, numPoints, p, numFaces; 795 PetscReal *n; 796 797 ierr = DMPlexGetLabel(dm, "boundary", &label);CHKERRQ(ierr); 798 ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr); 799 ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr); 800 ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr); 801 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 802 for (f = 0, cellDof = 0, numComponents = 0; f < Nf; ++f) { 803 PetscInt Nb, Nc; 804 805 ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr); 806 ierr = PetscFEGetNumComponents(feBd[f], &Nc);CHKERRQ(ierr); 807 cellDof += Nb*Nc; 808 numComponents += Nc; 809 } 810 for (p = 0, numFaces = 0; p < numPoints; ++p) { 811 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 812 if (dep == dim-1) ++numFaces; 813 } 814 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); 815 for (p = 0, f = 0; p < numPoints; ++p) { 816 const PetscInt point = points[p]; 817 PetscScalar *x = NULL; 818 PetscInt i; 819 820 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 821 if (dep != dim-1) continue; 822 ierr = DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr); 823 ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]); 824 if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point); 825 ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 826 for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i]; 827 ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 828 ++f; 829 } 830 for (f = 0; f < Nf; ++f) { 831 void (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdFuncs[f]; 832 void (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdFuncs[f]; 833 PetscInt numQuadPoints, Nb; 834 /* Conforming batches */ 835 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 836 /* Remainder */ 837 PetscInt Nr, offset; 838 839 ierr = PetscFEGetQuadrature(feBd[f], &q);CHKERRQ(ierr); 840 ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr); 841 ierr = PetscFEGetTileSizes(feBd[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 842 ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 843 blockSize = Nb*numQuadPoints; 844 batchSize = numBlocks * blockSize; 845 ierr = PetscFESetTileSizes(feBd[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 846 numChunks = numFaces / (numBatches*batchSize); 847 Ne = numChunks*numBatches*batchSize; 848 Nr = numFaces % (numBatches*batchSize); 849 offset = numFaces - Nr; 850 geom.v0 = v0; 851 geom.n = n; 852 geom.J = J; 853 geom.invJ = invJ; 854 geom.detJ = detJ; 855 ierr = PetscFEIntegrateBdResidual(feBd[f], Ne, Nf, feBd, f, geom, u, 0, NULL, NULL, f0, f1, elemVec);CHKERRQ(ierr); 856 geom.v0 = &v0[offset*dim]; 857 geom.n = &n[offset*dim]; 858 geom.J = &J[offset*dim*dim]; 859 geom.invJ = &invJ[offset*dim*dim]; 860 geom.detJ = &detJ[offset]; 861 ierr = PetscFEIntegrateBdResidual(feBd[f], Nr, Nf, feBd, f, geom, &u[offset*cellDof], 0, NULL, NULL, f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr); 862 } 863 for (p = 0, f = 0; p < numPoints; ++p) { 864 const PetscInt point = points[p]; 865 866 ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr); 867 if (dep != dim-1) continue; 868 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "BdResidual", cellDof, &elemVec[f*cellDof]);CHKERRQ(ierr);} 869 ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*cellDof], ADD_VALUES);CHKERRQ(ierr); 870 ++f; 871 } 872 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 873 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 874 ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr); 875 } 876 if (mesh->printFEM) {ierr = DMPrintLocalVec(dm, name, mesh->printTol, F);CHKERRQ(ierr);} 877 ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr); 878 PetscFunctionReturn(0); 879 } 880 881 #undef __FUNCT__ 882 #define __FUNCT__ "DMPlexComputeIFunctionFEM" 883 /*@ 884 DMPlexComputeIFunctionFEM - Form the local implicit function F from the local input X, X_t using pointwise functions specified by the user 885 886 Input Parameters: 887 + dm - The mesh 888 . time - The current time 889 . X - Local input vector 890 . X_t - Time derivative of the 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: DMPlexComputeResidualFEM() 905 @*/ 906 PetscErrorCode DMPlexComputeIFunctionFEM(DM dm, PetscReal time, Vec X, Vec X_t, 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, *u_t, *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 = 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); 958 if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);} 959 for (c = cStart; c < cEnd; ++c) { 960 PetscScalar *x = NULL, *x_t = 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 ierr = DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr); 969 for (i = 0; i < cellDof; ++i) u_t[c*cellDof+i] = x_t[i]; 970 ierr = DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr); 971 if (dmAux) { 972 PetscScalar *x_a = NULL; 973 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x_a);CHKERRQ(ierr); 974 for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x_a[i]; 975 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x_a);CHKERRQ(ierr); 976 } 977 } 978 for (f = 0; f < Nf; ++f) { 979 void (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0IFuncs[f]; 980 void (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1IFuncs[f]; 981 PetscInt numQuadPoints, Nb; 982 /* Conforming batches */ 983 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 984 /* Remainder */ 985 PetscInt Nr, offset; 986 987 ierr = PetscFEGetQuadrature(fe[f], &q);CHKERRQ(ierr); 988 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 989 ierr = PetscFEGetTileSizes(fe[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 990 ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 991 blockSize = Nb*numQuadPoints; 992 batchSize = numBlocks * blockSize; 993 ierr = PetscFESetTileSizes(fe[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 994 numChunks = numCells / (numBatches*batchSize); 995 Ne = numChunks*numBatches*batchSize; 996 Nr = numCells % (numBatches*batchSize); 997 offset = numCells - Nr; 998 geom.v0 = v0; 999 geom.J = J; 1000 geom.invJ = invJ; 1001 geom.detJ = detJ; 1002 ierr = PetscFEIntegrateIFunction(fe[f], Ne, Nf, fe, f, geom, u, u_t, NfAux, feAux, a, f0, f1, elemVec);CHKERRQ(ierr); 1003 geom.v0 = &v0[offset*dim]; 1004 geom.J = &J[offset*dim*dim]; 1005 geom.invJ = &invJ[offset*dim*dim]; 1006 geom.detJ = &detJ[offset]; 1007 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); 1008 } 1009 for (c = cStart; c < cEnd; ++c) { 1010 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);} 1011 ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr); 1012 } 1013 ierr = PetscFree7(u,u_t,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr); 1014 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1015 if (feBd) { 1016 DMLabel label, depth; 1017 IS pointIS; 1018 const PetscInt *points; 1019 PetscInt dep, numPoints, p, numFaces; 1020 PetscReal *n; 1021 1022 ierr = DMPlexGetLabel(dm, "boundary", &label);CHKERRQ(ierr); 1023 ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr); 1024 ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr); 1025 ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr); 1026 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 1027 for (f = 0, cellDof = 0, numComponents = 0; f < Nf; ++f) { 1028 PetscInt Nb, Nc; 1029 1030 ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr); 1031 ierr = PetscFEGetNumComponents(feBd[f], &Nc);CHKERRQ(ierr); 1032 cellDof += Nb*Nc; 1033 numComponents += Nc; 1034 } 1035 for (p = 0, numFaces = 0; p < numPoints; ++p) { 1036 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 1037 if (dep == dim-1) ++numFaces; 1038 } 1039 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); 1040 ierr = PetscMalloc1(numFaces*cellDof,&u_t);CHKERRQ(ierr); 1041 for (p = 0, f = 0; p < numPoints; ++p) { 1042 const PetscInt point = points[p]; 1043 PetscScalar *x = NULL; 1044 PetscInt i; 1045 1046 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 1047 if (dep != dim-1) continue; 1048 ierr = DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr); 1049 ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]); 1050 if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point); 1051 ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 1052 for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i]; 1053 ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 1054 ierr = DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr); 1055 for (i = 0; i < cellDof; ++i) u_t[f*cellDof+i] = x[i]; 1056 ierr = DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr); 1057 ++f; 1058 } 1059 for (f = 0; f < Nf; ++f) { 1060 void (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdIFuncs[f]; 1061 void (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdIFuncs[f]; 1062 PetscInt numQuadPoints, Nb; 1063 /* Conforming batches */ 1064 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 1065 /* Remainder */ 1066 PetscInt Nr, offset; 1067 1068 ierr = PetscFEGetQuadrature(feBd[f], &q);CHKERRQ(ierr); 1069 ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr); 1070 ierr = PetscFEGetTileSizes(feBd[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1071 ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 1072 blockSize = Nb*numQuadPoints; 1073 batchSize = numBlocks * blockSize; 1074 ierr = PetscFESetTileSizes(feBd[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1075 numChunks = numFaces / (numBatches*batchSize); 1076 Ne = numChunks*numBatches*batchSize; 1077 Nr = numFaces % (numBatches*batchSize); 1078 offset = numFaces - Nr; 1079 geom.v0 = v0; 1080 geom.n = n; 1081 geom.J = J; 1082 geom.invJ = invJ; 1083 geom.detJ = detJ; 1084 ierr = PetscFEIntegrateBdIFunction(feBd[f], Ne, Nf, feBd, f, geom, u, u_t, 0, NULL, NULL, f0, f1, elemVec);CHKERRQ(ierr); 1085 geom.v0 = &v0[offset*dim]; 1086 geom.n = &n[offset*dim]; 1087 geom.J = &J[offset*dim*dim]; 1088 geom.invJ = &invJ[offset*dim*dim]; 1089 geom.detJ = &detJ[offset]; 1090 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); 1091 } 1092 for (p = 0, f = 0; p < numPoints; ++p) { 1093 const PetscInt point = points[p]; 1094 1095 ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr); 1096 if (dep != dim-1) continue; 1097 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "BdResidual", cellDof, &elemVec[f*cellDof]);CHKERRQ(ierr);} 1098 ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*cellDof], ADD_VALUES);CHKERRQ(ierr); 1099 ++f; 1100 } 1101 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 1102 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 1103 ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr); 1104 ierr = PetscFree(u_t);CHKERRQ(ierr); 1105 } 1106 if (mesh->printFEM) {ierr = DMPrintLocalVec(dm, name, mesh->printTol, F);CHKERRQ(ierr);} 1107 ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr); 1108 PetscFunctionReturn(0); 1109 } 1110 1111 #undef __FUNCT__ 1112 #define __FUNCT__ "DMPlexComputeJacobianActionFEM" 1113 /*@C 1114 DMPlexComputeJacobianActionFEM - Form the local action of Jacobian J(u) on the local input X using pointwise functions specified by the user 1115 1116 Input Parameters: 1117 + dm - The mesh 1118 . J - The Jacobian shell matrix 1119 . X - Local input vector 1120 - user - The user context 1121 1122 Output Parameter: 1123 . F - Local output vector 1124 1125 Note: 1126 The first member of the user context must be an FEMContext. 1127 1128 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 1129 like a GPU, or vectorize on a multicore machine. 1130 1131 Level: developer 1132 1133 .seealso: DMPlexComputeResidualFEM() 1134 @*/ 1135 PetscErrorCode DMPlexComputeJacobianActionFEM(DM dm, Mat Jac, Vec X, Vec F, void *user) 1136 { 1137 DM_Plex *mesh = (DM_Plex *) dm->data; 1138 PetscFEM *fem = (PetscFEM *) user; 1139 PetscFE *fe = fem->fe; 1140 PetscQuadrature quad; 1141 PetscCellGeometry geom; 1142 PetscSection section; 1143 JacActionCtx *jctx; 1144 PetscReal *v0, *J, *invJ, *detJ; 1145 PetscScalar *elemVec, *u, *a; 1146 PetscInt dim, numFields, field, numCells, cStart, cEnd, c; 1147 PetscInt cellDof = 0; 1148 PetscErrorCode ierr; 1149 1150 PetscFunctionBegin; 1151 /* ierr = PetscLogEventBegin(DMPLEX_JacobianActionFEM,dm,0,0,0);CHKERRQ(ierr); */ 1152 ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr); 1153 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1154 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1155 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1156 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1157 numCells = cEnd - cStart; 1158 for (field = 0; field < numFields; ++field) { 1159 PetscInt Nb, Nc; 1160 1161 ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr); 1162 ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr); 1163 cellDof += Nb*Nc; 1164 } 1165 ierr = VecSet(F, 0.0);CHKERRQ(ierr); 1166 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); 1167 for (c = cStart; c < cEnd; ++c) { 1168 PetscScalar *x = NULL; 1169 PetscInt i; 1170 1171 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 1172 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 1173 ierr = DMPlexVecGetClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr); 1174 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 1175 ierr = DMPlexVecRestoreClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr); 1176 ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 1177 for (i = 0; i < cellDof; ++i) a[c*cellDof+i] = x[i]; 1178 ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 1179 } 1180 for (field = 0; field < numFields; ++field) { 1181 PetscInt numQuadPoints, Nb; 1182 /* Conforming batches */ 1183 PetscInt numBlocks = 1; 1184 PetscInt numBatches = 1; 1185 PetscInt numChunks, Ne, blockSize, batchSize; 1186 /* Remainder */ 1187 PetscInt Nr, offset; 1188 1189 ierr = PetscFEGetQuadrature(fe[field], &quad);CHKERRQ(ierr); 1190 ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr); 1191 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 1192 blockSize = Nb*numQuadPoints; 1193 batchSize = numBlocks * blockSize; 1194 numChunks = numCells / (numBatches*batchSize); 1195 Ne = numChunks*numBatches*batchSize; 1196 Nr = numCells % (numBatches*batchSize); 1197 offset = numCells - Nr; 1198 geom.v0 = v0; 1199 geom.J = J; 1200 geom.invJ = invJ; 1201 geom.detJ = detJ; 1202 ierr = PetscFEIntegrateJacobianAction(fe[field], Ne, numFields, fe, field, geom, u, a, fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, elemVec);CHKERRQ(ierr); 1203 geom.v0 = &v0[offset*dim]; 1204 geom.J = &J[offset*dim*dim]; 1205 geom.invJ = &invJ[offset*dim*dim]; 1206 geom.detJ = &detJ[offset]; 1207 ierr = PetscFEIntegrateJacobianAction(fe[field], Nr, numFields, fe, field, geom, &u[offset*cellDof], &a[offset*cellDof], 1208 fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, &elemVec[offset*cellDof]);CHKERRQ(ierr); 1209 } 1210 for (c = cStart; c < cEnd; ++c) { 1211 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Jacobian Action", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);} 1212 ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr); 1213 } 1214 ierr = PetscFree7(u,a,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr); 1215 if (mesh->printFEM) { 1216 PetscMPIInt rank, numProcs; 1217 PetscInt p; 1218 1219 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 1220 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr); 1221 ierr = PetscPrintf(PetscObjectComm((PetscObject)dm), "Jacobian Action:\n");CHKERRQ(ierr); 1222 for (p = 0; p < numProcs; ++p) { 1223 if (p == rank) {ierr = VecView(F, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);} 1224 ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr); 1225 } 1226 } 1227 /* ierr = PetscLogEventEnd(DMPLEX_JacobianActionFEM,dm,0,0,0);CHKERRQ(ierr); */ 1228 PetscFunctionReturn(0); 1229 } 1230 1231 #undef __FUNCT__ 1232 #define __FUNCT__ "DMPlexComputeJacobianFEM" 1233 /*@ 1234 DMPlexComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user. 1235 1236 Input Parameters: 1237 + dm - The mesh 1238 . X - Local input vector 1239 - user - The user context 1240 1241 Output Parameter: 1242 . Jac - Jacobian matrix 1243 1244 Note: 1245 The first member of the user context must be an FEMContext. 1246 1247 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 1248 like a GPU, or vectorize on a multicore machine. 1249 1250 Level: developer 1251 1252 .seealso: FormFunctionLocal() 1253 @*/ 1254 PetscErrorCode DMPlexComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user) 1255 { 1256 DM_Plex *mesh = (DM_Plex *) dm->data; 1257 PetscFEM *fem = (PetscFEM *) user; 1258 PetscFE *fe = fem->fe; 1259 PetscFE *feAux = fem->feAux; 1260 const char *name = "Jacobian"; 1261 DM dmAux; 1262 Vec A; 1263 PetscQuadrature quad; 1264 PetscCellGeometry geom; 1265 PetscSection section, globalSection, sectionAux; 1266 PetscReal *v0, *J, *invJ, *detJ; 1267 PetscScalar *elemMat, *u, *a; 1268 PetscInt dim, Nf, NfAux = 0, f, fieldI, fieldJ, numCells, cStart, cEnd, c; 1269 PetscInt cellDof = 0, numComponents = 0; 1270 PetscInt cellDofAux = 0, numComponentsAux = 0; 1271 PetscBool isShell; 1272 PetscErrorCode ierr; 1273 1274 PetscFunctionBegin; 1275 ierr = PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr); 1276 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1277 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1278 ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr); 1279 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1280 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1281 numCells = cEnd - cStart; 1282 for (f = 0; f < Nf; ++f) { 1283 PetscInt Nb, Nc; 1284 1285 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 1286 ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr); 1287 cellDof += Nb*Nc; 1288 numComponents += Nc; 1289 } 1290 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1291 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1292 if (dmAux) { 1293 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1294 ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr); 1295 } 1296 for (f = 0; f < NfAux; ++f) { 1297 PetscInt Nb, Nc; 1298 1299 ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr); 1300 ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr); 1301 cellDofAux += Nb*Nc; 1302 numComponentsAux += Nc; 1303 } 1304 ierr = DMPlexInsertBoundaryValuesFEM(dm, X);CHKERRQ(ierr); 1305 ierr = MatZeroEntries(JacP);CHKERRQ(ierr); 1306 ierr = PetscMalloc6(numCells*cellDof,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof*cellDof,&elemMat);CHKERRQ(ierr); 1307 if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);} 1308 for (c = cStart; c < cEnd; ++c) { 1309 PetscScalar *x = NULL; 1310 PetscInt i; 1311 1312 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 1313 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 1314 ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 1315 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 1316 ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 1317 if (dmAux) { 1318 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1319 for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x[i]; 1320 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1321 } 1322 } 1323 ierr = PetscMemzero(elemMat, numCells*cellDof*cellDof * sizeof(PetscScalar));CHKERRQ(ierr); 1324 for (fieldI = 0; fieldI < Nf; ++fieldI) { 1325 PetscInt numQuadPoints, Nb; 1326 /* Conforming batches */ 1327 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 1328 /* Remainder */ 1329 PetscInt Nr, offset; 1330 1331 ierr = PetscFEGetQuadrature(fe[fieldI], &quad);CHKERRQ(ierr); 1332 ierr = PetscFEGetDimension(fe[fieldI], &Nb);CHKERRQ(ierr); 1333 ierr = PetscFEGetTileSizes(fe[fieldI], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1334 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 1335 blockSize = Nb*numQuadPoints; 1336 batchSize = numBlocks * blockSize; 1337 ierr = PetscFESetTileSizes(fe[fieldI], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1338 numChunks = numCells / (numBatches*batchSize); 1339 Ne = numChunks*numBatches*batchSize; 1340 Nr = numCells % (numBatches*batchSize); 1341 offset = numCells - Nr; 1342 for (fieldJ = 0; fieldJ < Nf; ++fieldJ) { 1343 void (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g0Funcs[fieldI*Nf+fieldJ]; 1344 void (*g1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g1Funcs[fieldI*Nf+fieldJ]; 1345 void (*g2)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g2Funcs[fieldI*Nf+fieldJ]; 1346 void (*g3)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g3Funcs[fieldI*Nf+fieldJ]; 1347 1348 geom.v0 = v0; 1349 geom.J = J; 1350 geom.invJ = invJ; 1351 geom.detJ = detJ; 1352 ierr = PetscFEIntegrateJacobian(fe[fieldI], Ne, Nf, fe, fieldI, fieldJ, geom, u, NfAux, feAux, a, g0, g1, g2, g3, elemMat);CHKERRQ(ierr); 1353 geom.v0 = &v0[offset*dim]; 1354 geom.J = &J[offset*dim*dim]; 1355 geom.invJ = &invJ[offset*dim*dim]; 1356 geom.detJ = &detJ[offset]; 1357 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); 1358 } 1359 } 1360 for (c = cStart; c < cEnd; ++c) { 1361 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, name, cellDof, cellDof, &elemMat[c*cellDof*cellDof]);CHKERRQ(ierr);} 1362 ierr = DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMat[c*cellDof*cellDof], ADD_VALUES);CHKERRQ(ierr); 1363 } 1364 ierr = PetscFree6(u,v0,J,invJ,detJ,elemMat);CHKERRQ(ierr); 1365 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1366 ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1367 ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1368 if (mesh->printFEM) { 1369 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1370 ierr = MatChop(JacP, 1.0e-10);CHKERRQ(ierr); 1371 ierr = MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1372 } 1373 ierr = PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr); 1374 ierr = PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);CHKERRQ(ierr); 1375 if (isShell) { 1376 JacActionCtx *jctx; 1377 1378 ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr); 1379 ierr = VecCopy(X, jctx->u);CHKERRQ(ierr); 1380 } 1381 PetscFunctionReturn(0); 1382 } 1383 1384 #if 0 1385 1386 static void g0_identity_1d_static(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[]) 1387 { 1388 g3[0] = 1.0; 1389 } 1390 1391 static void g0_identity_2d_static(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[]) 1392 { 1393 g3[0] = g3[3] = 1.0; 1394 } 1395 1396 static void g0_identity_3d_static(const PetscScalar u[], const PetscScalar gradU[], const PetscScalar a[], const PetscScalar gradA[], const PetscReal x[], PetscScalar g3[]) 1397 { 1398 g3[0] = g3[4] = g3[8] = 1.0; 1399 } 1400 1401 #undef __FUNCT__ 1402 #define __FUNCT__ "DMPlexComputeInterpolatorFEMBroken" 1403 PetscErrorCode DMPlexComputeInterpolatorFEMBroken(DM dmc, DM dmf, Mat I, void *user) 1404 { 1405 DM_Plex *mesh = (DM_Plex *) dmc->data; 1406 PetscFEM *fem = (PetscFEM *) user; 1407 PetscFE *fe = fem->fe; 1408 const char *name = "Interpolator"; 1409 PetscFE *feRef; 1410 PetscQuadrature quad, quadOld; 1411 PetscCellGeometry geom; 1412 PetscSection fsection, fglobalSection; 1413 PetscSection csection, cglobalSection; 1414 PetscReal *v0, *J, *invJ, *detJ; 1415 PetscScalar *elemMat; 1416 PetscInt dim, Nf, f, fieldI, fieldJ, numCells, cStart, cEnd, c; 1417 PetscInt rCellDof = 0, cCellDof = 0, numComponents = 0; 1418 PetscErrorCode ierr; 1419 1420 PetscFunctionBegin; 1421 #if 0 1422 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1423 #endif 1424 ierr = DMPlexGetDimension(dmf, &dim);CHKERRQ(ierr); 1425 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1426 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1427 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1428 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1429 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1430 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1431 numCells = cEnd - cStart; 1432 ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr); 1433 for (fieldI = 0; fieldI < Nf; ++fieldI) { 1434 ierr = PetscFERefine(fe[fieldI], &feRef[fieldI]);CHKERRQ(ierr); 1435 } 1436 for (f = 0; f < Nf; ++f) { 1437 PetscInt rNb, cNb, Nc; 1438 1439 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1440 ierr = PetscFEGetDimension(fe[f], &cNb);CHKERRQ(ierr); 1441 ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr); 1442 numComponents += Nc; 1443 rCellDof += rNb*Nc; 1444 cCellDof += cNb*Nc; 1445 } 1446 ierr = MatZeroEntries(I);CHKERRQ(ierr); 1447 ierr = PetscMalloc5(numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*rCellDof*cCellDof,&elemMat);CHKERRQ(ierr); 1448 for (c = cStart; c < cEnd; ++c) { 1449 ierr = DMPlexComputeCellGeometry(dmc, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 1450 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 1451 } 1452 ierr = PetscMemzero(elemMat, numCells*rCellDof*cCellDof * sizeof(PetscScalar));CHKERRQ(ierr); 1453 for (fieldI = 0; fieldI < Nf; ++fieldI) { 1454 PetscInt numQuadPoints, Nb; 1455 /* Conforming batches */ 1456 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 1457 /* Remainder */ 1458 PetscInt Nr, offset; 1459 1460 /* Make new fine FE which refines the ref cell and the quadrature rule */ 1461 ierr = PetscFEGetQuadrature(feRef[fieldI], &quad);CHKERRQ(ierr); 1462 ierr = PetscFEGetDimension(feRef[fieldI], &Nb);CHKERRQ(ierr); 1463 ierr = PetscFEGetTileSizes(feRef[fieldI], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1464 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 1465 blockSize = Nb*numQuadPoints; 1466 batchSize = numBlocks * blockSize; 1467 ierr = PetscFESetTileSizes(feRef[fieldI], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1468 numChunks = numCells / (numBatches*batchSize); 1469 Ne = numChunks*numBatches*batchSize; 1470 Nr = numCells % (numBatches*batchSize); 1471 offset = numCells - Nr; 1472 1473 for (fieldJ = 0; fieldJ < Nf; ++fieldJ) { 1474 /* void (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g0Funcs[fieldI*Nf+fieldJ]; */ 1475 void (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = g0_identity_2d_static; 1476 1477 /* Replace quadrature in coarse FE with refined quadrature */ 1478 ierr = PetscFEGetQuadrature(fe[fieldJ], &quadOld);CHKERRQ(ierr); 1479 ierr = PetscObjectReference((PetscObject) quadOld);CHKERRQ(ierr); 1480 ierr = PetscFESetQuadrature(fe[fieldJ], quad);CHKERRQ(ierr); 1481 geom.v0 = v0; 1482 geom.J = J; 1483 geom.invJ = invJ; 1484 geom.detJ = detJ; 1485 ierr = PetscFEIntegrateInterpolator_Basic(feRef[fieldI], Ne, Nf, feRef, fieldI, fe, fieldJ, geom, g0, elemMat);CHKERRQ(ierr); 1486 geom.v0 = &v0[offset*dim]; 1487 geom.J = &J[offset*dim*dim]; 1488 geom.invJ = &invJ[offset*dim*dim]; 1489 geom.detJ = &detJ[offset]; 1490 ierr = PetscFEIntegrateInterpolator_Basic(feRef[fieldI], Nr, Nf, feRef, fieldI, fe, fieldJ, geom, g0, &elemMat[offset*rCellDof*cCellDof]);CHKERRQ(ierr); 1491 ierr = PetscFESetQuadrature(fe[fieldJ], quadOld);CHKERRQ(ierr); 1492 } 1493 } 1494 for (c = cStart; c < cEnd; ++c) { 1495 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, name, rCellDof, cCellDof, &elemMat[c*rCellDof*cCellDof]);CHKERRQ(ierr);} 1496 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, I, c, &elemMat[c*rCellDof*cCellDof], ADD_VALUES);CHKERRQ(ierr); 1497 } 1498 ierr = PetscFree5(v0,J,invJ,detJ,elemMat);CHKERRQ(ierr); 1499 ierr = PetscFree(feRef);CHKERRQ(ierr); 1500 ierr = MatAssemblyBegin(I, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1501 ierr = MatAssemblyEnd(I, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1502 if (mesh->printFEM) { 1503 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1504 ierr = MatChop(I, 1.0e-10);CHKERRQ(ierr); 1505 ierr = MatView(I, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1506 } 1507 #if 0 1508 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1509 #endif 1510 PetscFunctionReturn(0); 1511 } 1512 #endif 1513 1514 #undef __FUNCT__ 1515 #define __FUNCT__ "DMPlexComputeInterpolatorFEM" 1516 /*@ 1517 DMPlexComputeInterpolatorFEM - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1518 1519 Input Parameters: 1520 + dmf - The fine mesh 1521 . dmc - The coarse mesh 1522 - user - The user context 1523 1524 Output Parameter: 1525 . In - The interpolation matrix 1526 1527 Note: 1528 The first member of the user context must be an FEMContext. 1529 1530 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 1531 like a GPU, or vectorize on a multicore machine. 1532 1533 Level: developer 1534 1535 .seealso: DMPlexComputeJacobianFEM() 1536 @*/ 1537 PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user) 1538 { 1539 DM_Plex *mesh = (DM_Plex *) dmc->data; 1540 PetscFEM *fem = (PetscFEM *) user; 1541 PetscFE *fe = fem->fe; 1542 const char *name = "Interpolator"; 1543 PetscFE *feRef; 1544 PetscSection fsection, fglobalSection; 1545 PetscSection csection, cglobalSection; 1546 PetscScalar *elemMat; 1547 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, c; 1548 PetscInt rCellDof = 0, cCellDof = 0; 1549 PetscErrorCode ierr; 1550 1551 PetscFunctionBegin; 1552 #if 0 1553 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1554 #endif 1555 ierr = DMPlexGetDimension(dmf, &dim);CHKERRQ(ierr); 1556 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1557 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1558 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1559 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1560 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1561 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1562 ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr); 1563 for (f = 0; f < Nf; ++f) { 1564 PetscInt rNb, cNb, Nc; 1565 1566 ierr = PetscFERefine(fe[f], &feRef[f]);CHKERRQ(ierr); 1567 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1568 ierr = PetscFEGetDimension(fe[f], &cNb);CHKERRQ(ierr); 1569 ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr); 1570 rCellDof += rNb*Nc; 1571 cCellDof += cNb*Nc; 1572 } 1573 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1574 ierr = PetscMalloc1(rCellDof*cCellDof,&elemMat);CHKERRQ(ierr); 1575 ierr = PetscMemzero(elemMat, rCellDof*cCellDof * sizeof(PetscScalar));CHKERRQ(ierr); 1576 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1577 PetscDualSpace Qref; 1578 PetscQuadrature f; 1579 const PetscReal *qpoints, *qweights; 1580 PetscReal *points; 1581 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1582 1583 /* Compose points from all dual basis functionals */ 1584 ierr = PetscFEGetNumComponents(fe[fieldI], &Nc);CHKERRQ(ierr); 1585 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1586 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1587 for (i = 0; i < fpdim; ++i) { 1588 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1589 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1590 npoints += Np; 1591 } 1592 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1593 for (i = 0, k = 0; i < fpdim; ++i) { 1594 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1595 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1596 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1597 } 1598 1599 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1600 PetscReal *B; 1601 PetscInt NcJ, cpdim, j; 1602 1603 /* Evaluate basis at points */ 1604 ierr = PetscFEGetNumComponents(fe[fieldJ], &NcJ);CHKERRQ(ierr); 1605 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); 1606 ierr = PetscFEGetDimension(fe[fieldJ], &cpdim);CHKERRQ(ierr); 1607 /* For now, fields only interpolate themselves */ 1608 if (fieldI == fieldJ) { 1609 ierr = PetscFEGetTabulation(fe[fieldJ], npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1610 for (i = 0, k = 0; i < fpdim; ++i) { 1611 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1612 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1613 for (p = 0; p < Np; ++p, ++k) { 1614 for (j = 0; j < cpdim; ++j) { 1615 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]; 1616 } 1617 } 1618 } 1619 ierr = PetscFERestoreTabulation(fe[fieldJ], npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1620 } 1621 offsetJ += cpdim*NcJ; 1622 } 1623 offsetI += fpdim*Nc; 1624 ierr = PetscFree(points);CHKERRQ(ierr); 1625 } 1626 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rCellDof, cCellDof, elemMat);CHKERRQ(ierr);} 1627 for (c = cStart; c < cEnd; ++c) { 1628 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1629 } 1630 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1631 ierr = PetscFree(feRef);CHKERRQ(ierr); 1632 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1633 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1634 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1635 if (mesh->printFEM) { 1636 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1637 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1638 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1639 } 1640 #if 0 1641 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1642 #endif 1643 PetscFunctionReturn(0); 1644 } 1645 1646 #undef __FUNCT__ 1647 #define __FUNCT__ "DMPlexAddBoundary" 1648 /* The ids can be overridden by the command line option -bc_<boundary name> */ 1649 PetscErrorCode DMPlexAddBoundary(DM dm, PetscBool isEssential, const char name[], PetscInt field, void (*bcFunc)(), PetscInt numids, const PetscInt *ids, void *ctx) 1650 { 1651 DM_Plex *mesh = (DM_Plex *) dm->data; 1652 DMBoundary b; 1653 PetscErrorCode ierr; 1654 1655 PetscFunctionBegin; 1656 ierr = PetscNew(&b);CHKERRQ(ierr); 1657 ierr = PetscStrallocpy(name, (char **) &b->name);CHKERRQ(ierr); 1658 ierr = PetscMalloc1(numids, &b->ids);CHKERRQ(ierr); 1659 ierr = PetscMemcpy(b->ids, ids, numids*sizeof(PetscInt));CHKERRQ(ierr); 1660 b->essential = isEssential; 1661 b->field = field; 1662 b->func = bcFunc; 1663 b->numids = numids; 1664 b->ctx = ctx; 1665 b->next = mesh->boundary; 1666 mesh->boundary = b; 1667 PetscFunctionReturn(0); 1668 } 1669 1670 #undef __FUNCT__ 1671 #define __FUNCT__ "DMPlexGetNumBoundary" 1672 PetscErrorCode DMPlexGetNumBoundary(DM dm, PetscInt *numBd) 1673 { 1674 DM_Plex *mesh = (DM_Plex *) dm->data; 1675 DMBoundary b = mesh->boundary; 1676 1677 PetscFunctionBegin; 1678 *numBd = 0; 1679 while (b) {++(*numBd); b = b->next;} 1680 PetscFunctionReturn(0); 1681 } 1682 1683 #undef __FUNCT__ 1684 #define __FUNCT__ "DMPlexGetBoundary" 1685 PetscErrorCode DMPlexGetBoundary(DM dm, PetscInt bd, PetscBool *isEssential, const char **name, PetscInt *field, void (**func)(), PetscInt *numids, const PetscInt **ids, void **ctx) 1686 { 1687 DM_Plex *mesh = (DM_Plex *) dm->data; 1688 DMBoundary b = mesh->boundary; 1689 PetscInt n = 0; 1690 1691 PetscFunctionBegin; 1692 while (b) { 1693 if (n == bd) break; 1694 b = b->next; 1695 ++n; 1696 } 1697 if (n != bd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n); 1698 if (isEssential) { 1699 PetscValidPointer(isEssential, 3); 1700 *isEssential = b->essential; 1701 } 1702 if (name) { 1703 PetscValidPointer(name, 4); 1704 *name = b->name; 1705 } 1706 if (field) { 1707 PetscValidPointer(field, 5); 1708 *field = b->field; 1709 } 1710 if (func) { 1711 PetscValidPointer(func, 6); 1712 *func = b->func; 1713 } 1714 if (numids) { 1715 PetscValidPointer(numids, 7); 1716 *numids = b->numids; 1717 } 1718 if (ids) { 1719 PetscValidPointer(ids, 8); 1720 *ids = b->ids; 1721 } 1722 if (ctx) { 1723 PetscValidPointer(ctx, 9); 1724 *ctx = b->ctx; 1725 } 1726 PetscFunctionReturn(0); 1727 } 1728