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