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 *labelname; 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, NULL, &labelname, &field, &func, &numids, &ids, &ctx);CHKERRQ(ierr); 368 ierr = DMPlexGetLabel(dm, labelname, &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__ "DMPlexComputeL2Diff" 381 /*@C 382 DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h. 383 384 Input Parameters: 385 + dm - The DM 386 . fe - The PetscFE object for each field 387 . funcs - The functions to evaluate for each field component 388 . ctxs - Optional array of contexts to pass to each function, or NULL. 389 - X - The coefficient vector u_h 390 391 Output Parameter: 392 . diff - The diff ||u - u_h||_2 393 394 Level: developer 395 396 .seealso: DMPlexProjectFunction(), DMPlexComputeL2GradientDiff() 397 @*/ 398 PetscErrorCode DMPlexComputeL2Diff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 399 { 400 const PetscInt debug = 0; 401 PetscSection section; 402 PetscQuadrature quad; 403 Vec localX; 404 PetscScalar *funcVal; 405 PetscReal *coords, *v0, *J, *invJ, detJ; 406 PetscReal localDiff = 0.0; 407 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp; 408 PetscErrorCode ierr; 409 410 PetscFunctionBegin; 411 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 412 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 413 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 414 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 415 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 416 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 417 for (field = 0; field < numFields; ++field) { 418 PetscInt Nc; 419 420 ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr); 421 numComponents += Nc; 422 } 423 ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 424 ierr = PetscMalloc5(numComponents,&funcVal,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 425 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 426 ierr = PetscFEGetQuadrature(fe[0], &quad);CHKERRQ(ierr); 427 for (c = cStart; c < cEnd; ++c) { 428 PetscScalar *x = NULL; 429 PetscReal elemDiff = 0.0; 430 431 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr); 432 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 433 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 434 435 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 436 void * const ctx = ctxs ? ctxs[field] : NULL; 437 const PetscReal *quadPoints, *quadWeights; 438 PetscReal *basis; 439 PetscInt numQuadPoints, numBasisFuncs, numBasisComps, q, d, e, fc, f; 440 441 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 442 ierr = PetscFEGetDimension(fe[field], &numBasisFuncs);CHKERRQ(ierr); 443 ierr = PetscFEGetNumComponents(fe[field], &numBasisComps);CHKERRQ(ierr); 444 ierr = PetscFEGetDefaultTabulation(fe[field], &basis, NULL, NULL);CHKERRQ(ierr); 445 if (debug) { 446 char title[1024]; 447 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 448 ierr = DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);CHKERRQ(ierr); 449 } 450 for (q = 0; q < numQuadPoints; ++q) { 451 for (d = 0; d < dim; d++) { 452 coords[d] = v0[d]; 453 for (e = 0; e < dim; e++) { 454 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 455 } 456 } 457 (*funcs[field])(coords, funcVal, ctx); 458 for (fc = 0; fc < numBasisComps; ++fc) { 459 PetscScalar interpolant = 0.0; 460 461 for (f = 0; f < numBasisFuncs; ++f) { 462 const PetscInt fidx = f*numBasisComps+fc; 463 interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx]; 464 } 465 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);} 466 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 467 } 468 } 469 comp += numBasisComps; 470 fieldOffset += numBasisFuncs*numBasisComps; 471 } 472 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 473 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 474 localDiff += elemDiff; 475 } 476 ierr = PetscFree5(funcVal,coords,v0,J,invJ);CHKERRQ(ierr); 477 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 478 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 479 *diff = PetscSqrtReal(*diff); 480 PetscFunctionReturn(0); 481 } 482 483 #undef __FUNCT__ 484 #define __FUNCT__ "DMPlexComputeL2GradientDiff" 485 /*@C 486 DMPlexComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h. 487 488 Input Parameters: 489 + dm - The DM 490 . fe - The PetscFE object for each field 491 . funcs - The gradient functions to evaluate for each field component 492 . ctxs - Optional array of contexts to pass to each function, or NULL. 493 . X - The coefficient vector u_h 494 - n - The vector to project along 495 496 Output Parameter: 497 . diff - The diff ||(grad u - grad u_h) . n||_2 498 499 Level: developer 500 501 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff() 502 @*/ 503 PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, PetscFE fe[], void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff) 504 { 505 const PetscInt debug = 0; 506 PetscSection section; 507 PetscQuadrature quad; 508 Vec localX; 509 PetscScalar *funcVal, *interpolantVec; 510 PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ; 511 PetscReal localDiff = 0.0; 512 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp; 513 PetscErrorCode ierr; 514 515 PetscFunctionBegin; 516 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 517 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 518 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 519 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 520 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 521 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 522 for (field = 0; field < numFields; ++field) { 523 PetscInt Nc; 524 525 ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr); 526 numComponents += Nc; 527 } 528 /* ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 529 ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr); 530 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 531 ierr = PetscFEGetQuadrature(fe[0], &quad);CHKERRQ(ierr); 532 for (c = cStart; c < cEnd; ++c) { 533 PetscScalar *x = NULL; 534 PetscReal elemDiff = 0.0; 535 536 ierr = DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);CHKERRQ(ierr); 537 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 538 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 539 540 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 541 void * const ctx = ctxs ? ctxs[field] : NULL; 542 const PetscReal *quadPoints, *quadWeights; 543 PetscReal *basisDer; 544 PetscInt numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g; 545 546 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 547 ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr); 548 ierr = PetscFEGetNumComponents(fe[field], &Ncomp);CHKERRQ(ierr); 549 ierr = PetscFEGetDefaultTabulation(fe[field], NULL, &basisDer, NULL);CHKERRQ(ierr); 550 if (debug) { 551 char title[1024]; 552 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 553 ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr); 554 } 555 for (q = 0; q < numQuadPoints; ++q) { 556 for (d = 0; d < dim; d++) { 557 coords[d] = v0[d]; 558 for (e = 0; e < dim; e++) { 559 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 560 } 561 } 562 (*funcs[field])(coords, n, funcVal, ctx); 563 for (fc = 0; fc < Ncomp; ++fc) { 564 PetscScalar interpolant = 0.0; 565 566 for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0; 567 for (f = 0; f < Nb; ++f) { 568 const PetscInt fidx = f*Ncomp+fc; 569 570 for (d = 0; d < dim; ++d) { 571 realSpaceDer[d] = 0.0; 572 for (g = 0; g < dim; ++g) { 573 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g]; 574 } 575 interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d]; 576 } 577 } 578 for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d]; 579 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);} 580 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 581 } 582 } 583 comp += Ncomp; 584 fieldOffset += Nb*Ncomp; 585 } 586 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 587 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 588 localDiff += elemDiff; 589 } 590 ierr = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr); 591 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 592 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 593 *diff = PetscSqrtReal(*diff); 594 PetscFunctionReturn(0); 595 } 596 597 #undef __FUNCT__ 598 #define __FUNCT__ "DMPlexComputeResidualFEM" 599 /*@ 600 DMPlexComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user 601 602 Input Parameters: 603 + dm - The mesh 604 . X - Local input vector 605 - user - The user context 606 607 Output Parameter: 608 . F - Local output vector 609 610 Note: 611 The first member of the user context must be an FEMContext. 612 613 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 614 like a GPU, or vectorize on a multicore machine. 615 616 Level: developer 617 618 .seealso: DMPlexComputeJacobianActionFEM() 619 @*/ 620 PetscErrorCode DMPlexComputeResidualFEM(DM dm, Vec X, Vec F, void *user) 621 { 622 DM_Plex *mesh = (DM_Plex *) dm->data; 623 PetscFEM *fem = (PetscFEM *) user; 624 PetscFE *fe = fem->fe; 625 PetscFE *feAux = fem->feAux; 626 PetscFE *feBd = fem->feBd; 627 const char *name = "Residual"; 628 DM dmAux; 629 Vec A; 630 PetscQuadrature q; 631 PetscCellGeometry geom; 632 PetscSection section, sectionAux; 633 PetscReal *v0, *J, *invJ, *detJ; 634 PetscScalar *elemVec, *u, *a = NULL; 635 PetscInt dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, c; 636 PetscInt cellDof = 0, numComponents = 0; 637 PetscInt cellDofAux = 0, numComponentsAux = 0; 638 PetscErrorCode ierr; 639 640 PetscFunctionBegin; 641 ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr); 642 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 643 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 644 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 645 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 646 numCells = cEnd - cStart; 647 for (f = 0; f < Nf; ++f) { 648 PetscInt Nb, Nc; 649 650 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 651 ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr); 652 cellDof += Nb*Nc; 653 numComponents += Nc; 654 } 655 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 656 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 657 if (dmAux) { 658 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 659 ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr); 660 } 661 for (f = 0; f < NfAux; ++f) { 662 PetscInt Nb, Nc; 663 664 ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr); 665 ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr); 666 cellDofAux += Nb*Nc; 667 numComponentsAux += Nc; 668 } 669 ierr = DMPlexInsertBoundaryValuesFEM(dm, X);CHKERRQ(ierr); 670 ierr = VecSet(F, 0.0);CHKERRQ(ierr); 671 ierr = PetscMalloc6(numCells*cellDof,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof,&elemVec);CHKERRQ(ierr); 672 if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);} 673 for (c = cStart; c < cEnd; ++c) { 674 PetscScalar *x = NULL; 675 PetscInt i; 676 677 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 678 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 679 ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 680 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 681 ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 682 if (dmAux) { 683 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 684 for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x[i]; 685 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 686 } 687 } 688 for (f = 0; f < Nf; ++f) { 689 void (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0Funcs[f]; 690 void (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1Funcs[f]; 691 PetscInt numQuadPoints, Nb; 692 /* Conforming batches */ 693 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 694 /* Remainder */ 695 PetscInt Nr, offset; 696 697 ierr = PetscFEGetQuadrature(fe[f], &q);CHKERRQ(ierr); 698 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 699 ierr = PetscFEGetTileSizes(fe[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 700 ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 701 blockSize = Nb*numQuadPoints; 702 batchSize = numBlocks * blockSize; 703 ierr = PetscFESetTileSizes(fe[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 704 numChunks = numCells / (numBatches*batchSize); 705 Ne = numChunks*numBatches*batchSize; 706 Nr = numCells % (numBatches*batchSize); 707 offset = numCells - Nr; 708 geom.v0 = v0; 709 geom.J = J; 710 geom.invJ = invJ; 711 geom.detJ = detJ; 712 ierr = PetscFEIntegrateResidual(fe[f], Ne, Nf, fe, f, geom, u, NfAux, feAux, a, f0, f1, elemVec);CHKERRQ(ierr); 713 geom.v0 = &v0[offset*dim]; 714 geom.J = &J[offset*dim*dim]; 715 geom.invJ = &invJ[offset*dim*dim]; 716 geom.detJ = &detJ[offset]; 717 ierr = PetscFEIntegrateResidual(fe[f], Nr, Nf, fe, f, geom, &u[offset*cellDof], NfAux, feAux, &a[offset*cellDofAux], f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr); 718 } 719 for (c = cStart; c < cEnd; ++c) { 720 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);} 721 ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr); 722 } 723 ierr = PetscFree6(u,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr); 724 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 725 if (feBd) { 726 DMLabel depth; 727 PetscInt numBd, bd; 728 729 for (f = 0, cellDof = 0, numComponents = 0; f < Nf; ++f) { 730 PetscInt Nb, Nc; 731 732 ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr); 733 ierr = PetscFEGetNumComponents(feBd[f], &Nc);CHKERRQ(ierr); 734 cellDof += Nb*Nc; 735 numComponents += Nc; 736 } 737 ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr); 738 ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 739 for (bd = 0; bd < numBd; ++bd) { 740 const char *bdLabel; 741 DMLabel label; 742 IS pointIS; 743 const PetscInt *points; 744 const PetscInt *values; 745 PetscReal *n; 746 PetscInt field, numValues, numPoints, p, dep, numFaces; 747 PetscBool isEssential; 748 749 ierr = DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, &numValues, &values, NULL);CHKERRQ(ierr); 750 if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this"); 751 ierr = DMPlexGetLabel(dm, bdLabel, &label);CHKERRQ(ierr); 752 ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr); 753 ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr); 754 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 755 for (p = 0, numFaces = 0; p < numPoints; ++p) { 756 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 757 if (dep == dim-1) ++numFaces; 758 } 759 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); 760 for (p = 0, f = 0; p < numPoints; ++p) { 761 const PetscInt point = points[p]; 762 PetscScalar *x = NULL; 763 PetscInt i; 764 765 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 766 if (dep != dim-1) continue; 767 ierr = DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr); 768 ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]); 769 if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point); 770 ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 771 for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i]; 772 ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 773 ++f; 774 } 775 for (f = 0; f < Nf; ++f) { 776 void (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdFuncs[f]; 777 void (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdFuncs[f]; 778 PetscInt numQuadPoints, Nb; 779 /* Conforming batches */ 780 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 781 /* Remainder */ 782 PetscInt Nr, offset; 783 784 ierr = PetscFEGetQuadrature(feBd[f], &q);CHKERRQ(ierr); 785 ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr); 786 ierr = PetscFEGetTileSizes(feBd[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 787 ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 788 blockSize = Nb*numQuadPoints; 789 batchSize = numBlocks * blockSize; 790 ierr = PetscFESetTileSizes(feBd[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 791 numChunks = numFaces / (numBatches*batchSize); 792 Ne = numChunks*numBatches*batchSize; 793 Nr = numFaces % (numBatches*batchSize); 794 offset = numFaces - Nr; 795 geom.v0 = v0; 796 geom.n = n; 797 geom.J = J; 798 geom.invJ = invJ; 799 geom.detJ = detJ; 800 ierr = PetscFEIntegrateBdResidual(feBd[f], Ne, Nf, feBd, f, geom, u, 0, NULL, NULL, f0, f1, elemVec);CHKERRQ(ierr); 801 geom.v0 = &v0[offset*dim]; 802 geom.n = &n[offset*dim]; 803 geom.J = &J[offset*dim*dim]; 804 geom.invJ = &invJ[offset*dim*dim]; 805 geom.detJ = &detJ[offset]; 806 ierr = PetscFEIntegrateBdResidual(feBd[f], Nr, Nf, feBd, f, geom, &u[offset*cellDof], 0, NULL, NULL, f0, f1, &elemVec[offset*cellDof]);CHKERRQ(ierr); 807 } 808 for (p = 0, f = 0; p < numPoints; ++p) { 809 const PetscInt point = points[p]; 810 811 ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr); 812 if (dep != dim-1) continue; 813 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "BdResidual", cellDof, &elemVec[f*cellDof]);CHKERRQ(ierr);} 814 ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*cellDof], ADD_VALUES);CHKERRQ(ierr); 815 ++f; 816 } 817 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 818 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 819 ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr); 820 } 821 } 822 if (mesh->printFEM) {ierr = DMPrintLocalVec(dm, name, mesh->printTol, F);CHKERRQ(ierr);} 823 ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr); 824 PetscFunctionReturn(0); 825 } 826 827 #undef __FUNCT__ 828 #define __FUNCT__ "DMPlexComputeIFunctionFEM" 829 /*@ 830 DMPlexComputeIFunctionFEM - Form the local implicit function F from the local input X, X_t using pointwise functions specified by the user 831 832 Input Parameters: 833 + dm - The mesh 834 . time - The current time 835 . X - Local input vector 836 . X_t - Time derivative of the local input vector 837 - user - The user context 838 839 Output Parameter: 840 . F - Local output vector 841 842 Note: 843 The first member of the user context must be an FEMContext. 844 845 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 846 like a GPU, or vectorize on a multicore machine. 847 848 Level: developer 849 850 .seealso: DMPlexComputeResidualFEM() 851 @*/ 852 PetscErrorCode DMPlexComputeIFunctionFEM(DM dm, PetscReal time, Vec X, Vec X_t, Vec F, void *user) 853 { 854 DM_Plex *mesh = (DM_Plex *) dm->data; 855 PetscFEM *fem = (PetscFEM *) user; 856 PetscFE *fe = fem->fe; 857 PetscFE *feAux = fem->feAux; 858 PetscFE *feBd = fem->feBd; 859 const char *name = "Residual"; 860 DM dmAux; 861 Vec A; 862 PetscQuadrature q; 863 PetscCellGeometry geom; 864 PetscSection section, sectionAux; 865 PetscReal *v0, *J, *invJ, *detJ; 866 PetscScalar *elemVec, *u, *u_t, *a = NULL; 867 PetscInt dim, Nf, NfAux = 0, f, numCells, cStart, cEnd, c; 868 PetscInt cellDof = 0, numComponents = 0; 869 PetscInt cellDofAux = 0, numComponentsAux = 0; 870 PetscErrorCode ierr; 871 872 PetscFunctionBegin; 873 ierr = PetscLogEventBegin(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr); 874 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 875 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 876 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 877 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 878 numCells = cEnd - cStart; 879 for (f = 0; f < Nf; ++f) { 880 PetscInt Nb, Nc; 881 882 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 883 ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr); 884 cellDof += Nb*Nc; 885 numComponents += Nc; 886 } 887 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 888 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 889 if (dmAux) { 890 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 891 ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr); 892 } 893 for (f = 0; f < NfAux; ++f) { 894 PetscInt Nb, Nc; 895 896 ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr); 897 ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr); 898 cellDofAux += Nb*Nc; 899 numComponentsAux += Nc; 900 } 901 ierr = DMPlexInsertBoundaryValuesFEM(dm, X);CHKERRQ(ierr); 902 ierr = VecSet(F, 0.0);CHKERRQ(ierr); 903 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); 904 if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);} 905 for (c = cStart; c < cEnd; ++c) { 906 PetscScalar *x = NULL, *x_t = NULL; 907 PetscInt i; 908 909 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 910 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 911 ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 912 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 913 ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 914 ierr = DMPlexVecGetClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr); 915 for (i = 0; i < cellDof; ++i) u_t[c*cellDof+i] = x_t[i]; 916 ierr = DMPlexVecRestoreClosure(dm, section, X_t, c, NULL, &x_t);CHKERRQ(ierr); 917 if (dmAux) { 918 PetscScalar *x_a = NULL; 919 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x_a);CHKERRQ(ierr); 920 for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x_a[i]; 921 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x_a);CHKERRQ(ierr); 922 } 923 } 924 for (f = 0; f < Nf; ++f) { 925 void (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0IFuncs[f]; 926 void (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1IFuncs[f]; 927 PetscInt numQuadPoints, Nb; 928 /* Conforming batches */ 929 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 930 /* Remainder */ 931 PetscInt Nr, offset; 932 933 ierr = PetscFEGetQuadrature(fe[f], &q);CHKERRQ(ierr); 934 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 935 ierr = PetscFEGetTileSizes(fe[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 936 ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 937 blockSize = Nb*numQuadPoints; 938 batchSize = numBlocks * blockSize; 939 ierr = PetscFESetTileSizes(fe[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 940 numChunks = numCells / (numBatches*batchSize); 941 Ne = numChunks*numBatches*batchSize; 942 Nr = numCells % (numBatches*batchSize); 943 offset = numCells - Nr; 944 geom.v0 = v0; 945 geom.J = J; 946 geom.invJ = invJ; 947 geom.detJ = detJ; 948 ierr = PetscFEIntegrateIFunction(fe[f], Ne, Nf, fe, f, geom, u, u_t, NfAux, feAux, a, f0, f1, elemVec);CHKERRQ(ierr); 949 geom.v0 = &v0[offset*dim]; 950 geom.J = &J[offset*dim*dim]; 951 geom.invJ = &invJ[offset*dim*dim]; 952 geom.detJ = &detJ[offset]; 953 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); 954 } 955 for (c = cStart; c < cEnd; ++c) { 956 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, name, cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);} 957 ierr = DMPlexVecSetClosure(dm, section, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr); 958 } 959 ierr = PetscFree7(u,u_t,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr); 960 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 961 if (feBd) { 962 DMLabel label, depth; 963 IS pointIS; 964 const PetscInt *points; 965 PetscInt dep, numPoints, p, numFaces; 966 PetscReal *n; 967 968 ierr = DMPlexGetLabel(dm, "boundary", &label);CHKERRQ(ierr); 969 ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr); 970 ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr); 971 ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr); 972 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 973 for (f = 0, cellDof = 0, numComponents = 0; f < Nf; ++f) { 974 PetscInt Nb, Nc; 975 976 ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr); 977 ierr = PetscFEGetNumComponents(feBd[f], &Nc);CHKERRQ(ierr); 978 cellDof += Nb*Nc; 979 numComponents += Nc; 980 } 981 for (p = 0, numFaces = 0; p < numPoints; ++p) { 982 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 983 if (dep == dim-1) ++numFaces; 984 } 985 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); 986 ierr = PetscMalloc1(numFaces*cellDof,&u_t);CHKERRQ(ierr); 987 for (p = 0, f = 0; p < numPoints; ++p) { 988 const PetscInt point = points[p]; 989 PetscScalar *x = NULL; 990 PetscInt i; 991 992 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 993 if (dep != dim-1) continue; 994 ierr = DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr); 995 ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]); 996 if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point); 997 ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 998 for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i]; 999 ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 1000 ierr = DMPlexVecGetClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr); 1001 for (i = 0; i < cellDof; ++i) u_t[f*cellDof+i] = x[i]; 1002 ierr = DMPlexVecRestoreClosure(dm, section, X_t, point, NULL, &x);CHKERRQ(ierr); 1003 ++f; 1004 } 1005 for (f = 0; f < Nf; ++f) { 1006 void (*f0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdIFuncs[f]; 1007 void (*f1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdIFuncs[f]; 1008 PetscInt numQuadPoints, Nb; 1009 /* Conforming batches */ 1010 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 1011 /* Remainder */ 1012 PetscInt Nr, offset; 1013 1014 ierr = PetscFEGetQuadrature(feBd[f], &q);CHKERRQ(ierr); 1015 ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr); 1016 ierr = PetscFEGetTileSizes(feBd[f], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1017 ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 1018 blockSize = Nb*numQuadPoints; 1019 batchSize = numBlocks * blockSize; 1020 ierr = PetscFESetTileSizes(feBd[f], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1021 numChunks = numFaces / (numBatches*batchSize); 1022 Ne = numChunks*numBatches*batchSize; 1023 Nr = numFaces % (numBatches*batchSize); 1024 offset = numFaces - Nr; 1025 geom.v0 = v0; 1026 geom.n = n; 1027 geom.J = J; 1028 geom.invJ = invJ; 1029 geom.detJ = detJ; 1030 ierr = PetscFEIntegrateBdIFunction(feBd[f], Ne, Nf, feBd, f, geom, u, u_t, 0, NULL, NULL, f0, f1, elemVec);CHKERRQ(ierr); 1031 geom.v0 = &v0[offset*dim]; 1032 geom.n = &n[offset*dim]; 1033 geom.J = &J[offset*dim*dim]; 1034 geom.invJ = &invJ[offset*dim*dim]; 1035 geom.detJ = &detJ[offset]; 1036 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); 1037 } 1038 for (p = 0, f = 0; p < numPoints; ++p) { 1039 const PetscInt point = points[p]; 1040 1041 ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr); 1042 if (dep != dim-1) continue; 1043 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(point, "BdResidual", cellDof, &elemVec[f*cellDof]);CHKERRQ(ierr);} 1044 ierr = DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[f*cellDof], ADD_VALUES);CHKERRQ(ierr); 1045 ++f; 1046 } 1047 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 1048 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 1049 ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemVec);CHKERRQ(ierr); 1050 ierr = PetscFree(u_t);CHKERRQ(ierr); 1051 } 1052 if (mesh->printFEM) {ierr = DMPrintLocalVec(dm, name, mesh->printTol, F);CHKERRQ(ierr);} 1053 ierr = PetscLogEventEnd(DMPLEX_ResidualFEM,dm,0,0,0);CHKERRQ(ierr); 1054 PetscFunctionReturn(0); 1055 } 1056 1057 #undef __FUNCT__ 1058 #define __FUNCT__ "DMPlexComputeJacobianActionFEM" 1059 /*@C 1060 DMPlexComputeJacobianActionFEM - Form the local action of Jacobian J(u) on the local input X using pointwise functions specified by the user 1061 1062 Input Parameters: 1063 + dm - The mesh 1064 . J - The Jacobian shell matrix 1065 . X - Local input vector 1066 - user - The user context 1067 1068 Output Parameter: 1069 . F - Local output vector 1070 1071 Note: 1072 The first member of the user context must be an FEMContext. 1073 1074 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 1075 like a GPU, or vectorize on a multicore machine. 1076 1077 Level: developer 1078 1079 .seealso: DMPlexComputeResidualFEM() 1080 @*/ 1081 PetscErrorCode DMPlexComputeJacobianActionFEM(DM dm, Mat Jac, Vec X, Vec F, void *user) 1082 { 1083 DM_Plex *mesh = (DM_Plex *) dm->data; 1084 PetscFEM *fem = (PetscFEM *) user; 1085 PetscFE *fe = fem->fe; 1086 PetscQuadrature quad; 1087 PetscCellGeometry geom; 1088 PetscSection section; 1089 JacActionCtx *jctx; 1090 PetscReal *v0, *J, *invJ, *detJ; 1091 PetscScalar *elemVec, *u, *a; 1092 PetscInt dim, numFields, field, numCells, cStart, cEnd, c; 1093 PetscInt cellDof = 0; 1094 PetscErrorCode ierr; 1095 1096 PetscFunctionBegin; 1097 /* ierr = PetscLogEventBegin(DMPLEX_JacobianActionFEM,dm,0,0,0);CHKERRQ(ierr); */ 1098 ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr); 1099 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1100 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1101 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 1102 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1103 numCells = cEnd - cStart; 1104 for (field = 0; field < numFields; ++field) { 1105 PetscInt Nb, Nc; 1106 1107 ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr); 1108 ierr = PetscFEGetNumComponents(fe[field], &Nc);CHKERRQ(ierr); 1109 cellDof += Nb*Nc; 1110 } 1111 ierr = VecSet(F, 0.0);CHKERRQ(ierr); 1112 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); 1113 for (c = cStart; c < cEnd; ++c) { 1114 PetscScalar *x = NULL; 1115 PetscInt i; 1116 1117 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 1118 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 1119 ierr = DMPlexVecGetClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr); 1120 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 1121 ierr = DMPlexVecRestoreClosure(dm, NULL, jctx->u, c, NULL, &x);CHKERRQ(ierr); 1122 ierr = DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 1123 for (i = 0; i < cellDof; ++i) a[c*cellDof+i] = x[i]; 1124 ierr = DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);CHKERRQ(ierr); 1125 } 1126 for (field = 0; field < numFields; ++field) { 1127 PetscInt numQuadPoints, Nb; 1128 /* Conforming batches */ 1129 PetscInt numBlocks = 1; 1130 PetscInt numBatches = 1; 1131 PetscInt numChunks, Ne, blockSize, batchSize; 1132 /* Remainder */ 1133 PetscInt Nr, offset; 1134 1135 ierr = PetscFEGetQuadrature(fe[field], &quad);CHKERRQ(ierr); 1136 ierr = PetscFEGetDimension(fe[field], &Nb);CHKERRQ(ierr); 1137 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 1138 blockSize = Nb*numQuadPoints; 1139 batchSize = numBlocks * blockSize; 1140 numChunks = numCells / (numBatches*batchSize); 1141 Ne = numChunks*numBatches*batchSize; 1142 Nr = numCells % (numBatches*batchSize); 1143 offset = numCells - Nr; 1144 geom.v0 = v0; 1145 geom.J = J; 1146 geom.invJ = invJ; 1147 geom.detJ = detJ; 1148 ierr = PetscFEIntegrateJacobianAction(fe[field], Ne, numFields, fe, field, geom, u, a, fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, elemVec);CHKERRQ(ierr); 1149 geom.v0 = &v0[offset*dim]; 1150 geom.J = &J[offset*dim*dim]; 1151 geom.invJ = &invJ[offset*dim*dim]; 1152 geom.detJ = &detJ[offset]; 1153 ierr = PetscFEIntegrateJacobianAction(fe[field], Nr, numFields, fe, field, geom, &u[offset*cellDof], &a[offset*cellDof], 1154 fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, &elemVec[offset*cellDof]);CHKERRQ(ierr); 1155 } 1156 for (c = cStart; c < cEnd; ++c) { 1157 if (mesh->printFEM > 1) {ierr = DMPrintCellVector(c, "Jacobian Action", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);} 1158 ierr = DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr); 1159 } 1160 ierr = PetscFree7(u,a,v0,J,invJ,detJ,elemVec);CHKERRQ(ierr); 1161 if (mesh->printFEM) { 1162 PetscMPIInt rank, numProcs; 1163 PetscInt p; 1164 1165 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr); 1166 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);CHKERRQ(ierr); 1167 ierr = PetscPrintf(PetscObjectComm((PetscObject)dm), "Jacobian Action:\n");CHKERRQ(ierr); 1168 for (p = 0; p < numProcs; ++p) { 1169 if (p == rank) {ierr = VecView(F, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);} 1170 ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr); 1171 } 1172 } 1173 /* ierr = PetscLogEventEnd(DMPLEX_JacobianActionFEM,dm,0,0,0);CHKERRQ(ierr); */ 1174 PetscFunctionReturn(0); 1175 } 1176 1177 #undef __FUNCT__ 1178 #define __FUNCT__ "DMPlexComputeJacobianFEM" 1179 /*@ 1180 DMPlexComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user. 1181 1182 Input Parameters: 1183 + dm - The mesh 1184 . X - Local input vector 1185 - user - The user context 1186 1187 Output Parameter: 1188 . Jac - Jacobian matrix 1189 1190 Note: 1191 The first member of the user context must be an FEMContext. 1192 1193 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 1194 like a GPU, or vectorize on a multicore machine. 1195 1196 Level: developer 1197 1198 .seealso: FormFunctionLocal() 1199 @*/ 1200 PetscErrorCode DMPlexComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP,void *user) 1201 { 1202 DM_Plex *mesh = (DM_Plex *) dm->data; 1203 PetscFEM *fem = (PetscFEM *) user; 1204 PetscFE *fe = fem->fe; 1205 PetscFE *feAux = fem->feAux; 1206 PetscFE *feBd = fem->feBd; 1207 const char *name = "Jacobian"; 1208 DM dmAux; 1209 Vec A; 1210 PetscQuadrature quad; 1211 PetscCellGeometry geom; 1212 PetscSection section, globalSection, sectionAux; 1213 PetscReal *v0, *J, *invJ, *detJ; 1214 PetscScalar *elemMat, *u, *a; 1215 PetscInt dim, Nf, NfAux = 0, f, fieldI, fieldJ, numCells, cStart, cEnd, c; 1216 PetscInt cellDof = 0, numComponents = 0; 1217 PetscInt cellDofAux = 0, numComponentsAux = 0; 1218 PetscBool isShell; 1219 PetscErrorCode ierr; 1220 1221 PetscFunctionBegin; 1222 ierr = PetscLogEventBegin(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr); 1223 ierr = DMPlexGetDimension(dm, &dim);CHKERRQ(ierr); 1224 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 1225 ierr = DMGetDefaultGlobalSection(dm, &globalSection);CHKERRQ(ierr); 1226 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 1227 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1228 numCells = cEnd - cStart; 1229 for (f = 0; f < Nf; ++f) { 1230 PetscInt Nb, Nc; 1231 1232 ierr = PetscFEGetDimension(fe[f], &Nb);CHKERRQ(ierr); 1233 ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr); 1234 cellDof += Nb*Nc; 1235 numComponents += Nc; 1236 } 1237 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 1238 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1239 if (dmAux) { 1240 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1241 ierr = PetscSectionGetNumFields(sectionAux, &NfAux);CHKERRQ(ierr); 1242 } 1243 for (f = 0; f < NfAux; ++f) { 1244 PetscInt Nb, Nc; 1245 1246 ierr = PetscFEGetDimension(feAux[f], &Nb);CHKERRQ(ierr); 1247 ierr = PetscFEGetNumComponents(feAux[f], &Nc);CHKERRQ(ierr); 1248 cellDofAux += Nb*Nc; 1249 numComponentsAux += Nc; 1250 } 1251 ierr = DMPlexInsertBoundaryValuesFEM(dm, X);CHKERRQ(ierr); 1252 ierr = MatZeroEntries(JacP);CHKERRQ(ierr); 1253 ierr = PetscMalloc6(numCells*cellDof,&u,numCells*dim,&v0,numCells*dim*dim,&J,numCells*dim*dim,&invJ,numCells,&detJ,numCells*cellDof*cellDof,&elemMat);CHKERRQ(ierr); 1254 if (dmAux) {ierr = PetscMalloc1(numCells*cellDofAux, &a);CHKERRQ(ierr);} 1255 for (c = cStart; c < cEnd; ++c) { 1256 PetscScalar *x = NULL; 1257 PetscInt i; 1258 1259 ierr = DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);CHKERRQ(ierr); 1260 if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c); 1261 ierr = DMPlexVecGetClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 1262 for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i]; 1263 ierr = DMPlexVecRestoreClosure(dm, section, X, c, NULL, &x);CHKERRQ(ierr); 1264 if (dmAux) { 1265 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1266 for (i = 0; i < cellDofAux; ++i) a[c*cellDofAux+i] = x[i]; 1267 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1268 } 1269 } 1270 ierr = PetscMemzero(elemMat, numCells*cellDof*cellDof * sizeof(PetscScalar));CHKERRQ(ierr); 1271 for (fieldI = 0; fieldI < Nf; ++fieldI) { 1272 PetscInt numQuadPoints, Nb; 1273 /* Conforming batches */ 1274 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 1275 /* Remainder */ 1276 PetscInt Nr, offset; 1277 1278 ierr = PetscFEGetQuadrature(fe[fieldI], &quad);CHKERRQ(ierr); 1279 ierr = PetscFEGetDimension(fe[fieldI], &Nb);CHKERRQ(ierr); 1280 ierr = PetscFEGetTileSizes(fe[fieldI], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1281 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 1282 blockSize = Nb*numQuadPoints; 1283 batchSize = numBlocks * blockSize; 1284 ierr = PetscFESetTileSizes(fe[fieldI], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1285 numChunks = numCells / (numBatches*batchSize); 1286 Ne = numChunks*numBatches*batchSize; 1287 Nr = numCells % (numBatches*batchSize); 1288 offset = numCells - Nr; 1289 for (fieldJ = 0; fieldJ < Nf; ++fieldJ) { 1290 void (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g0Funcs[fieldI*Nf+fieldJ]; 1291 void (*g1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g1Funcs[fieldI*Nf+fieldJ]; 1292 void (*g2)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g2Funcs[fieldI*Nf+fieldJ]; 1293 void (*g3)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g3Funcs[fieldI*Nf+fieldJ]; 1294 1295 geom.v0 = v0; 1296 geom.J = J; 1297 geom.invJ = invJ; 1298 geom.detJ = detJ; 1299 ierr = PetscFEIntegrateJacobian(fe[fieldI], Ne, Nf, fe, fieldI, fieldJ, geom, u, NfAux, feAux, a, g0, g1, g2, g3, elemMat);CHKERRQ(ierr); 1300 geom.v0 = &v0[offset*dim]; 1301 geom.J = &J[offset*dim*dim]; 1302 geom.invJ = &invJ[offset*dim*dim]; 1303 geom.detJ = &detJ[offset]; 1304 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); 1305 } 1306 } 1307 for (c = cStart; c < cEnd; ++c) { 1308 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(c, name, cellDof, cellDof, &elemMat[c*cellDof*cellDof]);CHKERRQ(ierr);} 1309 ierr = DMPlexMatSetClosure(dm, section, globalSection, JacP, c, &elemMat[c*cellDof*cellDof], ADD_VALUES);CHKERRQ(ierr); 1310 } 1311 ierr = PetscFree6(u,v0,J,invJ,detJ,elemMat);CHKERRQ(ierr); 1312 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1313 if (feBd) { 1314 DMLabel depth; 1315 PetscInt numBd, bd; 1316 1317 for (f = 0, cellDof = 0, numComponents = 0; f < Nf; ++f) { 1318 PetscInt Nb, Nc; 1319 1320 ierr = PetscFEGetDimension(feBd[f], &Nb);CHKERRQ(ierr); 1321 ierr = PetscFEGetNumComponents(feBd[f], &Nc);CHKERRQ(ierr); 1322 cellDof += Nb*Nc; 1323 numComponents += Nc; 1324 } 1325 ierr = DMPlexGetDepthLabel(dm, &depth);CHKERRQ(ierr); 1326 ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 1327 for (bd = 0; bd < numBd; ++bd) { 1328 const char *bdLabel; 1329 DMLabel label; 1330 IS pointIS; 1331 const PetscInt *points; 1332 const PetscInt *values; 1333 PetscReal *n; 1334 PetscInt field, numValues, numPoints, p, dep, numFaces; 1335 PetscBool isEssential; 1336 1337 ierr = DMPlexGetBoundary(dm, bd, &isEssential, NULL, &bdLabel, &field, NULL, &numValues, &values, NULL);CHKERRQ(ierr); 1338 if (numValues != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Bug me and I will fix this"); 1339 ierr = DMPlexGetLabel(dm, bdLabel, &label);CHKERRQ(ierr); 1340 ierr = DMLabelGetStratumSize(label, 1, &numPoints);CHKERRQ(ierr); 1341 ierr = DMLabelGetStratumIS(label, 1, &pointIS);CHKERRQ(ierr); 1342 ierr = ISGetIndices(pointIS, &points);CHKERRQ(ierr); 1343 for (p = 0, numFaces = 0; p < numPoints; ++p) { 1344 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 1345 if (dep == dim-1) ++numFaces; 1346 } 1347 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); 1348 for (p = 0, f = 0; p < numPoints; ++p) { 1349 const PetscInt point = points[p]; 1350 PetscScalar *x = NULL; 1351 PetscInt i; 1352 1353 ierr = DMLabelGetValue(depth, points[p], &dep);CHKERRQ(ierr); 1354 if (dep != dim-1) continue; 1355 ierr = DMPlexComputeCellGeometry(dm, point, &v0[f*dim], &J[f*dim*dim], &invJ[f*dim*dim], &detJ[f]);CHKERRQ(ierr); 1356 ierr = DMPlexComputeCellGeometryFVM(dm, point, NULL, NULL, &n[f*dim]); 1357 if (detJ[f] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[f], point); 1358 ierr = DMPlexVecGetClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 1359 for (i = 0; i < cellDof; ++i) u[f*cellDof+i] = x[i]; 1360 ierr = DMPlexVecRestoreClosure(dm, section, X, point, NULL, &x);CHKERRQ(ierr); 1361 ++f; 1362 } 1363 ierr = PetscMemzero(elemMat, numFaces*cellDof*cellDof * sizeof(PetscScalar));CHKERRQ(ierr); 1364 for (fieldI = 0; fieldI < Nf; ++fieldI) { 1365 PetscInt numQuadPoints, Nb; 1366 /* Conforming batches */ 1367 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 1368 /* Remainder */ 1369 PetscInt Nr, offset; 1370 1371 ierr = PetscFEGetQuadrature(feBd[fieldI], &quad);CHKERRQ(ierr); 1372 ierr = PetscFEGetDimension(feBd[fieldI], &Nb);CHKERRQ(ierr); 1373 ierr = PetscFEGetTileSizes(feBd[fieldI], NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1374 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 1375 blockSize = Nb*numQuadPoints; 1376 batchSize = numBlocks * blockSize; 1377 ierr = PetscFESetTileSizes(feBd[fieldI], blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1378 numChunks = numFaces / (numBatches*batchSize); 1379 Ne = numChunks*numBatches*batchSize; 1380 Nr = numFaces % (numBatches*batchSize); 1381 offset = numFaces - Nr; 1382 for (fieldJ = 0; fieldJ < Nf; ++fieldJ) { 1383 void (*g0)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->g0BdFuncs[fieldI*Nf+fieldJ]; 1384 void (*g1)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->g1BdFuncs[fieldI*Nf+fieldJ]; 1385 void (*g2)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->g2BdFuncs[fieldI*Nf+fieldJ]; 1386 void (*g3)(const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->g3BdFuncs[fieldI*Nf+fieldJ]; 1387 1388 geom.v0 = v0; 1389 geom.n = n; 1390 geom.J = J; 1391 geom.invJ = invJ; 1392 geom.detJ = detJ; 1393 ierr = PetscFEIntegrateBdJacobian(feBd[fieldI], Ne, Nf, feBd, fieldI, fieldJ, geom, u, 0, NULL, NULL, g0, g1, g2, g3, elemMat);CHKERRQ(ierr); 1394 geom.v0 = &v0[offset*dim]; 1395 geom.n = &n[offset*dim]; 1396 geom.J = &J[offset*dim*dim]; 1397 geom.invJ = &invJ[offset*dim*dim]; 1398 geom.detJ = &detJ[offset]; 1399 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); 1400 } 1401 } 1402 for (p = 0, f = 0; p < numPoints; ++p) { 1403 const PetscInt point = points[p]; 1404 1405 ierr = DMLabelGetValue(depth, point, &dep);CHKERRQ(ierr); 1406 if (dep != dim-1) continue; 1407 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(point, "BdJacobian", cellDof, cellDof, &elemMat[f*cellDof*cellDof]);CHKERRQ(ierr);} 1408 ierr = DMPlexMatSetClosure(dm, section, globalSection, JacP, point, &elemMat[f*cellDof*cellDof], ADD_VALUES);CHKERRQ(ierr); 1409 ++f; 1410 } 1411 ierr = ISRestoreIndices(pointIS, &points);CHKERRQ(ierr); 1412 ierr = ISDestroy(&pointIS);CHKERRQ(ierr); 1413 ierr = PetscFree7(u,v0,n,J,invJ,detJ,elemMat);CHKERRQ(ierr); 1414 } 1415 } 1416 ierr = MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1417 ierr = MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1418 if (mesh->printFEM) { 1419 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1420 ierr = MatChop(JacP, 1.0e-10);CHKERRQ(ierr); 1421 ierr = MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1422 } 1423 ierr = PetscLogEventEnd(DMPLEX_JacobianFEM,dm,0,0,0);CHKERRQ(ierr); 1424 ierr = PetscObjectTypeCompare((PetscObject) Jac, MATSHELL, &isShell);CHKERRQ(ierr); 1425 if (isShell) { 1426 JacActionCtx *jctx; 1427 1428 ierr = MatShellGetContext(Jac, &jctx);CHKERRQ(ierr); 1429 ierr = VecCopy(X, jctx->u);CHKERRQ(ierr); 1430 } 1431 PetscFunctionReturn(0); 1432 } 1433 1434 #undef __FUNCT__ 1435 #define __FUNCT__ "DMPlexComputeInterpolatorFEM" 1436 /*@ 1437 DMPlexComputeInterpolatorFEM - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1438 1439 Input Parameters: 1440 + dmf - The fine mesh 1441 . dmc - The coarse mesh 1442 - user - The user context 1443 1444 Output Parameter: 1445 . In - The interpolation matrix 1446 1447 Note: 1448 The first member of the user context must be an FEMContext. 1449 1450 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 1451 like a GPU, or vectorize on a multicore machine. 1452 1453 Level: developer 1454 1455 .seealso: DMPlexComputeJacobianFEM() 1456 @*/ 1457 PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user) 1458 { 1459 DM_Plex *mesh = (DM_Plex *) dmc->data; 1460 PetscFEM *fem = (PetscFEM *) user; 1461 PetscFE *fe = fem->fe; 1462 const char *name = "Interpolator"; 1463 PetscFE *feRef; 1464 PetscSection fsection, fglobalSection; 1465 PetscSection csection, cglobalSection; 1466 PetscScalar *elemMat; 1467 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, c; 1468 PetscInt rCellDof = 0, cCellDof = 0; 1469 PetscErrorCode ierr; 1470 1471 PetscFunctionBegin; 1472 #if 0 1473 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1474 #endif 1475 ierr = DMPlexGetDimension(dmf, &dim);CHKERRQ(ierr); 1476 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1477 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1478 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1479 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1480 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1481 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1482 ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr); 1483 for (f = 0; f < Nf; ++f) { 1484 PetscInt rNb, cNb, Nc; 1485 1486 ierr = PetscFERefine(fe[f], &feRef[f]);CHKERRQ(ierr); 1487 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1488 ierr = PetscFEGetDimension(fe[f], &cNb);CHKERRQ(ierr); 1489 ierr = PetscFEGetNumComponents(fe[f], &Nc);CHKERRQ(ierr); 1490 rCellDof += rNb*Nc; 1491 cCellDof += cNb*Nc; 1492 } 1493 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1494 ierr = PetscMalloc1(rCellDof*cCellDof,&elemMat);CHKERRQ(ierr); 1495 ierr = PetscMemzero(elemMat, rCellDof*cCellDof * sizeof(PetscScalar));CHKERRQ(ierr); 1496 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1497 PetscDualSpace Qref; 1498 PetscQuadrature f; 1499 const PetscReal *qpoints, *qweights; 1500 PetscReal *points; 1501 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1502 1503 /* Compose points from all dual basis functionals */ 1504 ierr = PetscFEGetNumComponents(fe[fieldI], &Nc);CHKERRQ(ierr); 1505 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1506 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1507 for (i = 0; i < fpdim; ++i) { 1508 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1509 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1510 npoints += Np; 1511 } 1512 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1513 for (i = 0, k = 0; i < fpdim; ++i) { 1514 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1515 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1516 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1517 } 1518 1519 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1520 PetscReal *B; 1521 PetscInt NcJ, cpdim, j; 1522 1523 /* Evaluate basis at points */ 1524 ierr = PetscFEGetNumComponents(fe[fieldJ], &NcJ);CHKERRQ(ierr); 1525 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); 1526 ierr = PetscFEGetDimension(fe[fieldJ], &cpdim);CHKERRQ(ierr); 1527 /* For now, fields only interpolate themselves */ 1528 if (fieldI == fieldJ) { 1529 ierr = PetscFEGetTabulation(fe[fieldJ], npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1530 for (i = 0, k = 0; i < fpdim; ++i) { 1531 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1532 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1533 for (p = 0; p < Np; ++p, ++k) { 1534 for (j = 0; j < cpdim; ++j) { 1535 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]; 1536 } 1537 } 1538 } 1539 ierr = PetscFERestoreTabulation(fe[fieldJ], npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1540 } 1541 offsetJ += cpdim*NcJ; 1542 } 1543 offsetI += fpdim*Nc; 1544 ierr = PetscFree(points);CHKERRQ(ierr); 1545 } 1546 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rCellDof, cCellDof, elemMat);CHKERRQ(ierr);} 1547 for (c = cStart; c < cEnd; ++c) { 1548 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1549 } 1550 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1551 ierr = PetscFree(feRef);CHKERRQ(ierr); 1552 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1553 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1554 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1555 if (mesh->printFEM) { 1556 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1557 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1558 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1559 } 1560 #if 0 1561 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1562 #endif 1563 PetscFunctionReturn(0); 1564 } 1565 1566 #undef __FUNCT__ 1567 #define __FUNCT__ "DMPlexAddBoundary" 1568 /* The ids can be overridden by the command line option -bc_<boundary name> */ 1569 PetscErrorCode DMPlexAddBoundary(DM dm, PetscBool isEssential, const char name[], const char labelname[], PetscInt field, void (*bcFunc)(), PetscInt numids, const PetscInt *ids, void *ctx) 1570 { 1571 DM_Plex *mesh = (DM_Plex *) dm->data; 1572 DMBoundary b; 1573 PetscErrorCode ierr; 1574 1575 PetscFunctionBegin; 1576 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1577 ierr = PetscNew(&b);CHKERRQ(ierr); 1578 ierr = PetscStrallocpy(name, (char **) &b->name);CHKERRQ(ierr); 1579 ierr = PetscStrallocpy(labelname, (char **) &b->labelname);CHKERRQ(ierr); 1580 ierr = PetscMalloc1(numids, &b->ids);CHKERRQ(ierr); 1581 ierr = PetscMemcpy(b->ids, ids, numids*sizeof(PetscInt));CHKERRQ(ierr); 1582 if (b->labelname) { 1583 ierr = DMPlexGetLabel(dm, b->labelname, &b->label);CHKERRQ(ierr); 1584 if (!b->label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Label %s does not exist in this DM", b->labelname); 1585 } 1586 b->essential = isEssential; 1587 b->field = field; 1588 b->func = bcFunc; 1589 b->numids = numids; 1590 b->ctx = ctx; 1591 b->next = mesh->boundary; 1592 mesh->boundary = b; 1593 PetscFunctionReturn(0); 1594 } 1595 1596 #undef __FUNCT__ 1597 #define __FUNCT__ "DMPlexGetNumBoundary" 1598 PetscErrorCode DMPlexGetNumBoundary(DM dm, PetscInt *numBd) 1599 { 1600 DM_Plex *mesh = (DM_Plex *) dm->data; 1601 DMBoundary b = mesh->boundary; 1602 1603 PetscFunctionBegin; 1604 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1605 PetscValidPointer(numBd, 2); 1606 *numBd = 0; 1607 while (b) {++(*numBd); b = b->next;} 1608 PetscFunctionReturn(0); 1609 } 1610 1611 #undef __FUNCT__ 1612 #define __FUNCT__ "DMPlexGetBoundary" 1613 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) 1614 { 1615 DM_Plex *mesh = (DM_Plex *) dm->data; 1616 DMBoundary b = mesh->boundary; 1617 PetscInt n = 0; 1618 1619 PetscFunctionBegin; 1620 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1621 while (b) { 1622 if (n == bd) break; 1623 b = b->next; 1624 ++n; 1625 } 1626 if (n != bd) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n); 1627 if (isEssential) { 1628 PetscValidPointer(isEssential, 3); 1629 *isEssential = b->essential; 1630 } 1631 if (name) { 1632 PetscValidPointer(name, 4); 1633 *name = b->name; 1634 } 1635 if (labelname) { 1636 PetscValidPointer(labelname, 5); 1637 *labelname = b->labelname; 1638 } 1639 if (field) { 1640 PetscValidPointer(field, 6); 1641 *field = b->field; 1642 } 1643 if (func) { 1644 PetscValidPointer(func, 7); 1645 *func = b->func; 1646 } 1647 if (numids) { 1648 PetscValidPointer(numids, 8); 1649 *numids = b->numids; 1650 } 1651 if (ids) { 1652 PetscValidPointer(ids, 9); 1653 *ids = b->ids; 1654 } 1655 if (ctx) { 1656 PetscValidPointer(ctx, 10); 1657 *ctx = b->ctx; 1658 } 1659 PetscFunctionReturn(0); 1660 } 1661 1662 #undef __FUNCT__ 1663 #define __FUNCT__ "DMPlexIsBoundaryPoint" 1664 PetscErrorCode DMPlexIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd) 1665 { 1666 DM_Plex *mesh = (DM_Plex *) dm->data; 1667 DMBoundary b = mesh->boundary; 1668 PetscErrorCode ierr; 1669 1670 PetscFunctionBegin; 1671 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 1672 PetscValidPointer(isBd, 3); 1673 *isBd = PETSC_FALSE; 1674 while (b && !(*isBd)) { 1675 if (b->label) { 1676 PetscInt i; 1677 1678 for (i = 0; i < b->numids && !(*isBd); ++i) { 1679 ierr = DMLabelStratumHasPoint(b->label, b->ids[i], point, isBd);CHKERRQ(ierr); 1680 } 1681 } 1682 b = b->next; 1683 } 1684 PetscFunctionReturn(0); 1685 } 1686