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