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[], void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX) 187 { 188 PetscFE fe; 189 PetscDualSpace *sp; 190 PetscSection section; 191 PetscScalar *values; 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 = PetscMalloc1(numFields,&sp);CHKERRQ(ierr); 201 for (f = 0; f < numFields; ++f) { 202 ierr = DMGetField(dm, f, (PetscObject *) &fe);CHKERRQ(ierr); 203 ierr = PetscFEGetDualSpace(fe, &sp[f]);CHKERRQ(ierr); 204 ierr = PetscFEGetNumComponents(fe, &numComp);CHKERRQ(ierr); 205 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 206 totDim += spDim*numComp; 207 } 208 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 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 PetscFECellGeom geom; 225 226 if ((point < cStart) || (point >= cEnd)) continue; 227 ierr = DMPlexComputeCellGeometryFEM(dm, point, NULL, geom.v0, geom.J, NULL, &geom.detJ);CHKERRQ(ierr); 228 for (f = 0, v = 0; f < numFields; ++f) { 229 void * const ctx = ctxs ? ctxs[f] : NULL; 230 231 ierr = DMGetField(dm, f, (PetscObject *) &fe);CHKERRQ(ierr); 232 ierr = PetscFEGetNumComponents(fe, &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 = PetscFree(sp);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 PetscInt numFields, dim, spDim, totDim = 0, numValues, cStart, cEnd, c, f, d, v, comp; 263 PetscErrorCode ierr; 264 265 PetscFunctionBegin; 266 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 267 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 268 ierr = PetscMalloc2(numFields, &sp, numFields, &numComp);CHKERRQ(ierr); 269 for (f = 0; f < numFields; ++f) { 270 PetscObject obj; 271 PetscClassId id; 272 273 ierr = DMGetField(dm, f, &obj);CHKERRQ(ierr); 274 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 275 if (id == PETSCFE_CLASSID) { 276 PetscFE fe = (PetscFE) obj; 277 278 ierr = PetscFEGetNumComponents(fe, &numComp[f]);CHKERRQ(ierr); 279 ierr = PetscFEGetDualSpace(fe, &sp[f]);CHKERRQ(ierr); 280 ierr = PetscObjectReference((PetscObject) sp[f]);CHKERRQ(ierr); 281 } else if (id == PETSCFV_CLASSID) { 282 PetscFV fv = (PetscFV) obj; 283 PetscQuadrature q; 284 285 ierr = PetscFVGetNumComponents(fv, &numComp[f]);CHKERRQ(ierr); 286 ierr = PetscFVGetQuadrature(fv, &q);CHKERRQ(ierr); 287 ierr = PetscDualSpaceCreate(PetscObjectComm((PetscObject) fv), &sp[f]);CHKERRQ(ierr); 288 ierr = PetscDualSpaceSetDM(sp[f], dm);CHKERRQ(ierr); 289 ierr = PetscDualSpaceSetType(sp[f], PETSCDUALSPACESIMPLE);CHKERRQ(ierr); 290 ierr = PetscDualSpaceSimpleSetDimension(sp[f], 1);CHKERRQ(ierr); 291 ierr = PetscDualSpaceSimpleSetFunctional(sp[f], 0, q);CHKERRQ(ierr); 292 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f); 293 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 294 totDim += spDim*numComp[f]; 295 } 296 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 297 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 298 ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 299 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 300 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 301 for (c = cStart; c < cEnd; ++c) { 302 PetscFECellGeom geom; 303 304 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, geom.v0, geom.J, NULL, &geom.detJ);CHKERRQ(ierr); 305 for (f = 0, v = 0; f < numFields; ++f) { 306 void *const ctx = ctxs ? ctxs[f] : NULL; 307 308 ierr = PetscDualSpaceGetDimension(sp[f], &spDim);CHKERRQ(ierr); 309 for (d = 0; d < spDim; ++d) { 310 if (funcs[f]) { 311 ierr = PetscDualSpaceApply(sp[f], d, &geom, numComp[f], funcs[f], ctx, &values[v]);CHKERRQ(ierr); 312 } else { 313 for (comp = 0; comp < numComp[f]; ++comp) values[v+comp] = 0.0; 314 } 315 v += numComp[f]; 316 } 317 } 318 ierr = DMPlexVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr); 319 } 320 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 321 for (f = 0; f < numFields; ++f) {ierr = PetscDualSpaceDestroy(&sp[f]);CHKERRQ(ierr);} 322 ierr = PetscFree2(sp, numComp);CHKERRQ(ierr); 323 PetscFunctionReturn(0); 324 } 325 326 #undef __FUNCT__ 327 #define __FUNCT__ "DMPlexProjectFunction" 328 /*@C 329 DMPlexProjectFunction - This projects the given function into the function space provided. 330 331 Input Parameters: 332 + dm - The DM 333 . funcs - The coordinate functions to evaluate, one per field 334 . ctxs - Optional array of contexts to pass to each coordinate function. ctxs itself may be null. 335 - mode - The insertion mode for values 336 337 Output Parameter: 338 . X - vector 339 340 Level: developer 341 342 .seealso: DMPlexComputeL2Diff() 343 @*/ 344 PetscErrorCode DMPlexProjectFunction(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X) 345 { 346 Vec localX; 347 PetscErrorCode ierr; 348 349 PetscFunctionBegin; 350 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 351 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 352 ierr = DMPlexProjectFunctionLocal(dm, funcs, ctxs, mode, localX);CHKERRQ(ierr); 353 ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr); 354 ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr); 355 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 356 PetscFunctionReturn(0); 357 } 358 359 #undef __FUNCT__ 360 #define __FUNCT__ "DMPlexProjectFieldLocal" 361 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) 362 { 363 DM dmAux; 364 PetscDS prob, probAux; 365 Vec A; 366 PetscSection section, sectionAux; 367 PetscScalar *values, *u, *u_x, *a, *a_x; 368 PetscReal *x, *v0, *J, *invJ, detJ, **basisField, **basisFieldDer, **basisFieldAux, **basisFieldDerAux; 369 PetscInt Nf, dim, spDim, totDim, numValues, cStart, cEnd, c, f, d, v, comp; 370 PetscErrorCode ierr; 371 372 PetscFunctionBegin; 373 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 374 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 375 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 376 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 377 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 378 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 379 ierr = PetscDSGetTabulation(prob, &basisField, &basisFieldDer);CHKERRQ(ierr); 380 ierr = PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);CHKERRQ(ierr); 381 ierr = PetscDSGetRefCoordArrays(prob, &x, NULL);CHKERRQ(ierr); 382 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 383 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 384 if (dmAux) { 385 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 386 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 387 ierr = PetscDSGetTabulation(prob, &basisFieldAux, &basisFieldDerAux);CHKERRQ(ierr); 388 ierr = PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);CHKERRQ(ierr); 389 } 390 ierr = DMPlexInsertBoundaryValues(dm, localU, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 391 ierr = DMPlexVecGetClosure(dm, section, localX, cStart, &numValues, NULL);CHKERRQ(ierr); 392 if (numValues != totDim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The section cell closure size %d != dual space dimension %d", numValues, totDim); 393 ierr = DMGetWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 394 ierr = PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 395 for (c = cStart; c < cEnd; ++c) { 396 PetscScalar *coefficients = NULL, *coefficientsAux = NULL; 397 398 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 399 ierr = DMPlexVecGetClosure(dm, section, localU, c, NULL, &coefficients);CHKERRQ(ierr); 400 if (dmAux) {ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);CHKERRQ(ierr);} 401 for (f = 0, v = 0; f < Nf; ++f) { 402 PetscFE fe; 403 PetscDualSpace sp; 404 PetscInt Ncf; 405 406 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 407 ierr = PetscFEGetDualSpace(fe, &sp);CHKERRQ(ierr); 408 ierr = PetscFEGetNumComponents(fe, &Ncf);CHKERRQ(ierr); 409 ierr = PetscDualSpaceGetDimension(sp, &spDim);CHKERRQ(ierr); 410 for (d = 0; d < spDim; ++d) { 411 PetscQuadrature quad; 412 const PetscReal *points, *weights; 413 PetscInt numPoints, q; 414 415 if (funcs[f]) { 416 ierr = PetscDualSpaceGetFunctional(sp, d, &quad);CHKERRQ(ierr); 417 ierr = PetscQuadratureGetData(quad, NULL, &numPoints, &points, &weights);CHKERRQ(ierr); 418 ierr = PetscFEGetTabulation(fe, numPoints, points, &basisField[f], &basisFieldDer[f], NULL);CHKERRQ(ierr); 419 for (q = 0; q < numPoints; ++q) { 420 CoordinatesRefToReal(dim, dim, v0, J, &points[q*dim], x); 421 ierr = EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, coefficients, NULL, u, u_x, NULL);CHKERRQ(ierr); 422 ierr = EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, coefficientsAux, NULL, a, a_x, NULL);CHKERRQ(ierr); 423 (*funcs[f])(u, NULL, u_x, a, NULL, a_x, x, &values[v]); 424 } 425 ierr = PetscFERestoreTabulation(fe, numPoints, points, &basisField[f], &basisFieldDer[f], NULL);CHKERRQ(ierr); 426 } else { 427 for (comp = 0; comp < Ncf; ++comp) values[v+comp] = 0.0; 428 } 429 v += Ncf; 430 } 431 } 432 ierr = DMPlexVecRestoreClosure(dm, section, localU, c, NULL, &coefficients);CHKERRQ(ierr); 433 if (dmAux) {ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &coefficientsAux);CHKERRQ(ierr);} 434 ierr = DMPlexVecSetClosure(dm, section, localX, c, values, mode);CHKERRQ(ierr); 435 } 436 ierr = DMRestoreWorkArray(dm, numValues, PETSC_SCALAR, &values);CHKERRQ(ierr); 437 ierr = PetscFree3(v0,J,invJ);CHKERRQ(ierr); 438 PetscFunctionReturn(0); 439 } 440 441 #undef __FUNCT__ 442 #define __FUNCT__ "DMPlexProjectField" 443 /*@C 444 DMPlexProjectField - This projects the given function of the fields into the function space provided. 445 446 Input Parameters: 447 + dm - The DM 448 . U - The input field vector 449 . funcs - The functions to evaluate, one per field 450 - mode - The insertion mode for values 451 452 Output Parameter: 453 . X - The output vector 454 455 Level: developer 456 457 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff() 458 @*/ 459 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) 460 { 461 Vec localX, localU; 462 PetscErrorCode ierr; 463 464 PetscFunctionBegin; 465 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 466 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 467 ierr = DMGetLocalVector(dm, &localU);CHKERRQ(ierr); 468 ierr = DMGlobalToLocalBegin(dm, U, INSERT_VALUES, localU);CHKERRQ(ierr); 469 ierr = DMGlobalToLocalEnd(dm, U, INSERT_VALUES, localU);CHKERRQ(ierr); 470 ierr = DMPlexProjectFieldLocal(dm, localU, funcs, mode, localX);CHKERRQ(ierr); 471 ierr = DMLocalToGlobalBegin(dm, localX, mode, X);CHKERRQ(ierr); 472 ierr = DMLocalToGlobalEnd(dm, localX, mode, X);CHKERRQ(ierr); 473 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 474 ierr = DMRestoreLocalVector(dm, &localU);CHKERRQ(ierr); 475 PetscFunctionReturn(0); 476 } 477 478 #undef __FUNCT__ 479 #define __FUNCT__ "DMPlexInsertBoundaryValues_FEM_Internal" 480 static PetscErrorCode DMPlexInsertBoundaryValues_FEM_Internal(DM dm, PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], void (*func)(const PetscReal[], PetscScalar *, void *), void *ctx, Vec locX) 481 { 482 void (**funcs)(const PetscReal x[], PetscScalar *u, void *ctx); 483 void **ctxs; 484 PetscInt numFields; 485 PetscErrorCode ierr; 486 487 PetscFunctionBegin; 488 ierr = DMGetNumFields(dm, &numFields);CHKERRQ(ierr); 489 ierr = PetscCalloc2(numFields,&funcs,numFields,&ctxs);CHKERRQ(ierr); 490 funcs[field] = func; 491 ctxs[field] = ctx; 492 ierr = DMPlexProjectFunctionLabelLocal(dm, label, numids, ids, funcs, ctxs, INSERT_BC_VALUES, locX);CHKERRQ(ierr); 493 ierr = PetscFree2(funcs,ctxs);CHKERRQ(ierr); 494 PetscFunctionReturn(0); 495 } 496 497 #undef __FUNCT__ 498 #define __FUNCT__ "DMPlexInsertBoundaryValues_FVM_Internal" 499 static PetscErrorCode DMPlexInsertBoundaryValues_FVM_Internal(DM dm, PetscReal time, Vec faceGeometry, Vec cellGeometry, Vec Grad, 500 PetscInt field, DMLabel label, PetscInt numids, const PetscInt ids[], PetscErrorCode (*func)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*), void *ctx, Vec locX) 501 { 502 PetscFV fv; 503 DM dmFace, dmCell, dmGrad; 504 const PetscScalar *facegeom, *cellgeom, *grad; 505 PetscScalar *x, *fx; 506 PetscInt dim, fStart, fEnd, pdim, i; 507 PetscErrorCode ierr; 508 509 PetscFunctionBegin; 510 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 511 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 512 ierr = DMGetField(dm, field, (PetscObject *) &fv);CHKERRQ(ierr); 513 ierr = VecGetDM(faceGeometry, &dmFace);CHKERRQ(ierr); 514 ierr = VecGetArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 515 ierr = VecGetDM(cellGeometry, &dmCell);CHKERRQ(ierr); 516 ierr = VecGetArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 517 if (Grad) { 518 ierr = VecGetDM(Grad, &dmGrad);CHKERRQ(ierr); 519 ierr = VecGetArrayRead(Grad, &grad);CHKERRQ(ierr); 520 ierr = PetscFVGetNumComponents(fv, &pdim);CHKERRQ(ierr); 521 ierr = DMGetWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 522 } 523 ierr = VecGetArray(locX, &x);CHKERRQ(ierr); 524 for (i = 0; i < numids; ++i) { 525 IS faceIS; 526 const PetscInt *faces; 527 PetscInt numFaces, f; 528 529 ierr = DMLabelGetStratumIS(label, ids[i], &faceIS);CHKERRQ(ierr); 530 if (!faceIS) continue; /* No points with that id on this process */ 531 ierr = ISGetLocalSize(faceIS, &numFaces);CHKERRQ(ierr); 532 ierr = ISGetIndices(faceIS, &faces);CHKERRQ(ierr); 533 for (f = 0; f < numFaces; ++f) { 534 const PetscInt face = faces[f], *cells; 535 const PetscFVFaceGeom *fg; 536 537 if ((face < fStart) || (face >= fEnd)) continue; /* Refinement adds non-faces to labels */ 538 ierr = DMPlexPointLocalRead(dmFace, face, facegeom, &fg);CHKERRQ(ierr); 539 ierr = DMPlexGetSupport(dm, face, &cells);CHKERRQ(ierr); 540 if (Grad) { 541 const PetscFVCellGeom *cg; 542 const PetscScalar *cx, *cgrad; 543 PetscScalar *xG; 544 PetscReal dx[3]; 545 PetscInt d; 546 547 ierr = DMPlexPointLocalRead(dmCell, cells[0], cellgeom, &cg);CHKERRQ(ierr); 548 ierr = DMPlexPointLocalRead(dm, cells[0], x, &cx);CHKERRQ(ierr); 549 ierr = DMPlexPointLocalRead(dmGrad, cells[0], grad, &cgrad);CHKERRQ(ierr); 550 ierr = DMPlexPointLocalRef(dm, cells[1], x, &xG);CHKERRQ(ierr); 551 DMPlex_WaxpyD_Internal(dim, -1, cg->centroid, fg->centroid, dx); 552 for (d = 0; d < pdim; ++d) fx[d] = cx[d] + DMPlex_DotD_Internal(dim, &cgrad[d*dim], dx); 553 ierr = (*func)(time, fg->centroid, fg->normal, fx, xG, ctx);CHKERRQ(ierr); 554 } else { 555 const PetscScalar *xI; 556 PetscScalar *xG; 557 558 ierr = DMPlexPointLocalRead(dm, cells[0], x, &xI);CHKERRQ(ierr); 559 ierr = DMPlexPointLocalRef(dm, cells[1], x, &xG);CHKERRQ(ierr); 560 ierr = (*func)(time, fg->centroid, fg->normal, xI, xG, ctx);CHKERRQ(ierr); 561 } 562 } 563 ierr = ISRestoreIndices(faceIS, &faces);CHKERRQ(ierr); 564 ierr = ISDestroy(&faceIS);CHKERRQ(ierr); 565 } 566 ierr = VecRestoreArray(locX, &x);CHKERRQ(ierr); 567 if (Grad) { 568 ierr = DMRestoreWorkArray(dm, pdim, PETSC_SCALAR, &fx);CHKERRQ(ierr); 569 ierr = VecRestoreArrayRead(Grad, &grad);CHKERRQ(ierr); 570 } 571 ierr = VecRestoreArrayRead(cellGeometry, &cellgeom);CHKERRQ(ierr); 572 ierr = VecRestoreArrayRead(faceGeometry, &facegeom);CHKERRQ(ierr); 573 PetscFunctionReturn(0); 574 } 575 576 #undef __FUNCT__ 577 #define __FUNCT__ "DMPlexInsertBoundaryValues" 578 PetscErrorCode DMPlexInsertBoundaryValues(DM dm, Vec locX, PetscReal time, Vec faceGeomFVM, Vec cellGeomFVM, Vec gradFVM) 579 { 580 PetscInt numBd, b; 581 PetscErrorCode ierr; 582 583 PetscFunctionBegin; 584 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 585 PetscValidHeaderSpecific(locX, VEC_CLASSID, 2); 586 if (faceGeomFVM) {PetscValidHeaderSpecific(faceGeomFVM, VEC_CLASSID, 4);} 587 if (cellGeomFVM) {PetscValidHeaderSpecific(cellGeomFVM, VEC_CLASSID, 5);} 588 if (gradFVM) {PetscValidHeaderSpecific(gradFVM, VEC_CLASSID, 6);} 589 ierr = DMPlexGetNumBoundary(dm, &numBd);CHKERRQ(ierr); 590 for (b = 0; b < numBd; ++b) { 591 PetscBool isEssential; 592 const char *labelname; 593 DMLabel label; 594 PetscInt field; 595 PetscObject obj; 596 PetscClassId id; 597 void (*func)(); 598 PetscInt numids; 599 const PetscInt *ids; 600 void *ctx; 601 602 ierr = DMPlexGetBoundary(dm, b, &isEssential, NULL, &labelname, &field, &func, &numids, &ids, &ctx);CHKERRQ(ierr); 603 if (!isEssential) continue; 604 ierr = DMPlexGetLabel(dm, labelname, &label);CHKERRQ(ierr); 605 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 606 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 607 if (id == PETSCFE_CLASSID) { 608 ierr = DMPlexInsertBoundaryValues_FEM_Internal(dm, field, label, numids, ids, (void (*)(const PetscReal[], PetscScalar *, void *)) func, ctx, locX);CHKERRQ(ierr); 609 } else if (id == PETSCFV_CLASSID) { 610 ierr = DMPlexInsertBoundaryValues_FVM_Internal(dm, time, faceGeomFVM, cellGeomFVM, gradFVM, 611 field, label, numids, ids, (PetscErrorCode (*)(PetscReal,const PetscReal*,const PetscReal*,const PetscScalar*,PetscScalar*,void*)) func, ctx, locX);CHKERRQ(ierr); 612 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 613 } 614 PetscFunctionReturn(0); 615 } 616 617 #undef __FUNCT__ 618 #define __FUNCT__ "DMPlexComputeL2Diff" 619 /*@C 620 DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h. 621 622 Input Parameters: 623 + dm - The DM 624 . funcs - The functions to evaluate for each field component 625 . ctxs - Optional array of contexts to pass to each function, or NULL. 626 - X - The coefficient vector u_h 627 628 Output Parameter: 629 . diff - The diff ||u - u_h||_2 630 631 Level: developer 632 633 .seealso: DMPlexProjectFunction(), DMPlexComputeL2FieldDiff(), DMPlexComputeL2GradientDiff() 634 @*/ 635 PetscErrorCode DMPlexComputeL2Diff(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff) 636 { 637 const PetscInt debug = 0; 638 PetscSection section; 639 PetscQuadrature quad; 640 Vec localX; 641 PetscScalar *funcVal, *interpolant; 642 PetscReal *coords, *v0, *J, *invJ, detJ; 643 PetscReal localDiff = 0.0; 644 const PetscReal *quadPoints, *quadWeights; 645 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, c, field, fieldOffset; 646 PetscErrorCode ierr; 647 648 PetscFunctionBegin; 649 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 650 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 651 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 652 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 653 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 654 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 655 for (field = 0; field < numFields; ++field) { 656 PetscObject obj; 657 PetscClassId id; 658 PetscInt Nc; 659 660 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 661 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 662 if (id == PETSCFE_CLASSID) { 663 PetscFE fe = (PetscFE) obj; 664 665 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 666 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 667 } else if (id == PETSCFV_CLASSID) { 668 PetscFV fv = (PetscFV) obj; 669 670 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 671 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 672 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 673 numComponents += Nc; 674 } 675 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 676 ierr = DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 677 ierr = PetscMalloc6(numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 678 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 679 for (c = cStart; c < cEnd; ++c) { 680 PetscScalar *x = NULL; 681 PetscReal elemDiff = 0.0; 682 683 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 684 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 685 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 686 687 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 688 PetscObject obj; 689 PetscClassId id; 690 void * const ctx = ctxs ? ctxs[field] : NULL; 691 PetscInt Nb, Nc, q, fc; 692 693 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 694 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 695 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 696 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 697 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 698 if (debug) { 699 char title[1024]; 700 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 701 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 702 } 703 for (q = 0; q < numQuadPoints; ++q) { 704 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 705 (*funcs[field])(coords, funcVal, ctx); 706 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 707 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 708 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 709 for (fc = 0; fc < Nc; ++fc) { 710 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);} 711 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 712 } 713 } 714 fieldOffset += Nb*Nc; 715 } 716 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 717 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 718 localDiff += elemDiff; 719 } 720 ierr = PetscFree6(funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 721 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 722 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 723 *diff = PetscSqrtReal(*diff); 724 PetscFunctionReturn(0); 725 } 726 727 #undef __FUNCT__ 728 #define __FUNCT__ "DMPlexComputeL2GradientDiff" 729 /*@C 730 DMPlexComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h. 731 732 Input Parameters: 733 + dm - The DM 734 . funcs - The gradient functions to evaluate for each field component 735 . ctxs - Optional array of contexts to pass to each function, or NULL. 736 . X - The coefficient vector u_h 737 - n - The vector to project along 738 739 Output Parameter: 740 . diff - The diff ||(grad u - grad u_h) . n||_2 741 742 Level: developer 743 744 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff() 745 @*/ 746 PetscErrorCode DMPlexComputeL2GradientDiff(DM dm, void (**funcs)(const PetscReal [], const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff) 747 { 748 const PetscInt debug = 0; 749 PetscSection section; 750 PetscQuadrature quad; 751 Vec localX; 752 PetscScalar *funcVal, *interpolantVec; 753 PetscReal *coords, *realSpaceDer, *v0, *J, *invJ, detJ; 754 PetscReal localDiff = 0.0; 755 PetscInt dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp; 756 PetscErrorCode ierr; 757 758 PetscFunctionBegin; 759 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 760 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 761 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 762 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 763 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 764 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 765 for (field = 0; field < numFields; ++field) { 766 PetscFE fe; 767 PetscInt Nc; 768 769 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 770 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 771 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 772 numComponents += Nc; 773 } 774 /* ierr = DMPlexProjectFunctionLocal(dm, fe, funcs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); */ 775 ierr = PetscMalloc7(numComponents,&funcVal,dim,&coords,dim,&realSpaceDer,dim,&v0,dim*dim,&J,dim*dim,&invJ,dim,&interpolantVec);CHKERRQ(ierr); 776 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 777 for (c = cStart; c < cEnd; ++c) { 778 PetscScalar *x = NULL; 779 PetscReal elemDiff = 0.0; 780 781 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 782 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 783 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 784 785 for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) { 786 PetscFE fe; 787 void * const ctx = ctxs ? ctxs[field] : NULL; 788 const PetscReal *quadPoints, *quadWeights; 789 PetscReal *basisDer; 790 PetscInt numQuadPoints, Nb, Ncomp, q, d, e, fc, f, g; 791 792 ierr = DMGetField(dm, field, (PetscObject *) &fe);CHKERRQ(ierr); 793 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 794 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 795 ierr = PetscFEGetNumComponents(fe, &Ncomp);CHKERRQ(ierr); 796 ierr = PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);CHKERRQ(ierr); 797 if (debug) { 798 char title[1024]; 799 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 800 ierr = DMPrintCellVector(c, title, Nb*Ncomp, &x[fieldOffset]);CHKERRQ(ierr); 801 } 802 for (q = 0; q < numQuadPoints; ++q) { 803 for (d = 0; d < dim; d++) { 804 coords[d] = v0[d]; 805 for (e = 0; e < dim; e++) { 806 coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0); 807 } 808 } 809 (*funcs[field])(coords, n, funcVal, ctx); 810 for (fc = 0; fc < Ncomp; ++fc) { 811 PetscScalar interpolant = 0.0; 812 813 for (d = 0; d < dim; ++d) interpolantVec[d] = 0.0; 814 for (f = 0; f < Nb; ++f) { 815 const PetscInt fidx = f*Ncomp+fc; 816 817 for (d = 0; d < dim; ++d) { 818 realSpaceDer[d] = 0.0; 819 for (g = 0; g < dim; ++g) { 820 realSpaceDer[d] += invJ[g*dim+d]*basisDer[(q*Nb*Ncomp+fidx)*dim+g]; 821 } 822 interpolantVec[d] += x[fieldOffset+fidx]*realSpaceDer[d]; 823 } 824 } 825 for (d = 0; d < dim; ++d) interpolant += interpolantVec[d]*n[d]; 826 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);} 827 elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal[fc]))*quadWeights[q]*detJ; 828 } 829 } 830 comp += Ncomp; 831 fieldOffset += Nb*Ncomp; 832 } 833 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 834 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);CHKERRQ(ierr);} 835 localDiff += elemDiff; 836 } 837 ierr = PetscFree7(funcVal,coords,realSpaceDer,v0,J,invJ,interpolantVec);CHKERRQ(ierr); 838 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 839 ierr = MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 840 *diff = PetscSqrtReal(*diff); 841 PetscFunctionReturn(0); 842 } 843 844 #undef __FUNCT__ 845 #define __FUNCT__ "DMPlexComputeL2FieldDiff" 846 /*@C 847 DMPlexComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components. 848 849 Input Parameters: 850 + dm - The DM 851 . funcs - The functions to evaluate for each field component 852 . ctxs - Optional array of contexts to pass to each function, or NULL. 853 - X - The coefficient vector u_h 854 855 Output Parameter: 856 . diff - The array of differences, ||u^f - u^f_h||_2 857 858 Level: developer 859 860 .seealso: DMPlexProjectFunction(), DMPlexComputeL2Diff(), DMPlexComputeL2GradientDiff() 861 @*/ 862 PetscErrorCode DMPlexComputeL2FieldDiff(DM dm, void (**funcs)(const PetscReal [], PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[]) 863 { 864 const PetscInt debug = 0; 865 PetscSection section; 866 PetscQuadrature quad; 867 Vec localX; 868 PetscScalar *funcVal, *interpolant; 869 PetscReal *coords, *v0, *J, *invJ, detJ; 870 PetscReal *localDiff; 871 const PetscReal *quadPoints, *quadWeights; 872 PetscInt dim, numFields, numComponents = 0, numQuadPoints, cStart, cEnd, c, field, fieldOffset; 873 PetscErrorCode ierr; 874 875 PetscFunctionBegin; 876 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 877 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 878 ierr = PetscSectionGetNumFields(section, &numFields);CHKERRQ(ierr); 879 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 880 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 881 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 882 for (field = 0; field < numFields; ++field) { 883 PetscObject obj; 884 PetscClassId id; 885 PetscInt Nc; 886 887 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 888 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 889 if (id == PETSCFE_CLASSID) { 890 PetscFE fe = (PetscFE) obj; 891 892 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr); 893 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 894 } else if (id == PETSCFV_CLASSID) { 895 PetscFV fv = (PetscFV) obj; 896 897 ierr = PetscFVGetQuadrature(fv, &quad);CHKERRQ(ierr); 898 ierr = PetscFVGetNumComponents(fv, &Nc);CHKERRQ(ierr); 899 } else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 900 numComponents += Nc; 901 } 902 ierr = PetscQuadratureGetData(quad, NULL, &numQuadPoints, &quadPoints, &quadWeights);CHKERRQ(ierr); 903 ierr = DMPlexProjectFunctionLocal(dm, funcs, ctxs, INSERT_BC_VALUES, localX);CHKERRQ(ierr); 904 ierr = PetscCalloc7(numFields,&localDiff,numComponents,&funcVal,numComponents,&interpolant,dim,&coords,dim,&v0,dim*dim,&J,dim*dim,&invJ);CHKERRQ(ierr); 905 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 906 for (c = cStart; c < cEnd; ++c) { 907 PetscScalar *x = NULL; 908 909 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ);CHKERRQ(ierr); 910 if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c); 911 ierr = DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 912 913 for (field = 0, fieldOffset = 0; field < numFields; ++field) { 914 PetscObject obj; 915 PetscClassId id; 916 void * const ctx = ctxs ? ctxs[field] : NULL; 917 PetscInt Nb, Nc, q, fc; 918 919 PetscReal elemDiff = 0.0; 920 921 ierr = DMGetField(dm, field, &obj);CHKERRQ(ierr); 922 ierr = PetscObjectGetClassId(obj, &id);CHKERRQ(ierr); 923 if (id == PETSCFE_CLASSID) {ierr = PetscFEGetNumComponents((PetscFE) obj, &Nc);CHKERRQ(ierr);ierr = PetscFEGetDimension((PetscFE) obj, &Nb);CHKERRQ(ierr);} 924 else if (id == PETSCFV_CLASSID) {ierr = PetscFVGetNumComponents((PetscFV) obj, &Nc);CHKERRQ(ierr);Nb = 1;} 925 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 926 if (debug) { 927 char title[1024]; 928 ierr = PetscSNPrintf(title, 1023, "Solution for Field %d", field);CHKERRQ(ierr); 929 ierr = DMPrintCellVector(c, title, Nb*Nc, &x[fieldOffset]);CHKERRQ(ierr); 930 } 931 for (q = 0; q < numQuadPoints; ++q) { 932 CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], coords); 933 (*funcs[field])(coords, funcVal, ctx); 934 if (id == PETSCFE_CLASSID) {ierr = PetscFEInterpolate_Static((PetscFE) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 935 else if (id == PETSCFV_CLASSID) {ierr = PetscFVInterpolate_Static((PetscFV) obj, &x[fieldOffset], q, interpolant);CHKERRQ(ierr);} 936 else SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", field); 937 for (fc = 0; fc < Nc; ++fc) { 938 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);} 939 elemDiff += PetscSqr(PetscRealPart(interpolant[fc] - funcVal[fc]))*quadWeights[q]*detJ; 940 } 941 } 942 fieldOffset += Nb*Nc; 943 localDiff[field] += elemDiff; 944 } 945 ierr = DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);CHKERRQ(ierr); 946 } 947 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 948 ierr = MPI_Allreduce(localDiff, diff, numFields, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));CHKERRQ(ierr); 949 for (field = 0; field < numFields; ++field) diff[field] = PetscSqrtReal(diff[field]); 950 ierr = PetscFree7(localDiff,funcVal,interpolant,coords,v0,J,invJ);CHKERRQ(ierr); 951 PetscFunctionReturn(0); 952 } 953 954 #undef __FUNCT__ 955 #define __FUNCT__ "DMPlexComputeIntegralFEM" 956 /*@ 957 DMPlexComputeIntegralFEM - Form the local integral F from the local input X using pointwise functions specified by the user 958 959 Input Parameters: 960 + dm - The mesh 961 . X - Local input vector 962 - user - The user context 963 964 Output Parameter: 965 . integral - Local integral for each field 966 967 Level: developer 968 969 .seealso: DMPlexComputeResidualFEM() 970 @*/ 971 PetscErrorCode DMPlexComputeIntegralFEM(DM dm, Vec X, PetscReal *integral, void *user) 972 { 973 DM_Plex *mesh = (DM_Plex *) dm->data; 974 DM dmAux; 975 Vec localX, A; 976 PetscDS prob, probAux = NULL; 977 PetscQuadrature q; 978 PetscSection section, sectionAux; 979 PetscFECellGeom *cgeom; 980 PetscScalar *u, *a = NULL; 981 PetscInt dim, Nf, f, numCells, cStart, cEnd, c; 982 PetscInt totDim, totDimAux; 983 PetscErrorCode ierr; 984 985 PetscFunctionBegin; 986 /*ierr = PetscLogEventBegin(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr);*/ 987 ierr = DMGetLocalVector(dm, &localX);CHKERRQ(ierr); 988 ierr = DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 989 ierr = DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);CHKERRQ(ierr); 990 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 991 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 992 ierr = DMGetDS(dm, &prob);CHKERRQ(ierr); 993 ierr = PetscDSGetTotalDimension(prob, &totDim);CHKERRQ(ierr); 994 ierr = PetscSectionGetNumFields(section, &Nf);CHKERRQ(ierr); 995 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 996 numCells = cEnd - cStart; 997 for (f = 0; f < Nf; ++f) {integral[f] = 0.0;} 998 ierr = PetscObjectQuery((PetscObject) dm, "dmAux", (PetscObject *) &dmAux);CHKERRQ(ierr); 999 ierr = PetscObjectQuery((PetscObject) dm, "A", (PetscObject *) &A);CHKERRQ(ierr); 1000 if (dmAux) { 1001 ierr = DMGetDefaultSection(dmAux, §ionAux);CHKERRQ(ierr); 1002 ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr); 1003 ierr = PetscDSGetTotalDimension(probAux, &totDimAux);CHKERRQ(ierr); 1004 } 1005 ierr = DMPlexInsertBoundaryValues(dm, localX, 0.0, NULL, NULL, NULL);CHKERRQ(ierr); 1006 ierr = PetscMalloc2(numCells*totDim,&u,numCells,&cgeom);CHKERRQ(ierr); 1007 if (dmAux) {ierr = PetscMalloc1(numCells*totDimAux, &a);CHKERRQ(ierr);} 1008 for (c = cStart; c < cEnd; ++c) { 1009 PetscScalar *x = NULL; 1010 PetscInt i; 1011 1012 ierr = DMPlexComputeCellGeometryFEM(dm, c, NULL, cgeom[c].v0, cgeom[c].J, cgeom[c].invJ, &cgeom[c].detJ);CHKERRQ(ierr); 1013 if (cgeom[c].detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", cgeom[c].detJ, c); 1014 ierr = DMPlexVecGetClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1015 for (i = 0; i < totDim; ++i) u[c*totDim+i] = x[i]; 1016 ierr = DMPlexVecRestoreClosure(dm, section, localX, c, NULL, &x);CHKERRQ(ierr); 1017 if (dmAux) { 1018 ierr = DMPlexVecGetClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1019 for (i = 0; i < totDimAux; ++i) a[c*totDimAux+i] = x[i]; 1020 ierr = DMPlexVecRestoreClosure(dmAux, sectionAux, A, c, NULL, &x);CHKERRQ(ierr); 1021 } 1022 } 1023 for (f = 0; f < Nf; ++f) { 1024 PetscFE fe; 1025 PetscInt numQuadPoints, Nb; 1026 /* Conforming batches */ 1027 PetscInt numChunks, numBatches, numBlocks, Ne, blockSize, batchSize; 1028 /* Remainder */ 1029 PetscInt Nr, offset; 1030 1031 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 1032 ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr); 1033 ierr = PetscFEGetDimension(fe, &Nb);CHKERRQ(ierr); 1034 ierr = PetscFEGetTileSizes(fe, NULL, &numBlocks, NULL, &numBatches);CHKERRQ(ierr); 1035 ierr = PetscQuadratureGetData(q, NULL, &numQuadPoints, NULL, NULL);CHKERRQ(ierr); 1036 blockSize = Nb*numQuadPoints; 1037 batchSize = numBlocks * blockSize; 1038 ierr = PetscFESetTileSizes(fe, blockSize, numBlocks, batchSize, numBatches);CHKERRQ(ierr); 1039 numChunks = numCells / (numBatches*batchSize); 1040 Ne = numChunks*numBatches*batchSize; 1041 Nr = numCells % (numBatches*batchSize); 1042 offset = numCells - Nr; 1043 ierr = PetscFEIntegrate(fe, prob, f, Ne, cgeom, u, probAux, a, integral);CHKERRQ(ierr); 1044 ierr = PetscFEIntegrate(fe, prob, f, Nr, &cgeom[offset], &u[offset*totDim], probAux, &a[offset*totDimAux], integral);CHKERRQ(ierr); 1045 } 1046 ierr = PetscFree2(u,cgeom);CHKERRQ(ierr); 1047 if (dmAux) {ierr = PetscFree(a);CHKERRQ(ierr);} 1048 if (mesh->printFEM) { 1049 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "Local integral:");CHKERRQ(ierr); 1050 for (f = 0; f < Nf; ++f) {ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), " %g", integral[f]);CHKERRQ(ierr);} 1051 ierr = PetscPrintf(PetscObjectComm((PetscObject) dm), "\n");CHKERRQ(ierr); 1052 } 1053 ierr = DMRestoreLocalVector(dm, &localX);CHKERRQ(ierr); 1054 /* TODO: Allreduce for integral */ 1055 /*ierr = PetscLogEventEnd(DMPLEX_IntegralFEM,dm,0,0,0);CHKERRQ(ierr);*/ 1056 PetscFunctionReturn(0); 1057 } 1058 1059 #undef __FUNCT__ 1060 #define __FUNCT__ "DMPlexComputeInterpolatorFEM" 1061 /*@ 1062 DMPlexComputeInterpolatorFEM - Form the local portion of the interpolation matrix I from the coarse DM to the uniformly refined DM. 1063 1064 Input Parameters: 1065 + dmf - The fine mesh 1066 . dmc - The coarse mesh 1067 - user - The user context 1068 1069 Output Parameter: 1070 . In - The interpolation matrix 1071 1072 Note: 1073 The first member of the user context must be an FEMContext. 1074 1075 We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator, 1076 like a GPU, or vectorize on a multicore machine. 1077 1078 Level: developer 1079 1080 .seealso: DMPlexComputeJacobianFEM() 1081 @*/ 1082 PetscErrorCode DMPlexComputeInterpolatorFEM(DM dmc, DM dmf, Mat In, void *user) 1083 { 1084 DM_Plex *mesh = (DM_Plex *) dmc->data; 1085 const char *name = "Interpolator"; 1086 PetscDS prob; 1087 PetscFE *feRef; 1088 PetscSection fsection, fglobalSection; 1089 PetscSection csection, cglobalSection; 1090 PetscScalar *elemMat; 1091 PetscInt dim, Nf, f, fieldI, fieldJ, offsetI, offsetJ, cStart, cEnd, c; 1092 PetscInt cTotDim, rTotDim = 0; 1093 PetscErrorCode ierr; 1094 1095 PetscFunctionBegin; 1096 ierr = PetscLogEventBegin(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1097 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1098 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1099 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1100 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1101 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1102 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1103 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1104 ierr = DMGetDS(dmf, &prob);CHKERRQ(ierr); 1105 ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr); 1106 for (f = 0; f < Nf; ++f) { 1107 PetscFE fe; 1108 PetscInt rNb, Nc; 1109 1110 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 1111 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1112 ierr = PetscFEGetDimension(feRef[f], &rNb);CHKERRQ(ierr); 1113 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1114 rTotDim += rNb*Nc; 1115 } 1116 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1117 ierr = PetscMalloc1(rTotDim*cTotDim,&elemMat);CHKERRQ(ierr); 1118 ierr = PetscMemzero(elemMat, rTotDim*cTotDim * sizeof(PetscScalar));CHKERRQ(ierr); 1119 for (fieldI = 0, offsetI = 0; fieldI < Nf; ++fieldI) { 1120 PetscDualSpace Qref; 1121 PetscQuadrature f; 1122 const PetscReal *qpoints, *qweights; 1123 PetscReal *points; 1124 PetscInt npoints = 0, Nc, Np, fpdim, i, k, p, d; 1125 1126 /* Compose points from all dual basis functionals */ 1127 ierr = PetscFEGetDualSpace(feRef[fieldI], &Qref);CHKERRQ(ierr); 1128 ierr = PetscFEGetNumComponents(feRef[fieldI], &Nc);CHKERRQ(ierr); 1129 ierr = PetscDualSpaceGetDimension(Qref, &fpdim);CHKERRQ(ierr); 1130 for (i = 0; i < fpdim; ++i) { 1131 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1132 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, NULL);CHKERRQ(ierr); 1133 npoints += Np; 1134 } 1135 ierr = PetscMalloc1(npoints*dim,&points);CHKERRQ(ierr); 1136 for (i = 0, k = 0; i < fpdim; ++i) { 1137 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1138 ierr = PetscQuadratureGetData(f, NULL, &Np, &qpoints, NULL);CHKERRQ(ierr); 1139 for (p = 0; p < Np; ++p, ++k) for (d = 0; d < dim; ++d) points[k*dim+d] = qpoints[p*dim+d]; 1140 } 1141 1142 for (fieldJ = 0, offsetJ = 0; fieldJ < Nf; ++fieldJ) { 1143 PetscFE fe; 1144 PetscReal *B; 1145 PetscInt NcJ, cpdim, j; 1146 1147 /* Evaluate basis at points */ 1148 ierr = PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);CHKERRQ(ierr); 1149 ierr = PetscFEGetNumComponents(fe, &NcJ);CHKERRQ(ierr); 1150 ierr = PetscFEGetDimension(fe, &cpdim);CHKERRQ(ierr); 1151 /* For now, fields only interpolate themselves */ 1152 if (fieldI == fieldJ) { 1153 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); 1154 ierr = PetscFEGetTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr); 1155 for (i = 0, k = 0; i < fpdim; ++i) { 1156 ierr = PetscDualSpaceGetFunctional(Qref, i, &f);CHKERRQ(ierr); 1157 ierr = PetscQuadratureGetData(f, NULL, &Np, NULL, &qweights);CHKERRQ(ierr); 1158 for (p = 0; p < Np; ++p, ++k) { 1159 for (j = 0; j < cpdim; ++j) { 1160 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]; 1161 } 1162 } 1163 } 1164 ierr = PetscFERestoreTabulation(fe, npoints, points, &B, NULL, NULL);CHKERRQ(ierr);CHKERRQ(ierr); 1165 } 1166 offsetJ += cpdim*NcJ; 1167 } 1168 offsetI += fpdim*Nc; 1169 ierr = PetscFree(points);CHKERRQ(ierr); 1170 } 1171 if (mesh->printFEM > 1) {ierr = DMPrintCellMatrix(0, name, rTotDim, cTotDim, elemMat);CHKERRQ(ierr);} 1172 /* Preallocate matrix */ 1173 { 1174 PetscHashJK ht; 1175 PetscLayout rLayout; 1176 PetscInt *dnz, *onz, *cellCIndices, *cellFIndices; 1177 PetscInt locRows, rStart, rEnd, cell, r; 1178 1179 ierr = MatGetLocalSize(In, &locRows, NULL);CHKERRQ(ierr); 1180 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) In), &rLayout);CHKERRQ(ierr); 1181 ierr = PetscLayoutSetLocalSize(rLayout, locRows);CHKERRQ(ierr); 1182 ierr = PetscLayoutSetBlockSize(rLayout, 1);CHKERRQ(ierr); 1183 ierr = PetscLayoutSetUp(rLayout);CHKERRQ(ierr); 1184 ierr = PetscLayoutGetRange(rLayout, &rStart, &rEnd);CHKERRQ(ierr); 1185 ierr = PetscLayoutDestroy(&rLayout);CHKERRQ(ierr); 1186 ierr = PetscCalloc4(locRows,&dnz,locRows,&onz,cTotDim,&cellCIndices,rTotDim,&cellFIndices);CHKERRQ(ierr); 1187 ierr = PetscHashJKCreate(&ht);CHKERRQ(ierr); 1188 for (cell = cStart; cell < cEnd; ++cell) { 1189 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, cell, cellCIndices, cellFIndices);CHKERRQ(ierr); 1190 for (r = 0; r < rTotDim; ++r) { 1191 PetscHashJKKey key; 1192 PetscHashJKIter missing, iter; 1193 1194 key.j = cellFIndices[r]; 1195 if (key.j < 0) continue; 1196 for (c = 0; c < cTotDim; ++c) { 1197 key.k = cellCIndices[c]; 1198 if (key.k < 0) continue; 1199 ierr = PetscHashJKPut(ht, key, &missing, &iter);CHKERRQ(ierr); 1200 if (missing) { 1201 ierr = PetscHashJKSet(ht, iter, 1);CHKERRQ(ierr); 1202 if ((key.k >= rStart) && (key.k < rEnd)) ++dnz[key.j-rStart]; 1203 else ++onz[key.j-rStart]; 1204 } 1205 } 1206 } 1207 } 1208 ierr = PetscHashJKDestroy(&ht);CHKERRQ(ierr); 1209 ierr = MatXAIJSetPreallocation(In, 1, dnz, onz, NULL, NULL);CHKERRQ(ierr); 1210 ierr = MatSetOption(In, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1211 ierr = PetscFree4(dnz,onz,cellCIndices,cellFIndices);CHKERRQ(ierr); 1212 } 1213 /* Fill matrix */ 1214 ierr = MatZeroEntries(In);CHKERRQ(ierr); 1215 for (c = cStart; c < cEnd; ++c) { 1216 ierr = DMPlexMatSetClosureRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, In, c, elemMat, INSERT_VALUES);CHKERRQ(ierr); 1217 } 1218 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1219 ierr = PetscFree(feRef);CHKERRQ(ierr); 1220 ierr = PetscFree(elemMat);CHKERRQ(ierr); 1221 ierr = MatAssemblyBegin(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1222 ierr = MatAssemblyEnd(In, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1223 if (mesh->printFEM) { 1224 ierr = PetscPrintf(PETSC_COMM_WORLD, "%s:\n", name);CHKERRQ(ierr); 1225 ierr = MatChop(In, 1.0e-10);CHKERRQ(ierr); 1226 ierr = MatView(In, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 1227 } 1228 ierr = PetscLogEventEnd(DMPLEX_InterpolatorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1229 PetscFunctionReturn(0); 1230 } 1231 1232 #undef __FUNCT__ 1233 #define __FUNCT__ "DMPlexComputeInjectorFEM" 1234 PetscErrorCode DMPlexComputeInjectorFEM(DM dmc, DM dmf, VecScatter *sc, void *user) 1235 { 1236 PetscDS prob; 1237 PetscFE *feRef; 1238 Vec fv, cv; 1239 IS fis, cis; 1240 PetscSection fsection, fglobalSection, csection, cglobalSection; 1241 PetscInt *cmap, *cellCIndices, *cellFIndices, *cindices, *findices; 1242 PetscInt cTotDim, fTotDim = 0, Nf, f, field, cStart, cEnd, c, dim, d, startC, offsetC, offsetF, m; 1243 PetscErrorCode ierr; 1244 1245 PetscFunctionBegin; 1246 ierr = PetscLogEventBegin(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1247 ierr = DMGetDimension(dmf, &dim);CHKERRQ(ierr); 1248 ierr = DMGetDefaultSection(dmf, &fsection);CHKERRQ(ierr); 1249 ierr = DMGetDefaultGlobalSection(dmf, &fglobalSection);CHKERRQ(ierr); 1250 ierr = DMGetDefaultSection(dmc, &csection);CHKERRQ(ierr); 1251 ierr = DMGetDefaultGlobalSection(dmc, &cglobalSection);CHKERRQ(ierr); 1252 ierr = PetscSectionGetNumFields(fsection, &Nf);CHKERRQ(ierr); 1253 ierr = DMPlexGetHeightStratum(dmc, 0, &cStart, &cEnd);CHKERRQ(ierr); 1254 ierr = DMGetDS(dmc, &prob);CHKERRQ(ierr); 1255 ierr = PetscMalloc1(Nf,&feRef);CHKERRQ(ierr); 1256 for (f = 0; f < Nf; ++f) { 1257 PetscFE fe; 1258 PetscInt fNb, Nc; 1259 1260 ierr = PetscDSGetDiscretization(prob, f, (PetscObject *) &fe);CHKERRQ(ierr); 1261 ierr = PetscFERefine(fe, &feRef[f]);CHKERRQ(ierr); 1262 ierr = PetscFEGetDimension(feRef[f], &fNb);CHKERRQ(ierr); 1263 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr); 1264 fTotDim += fNb*Nc; 1265 } 1266 ierr = PetscDSGetTotalDimension(prob, &cTotDim);CHKERRQ(ierr); 1267 ierr = PetscMalloc1(cTotDim,&cmap);CHKERRQ(ierr); 1268 for (field = 0, offsetC = 0, offsetF = 0; field < Nf; ++field) { 1269 PetscFE feC; 1270 PetscDualSpace QF, QC; 1271 PetscInt NcF, NcC, fpdim, cpdim; 1272 1273 ierr = PetscDSGetDiscretization(prob, field, (PetscObject *) &feC);CHKERRQ(ierr); 1274 ierr = PetscFEGetNumComponents(feC, &NcC);CHKERRQ(ierr); 1275 ierr = PetscFEGetNumComponents(feRef[field], &NcF);CHKERRQ(ierr); 1276 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); 1277 ierr = PetscFEGetDualSpace(feRef[field], &QF);CHKERRQ(ierr); 1278 ierr = PetscDualSpaceGetDimension(QF, &fpdim);CHKERRQ(ierr); 1279 ierr = PetscFEGetDualSpace(feC, &QC);CHKERRQ(ierr); 1280 ierr = PetscDualSpaceGetDimension(QC, &cpdim);CHKERRQ(ierr); 1281 for (c = 0; c < cpdim; ++c) { 1282 PetscQuadrature cfunc; 1283 const PetscReal *cqpoints; 1284 PetscInt NpC; 1285 1286 ierr = PetscDualSpaceGetFunctional(QC, c, &cfunc);CHKERRQ(ierr); 1287 ierr = PetscQuadratureGetData(cfunc, NULL, &NpC, &cqpoints, NULL);CHKERRQ(ierr); 1288 if (NpC != 1) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Do not know how to do injection for moments"); 1289 for (f = 0; f < fpdim; ++f) { 1290 PetscQuadrature ffunc; 1291 const PetscReal *fqpoints; 1292 PetscReal sum = 0.0; 1293 PetscInt NpF, comp; 1294 1295 ierr = PetscDualSpaceGetFunctional(QF, f, &ffunc);CHKERRQ(ierr); 1296 ierr = PetscQuadratureGetData(ffunc, NULL, &NpF, &fqpoints, NULL);CHKERRQ(ierr); 1297 if (NpC != NpF) continue; 1298 for (d = 0; d < dim; ++d) sum += PetscAbsReal(cqpoints[d] - fqpoints[d]); 1299 if (sum > 1.0e-9) continue; 1300 for (comp = 0; comp < NcC; ++comp) { 1301 cmap[(offsetC+c)*NcC+comp] = (offsetF+f)*NcF+comp; 1302 } 1303 break; 1304 } 1305 } 1306 offsetC += cpdim*NcC; 1307 offsetF += fpdim*NcF; 1308 } 1309 for (f = 0; f < Nf; ++f) {ierr = PetscFEDestroy(&feRef[f]);CHKERRQ(ierr);} 1310 ierr = PetscFree(feRef);CHKERRQ(ierr); 1311 1312 ierr = DMGetGlobalVector(dmf, &fv);CHKERRQ(ierr); 1313 ierr = DMGetGlobalVector(dmc, &cv);CHKERRQ(ierr); 1314 ierr = VecGetOwnershipRange(cv, &startC, NULL);CHKERRQ(ierr); 1315 ierr = PetscSectionGetConstrainedStorageSize(cglobalSection, &m);CHKERRQ(ierr); 1316 ierr = PetscMalloc2(cTotDim,&cellCIndices,fTotDim,&cellFIndices);CHKERRQ(ierr); 1317 ierr = PetscMalloc1(m,&cindices);CHKERRQ(ierr); 1318 ierr = PetscMalloc1(m,&findices);CHKERRQ(ierr); 1319 for (d = 0; d < m; ++d) cindices[d] = findices[d] = -1; 1320 for (c = cStart; c < cEnd; ++c) { 1321 ierr = DMPlexMatGetClosureIndicesRefined(dmf, fsection, fglobalSection, dmc, csection, cglobalSection, c, cellCIndices, cellFIndices);CHKERRQ(ierr); 1322 for (d = 0; d < cTotDim; ++d) { 1323 if (cellCIndices[d] < 0) continue; 1324 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]]); 1325 cindices[cellCIndices[d]-startC] = cellCIndices[d]; 1326 findices[cellCIndices[d]-startC] = cellFIndices[cmap[d]]; 1327 } 1328 } 1329 ierr = PetscFree(cmap);CHKERRQ(ierr); 1330 ierr = PetscFree2(cellCIndices,cellFIndices);CHKERRQ(ierr); 1331 1332 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, cindices, PETSC_OWN_POINTER, &cis);CHKERRQ(ierr); 1333 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, findices, PETSC_OWN_POINTER, &fis);CHKERRQ(ierr); 1334 ierr = VecScatterCreate(cv, cis, fv, fis, sc);CHKERRQ(ierr); 1335 ierr = ISDestroy(&cis);CHKERRQ(ierr); 1336 ierr = ISDestroy(&fis);CHKERRQ(ierr); 1337 ierr = DMRestoreGlobalVector(dmf, &fv);CHKERRQ(ierr); 1338 ierr = DMRestoreGlobalVector(dmc, &cv);CHKERRQ(ierr); 1339 ierr = PetscLogEventEnd(DMPLEX_InjectorFEM,dmc,dmf,0,0);CHKERRQ(ierr); 1340 PetscFunctionReturn(0); 1341 } 1342