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