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