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