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