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