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