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