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