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