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