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