Actual source code: plexfem.c
petsc-3.4.5 2014-06-29
1: #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/
5: PetscErrorCode DMPlexGetScale(DM dm, PetscUnit unit, PetscReal *scale)
6: {
7: DM_Plex *mesh = (DM_Plex*) dm->data;
12: *scale = mesh->scale[unit];
13: return(0);
14: }
18: PetscErrorCode DMPlexSetScale(DM dm, PetscUnit unit, PetscReal scale)
19: {
20: DM_Plex *mesh = (DM_Plex*) dm->data;
24: mesh->scale[unit] = scale;
25: return(0);
26: }
28: PETSC_STATIC_INLINE PetscInt epsilon(PetscInt i, PetscInt j, PetscInt k)
29: {
30: switch (i) {
31: case 0:
32: switch (j) {
33: case 0: return 0;
34: case 1:
35: switch (k) {
36: case 0: return 0;
37: case 1: return 0;
38: case 2: return 1;
39: }
40: case 2:
41: switch (k) {
42: case 0: return 0;
43: case 1: return -1;
44: case 2: return 0;
45: }
46: }
47: case 1:
48: switch (j) {
49: case 0:
50: switch (k) {
51: case 0: return 0;
52: case 1: return 0;
53: case 2: return -1;
54: }
55: case 1: return 0;
56: case 2:
57: switch (k) {
58: case 0: return 1;
59: case 1: return 0;
60: case 2: return 0;
61: }
62: }
63: case 2:
64: switch (j) {
65: case 0:
66: switch (k) {
67: case 0: return 0;
68: case 1: return 1;
69: case 2: return 0;
70: }
71: case 1:
72: switch (k) {
73: case 0: return -1;
74: case 1: return 0;
75: case 2: return 0;
76: }
77: case 2: return 0;
78: }
79: }
80: return 0;
81: }
85: /*@C
86: DMPlexCreateRigidBody - create rigid body modes from coordinates
88: Collective on DM
90: Input Arguments:
91: + dm - the DM
92: . section - the local section associated with the rigid field, or NULL for the default section
93: - globalSection - the global section associated with the rigid field, or NULL for the default section
95: Output Argument:
96: . sp - the null space
98: Note: This is necessary to take account of Dirichlet conditions on the displacements
100: Level: advanced
102: .seealso: MatNullSpaceCreate()
103: @*/
104: PetscErrorCode DMPlexCreateRigidBody(DM dm, PetscSection section, PetscSection globalSection, MatNullSpace *sp)
105: {
106: MPI_Comm comm;
107: Vec coordinates, localMode, mode[6];
108: PetscSection coordSection;
109: PetscScalar *coords;
110: PetscInt dim, vStart, vEnd, v, n, m, d, i, j;
114: PetscObjectGetComm((PetscObject)dm,&comm);
115: DMPlexGetDimension(dm, &dim);
116: if (dim == 1) {
117: MatNullSpaceCreate(comm, PETSC_TRUE, 0, NULL, sp);
118: return(0);
119: }
120: if (!section) {DMGetDefaultSection(dm, §ion);}
121: if (!globalSection) {DMGetDefaultGlobalSection(dm, &globalSection);}
122: PetscSectionGetConstrainedStorageSize(globalSection, &n);
123: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
124: DMPlexGetCoordinateSection(dm, &coordSection);
125: DMGetCoordinatesLocal(dm, &coordinates);
126: m = (dim*(dim+1))/2;
127: VecCreate(comm, &mode[0]);
128: VecSetSizes(mode[0], n, PETSC_DETERMINE);
129: VecSetUp(mode[0]);
130: for (i = 1; i < m; ++i) {VecDuplicate(mode[0], &mode[i]);}
131: /* Assume P1 */
132: DMGetLocalVector(dm, &localMode);
133: for (d = 0; d < dim; ++d) {
134: PetscScalar values[3] = {0.0, 0.0, 0.0};
136: values[d] = 1.0;
137: VecSet(localMode, 0.0);
138: for (v = vStart; v < vEnd; ++v) {
139: DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);
140: }
141: DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);
142: DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);
143: }
144: VecGetArray(coordinates, &coords);
145: for (d = dim; d < dim*(dim+1)/2; ++d) {
146: PetscInt i, j, k = dim > 2 ? d - dim : d;
148: VecSet(localMode, 0.0);
149: for (v = vStart; v < vEnd; ++v) {
150: PetscScalar values[3] = {0.0, 0.0, 0.0};
151: PetscInt off;
153: PetscSectionGetOffset(coordSection, v, &off);
154: for (i = 0; i < dim; ++i) {
155: for (j = 0; j < dim; ++j) {
156: values[j] += epsilon(i, j, k)*PetscRealPart(coords[off+i]);
157: }
158: }
159: DMPlexVecSetClosure(dm, section, localMode, v, values, INSERT_VALUES);
160: }
161: DMLocalToGlobalBegin(dm, localMode, INSERT_VALUES, mode[d]);
162: DMLocalToGlobalEnd(dm, localMode, INSERT_VALUES, mode[d]);
163: }
164: VecRestoreArray(coordinates, &coords);
165: DMRestoreLocalVector(dm, &localMode);
166: for (i = 0; i < dim; ++i) {VecNormalize(mode[i], NULL);}
167: /* Orthonormalize system */
168: for (i = dim; i < m; ++i) {
169: PetscScalar dots[6];
171: VecMDot(mode[i], i, mode, dots);
172: for (j = 0; j < i; ++j) dots[j] *= -1.0;
173: VecMAXPY(mode[i], i, dots, mode);
174: VecNormalize(mode[i], NULL);
175: }
176: MatNullSpaceCreate(comm, PETSC_FALSE, m, mode, sp);
177: for (i = 0; i< m; ++i) {VecDestroy(&mode[i]);}
178: return(0);
179: }
180: /*******************************************************************************
181: This should be in a separate Discretization object, but I am not sure how to lay
182: it out yet, so I am stuffing things here while I experiment.
183: *******************************************************************************/
186: PetscErrorCode DMPlexSetFEMIntegration(DM dm,
187: PetscErrorCode (*integrateResidualFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[],
188: const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
189: void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
190: void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]),
191: PetscErrorCode (*integrateBdResidualFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[],
192: const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
193: void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]),
194: void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]), PetscScalar[]),
195: PetscErrorCode (*integrateJacobianActionFEM)(PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[], const PetscScalar[],
196: const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
197: void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
198: void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
199: void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
200: void (**)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]),
201: PetscErrorCode (*integrateJacobianFEM)(PetscInt, PetscInt, PetscInt, PetscInt, PetscQuadrature[], const PetscScalar[],
202: const PetscReal[], const PetscReal[], const PetscReal[], const PetscReal[],
203: void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
204: void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
205: void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]),
206: void (*)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]), PetscScalar[]))
207: {
208: DM_Plex *mesh = (DM_Plex*) dm->data;
212: mesh->integrateResidualFEM = integrateResidualFEM;
213: mesh->integrateBdResidualFEM = integrateBdResidualFEM;
214: mesh->integrateJacobianActionFEM = integrateJacobianActionFEM;
215: mesh->integrateJacobianFEM = integrateJacobianFEM;
216: return(0);
217: }
221: PetscErrorCode DMPlexProjectFunctionLocal(DM dm, PetscInt numComp, void (**funcs)(const PetscReal [], PetscScalar *), InsertMode mode, Vec localX)
222: {
223: Vec coordinates;
224: PetscSection section, cSection;
225: PetscInt dim, vStart, vEnd, v, c, d;
226: PetscScalar *values, *cArray;
227: PetscReal *coords;
231: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
232: DMGetDefaultSection(dm, §ion);
233: DMPlexGetCoordinateSection(dm, &cSection);
234: DMGetCoordinatesLocal(dm, &coordinates);
235: PetscMalloc(numComp * sizeof(PetscScalar), &values);
236: VecGetArray(coordinates, &cArray);
237: PetscSectionGetDof(cSection, vStart, &dim);
238: PetscMalloc(dim * sizeof(PetscReal),&coords);
239: for (v = vStart; v < vEnd; ++v) {
240: PetscInt dof, off;
242: PetscSectionGetDof(cSection, v, &dof);
243: PetscSectionGetOffset(cSection, v, &off);
244: if (dof > dim) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Cannot have more coordinates %d then dimensions %d", dof, dim);
245: for (d = 0; d < dof; ++d) coords[d] = PetscRealPart(cArray[off+d]);
246: for (c = 0; c < numComp; ++c) (*funcs[c])(coords, &values[c]);
247: VecSetValuesSection(localX, section, v, values, mode);
248: }
249: VecRestoreArray(coordinates, &cArray);
250: /* Temporary, must be replaced by a projection on the finite element basis */
251: {
252: PetscInt eStart = 0, eEnd = 0, e, depth;
254: DMPlexGetLabelSize(dm, "depth", &depth);
255: --depth;
256: if (depth > 1) {DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);}
257: for (e = eStart; e < eEnd; ++e) {
258: const PetscInt *cone = NULL;
259: PetscInt coneSize, d;
260: PetscScalar *coordsA, *coordsB;
262: DMPlexGetConeSize(dm, e, &coneSize);
263: DMPlexGetCone(dm, e, &cone);
264: if (coneSize != 2) SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_SIZ, "Cone size %d for point %d should be 2", coneSize, e);
265: VecGetValuesSection(coordinates, cSection, cone[0], &coordsA);
266: VecGetValuesSection(coordinates, cSection, cone[1], &coordsB);
267: for (d = 0; d < dim; ++d) {
268: coords[d] = 0.5*(PetscRealPart(coordsA[d]) + PetscRealPart(coordsB[d]));
269: }
270: for (c = 0; c < numComp; ++c) (*funcs[c])(coords, &values[c]);
271: VecSetValuesSection(localX, section, e, values, mode);
272: }
273: }
275: PetscFree(coords);
276: PetscFree(values);
277: #if 0
278: const PetscInt localDof = this->_mesh->sizeWithBC(s, *cells->begin());
279: PetscReal detJ;
281: PetscMalloc(localDof * sizeof(PetscScalar), &values);
282: PetscMalloc2(dim,PetscReal,&v0,dim*dim,PetscReal,&J);
283: ALE::ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type> pV(PetscPowInt(this->_mesh->getSieve()->getMaxConeSize(),dim+1), true);
285: for (PetscInt c = cStart; c < cEnd; ++c) {
286: ALE::ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*this->_mesh->getSieve(), c, pV);
287: const PETSC_MESH_TYPE::point_type *oPoints = pV.getPoints();
288: const int oSize = pV.getSize();
289: int v = 0;
291: DMPlexComputeCellGeometry(dm, c, v0, J, NULL, &detJ);
292: for (PetscInt cl = 0; cl < oSize; ++cl) {
293: const PetscInt fDim;
295: PetscSectionGetDof(oPoints[cl], &fDim);
296: if (pointDim) {
297: for (PetscInt d = 0; d < fDim; ++d, ++v) {
298: values[v] = (*this->_options.integrate)(v0, J, v, initFunc);
299: }
300: }
301: }
302: DMPlexVecSetClosure(dm, NULL, localX, c, values);
303: pV.clear();
304: }
305: PetscFree2(v0,J);
306: PetscFree(values);
307: #endif
308: return(0);
309: }
313: /*@C
314: DMPlexProjectFunction - This projects the given function into the function space provided.
316: Input Parameters:
317: + dm - The DM
318: . numComp - The number of components (functions)
319: . funcs - The coordinate functions to evaluate
320: - mode - The insertion mode for values
322: Output Parameter:
323: . X - vector
325: Level: developer
327: Note:
328: This currently just calls the function with the coordinates of each vertex and edge midpoint, and stores the result in a vector.
329: We will eventually fix it.
331: .seealso: DMPlexComputeL2Diff()
332: @*/
333: PetscErrorCode DMPlexProjectFunction(DM dm, PetscInt numComp, void (**funcs)(const PetscReal [], PetscScalar *), InsertMode mode, Vec X)
334: {
335: Vec localX;
339: DMGetLocalVector(dm, &localX);
340: DMPlexProjectFunctionLocal(dm, numComp, funcs, mode, localX);
341: DMLocalToGlobalBegin(dm, localX, mode, X);
342: DMLocalToGlobalEnd(dm, localX, mode, X);
343: DMRestoreLocalVector(dm, &localX);
344: return(0);
345: }
349: /*@C
350: DMPlexComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.
352: Input Parameters:
353: + dm - The DM
354: . quad - The PetscQuadrature object for each field
355: . funcs - The functions to evaluate for each field component
356: - X - The coefficient vector u_h
358: Output Parameter:
359: . diff - The diff ||u - u_h||_2
361: Level: developer
363: .seealso: DMPlexProjectFunction()
364: @*/
365: PetscErrorCode DMPlexComputeL2Diff(DM dm, PetscQuadrature quad[], void (**funcs)(const PetscReal [], PetscScalar *), Vec X, PetscReal *diff)
366: {
367: const PetscInt debug = 0;
368: PetscSection section;
369: Vec localX;
370: PetscReal *coords, *v0, *J, *invJ, detJ;
371: PetscReal localDiff = 0.0;
372: PetscInt dim, numFields, numComponents = 0, cStart, cEnd, c, field, fieldOffset, comp;
376: DMPlexGetDimension(dm, &dim);
377: DMGetDefaultSection(dm, §ion);
378: PetscSectionGetNumFields(section, &numFields);
379: DMGetLocalVector(dm, &localX);
380: DMGlobalToLocalBegin(dm, X, INSERT_VALUES, localX);
381: DMGlobalToLocalEnd(dm, X, INSERT_VALUES, localX);
382: for (field = 0; field < numFields; ++field) {
383: numComponents += quad[field].numComponents;
384: }
385: DMPlexProjectFunctionLocal(dm, numComponents, funcs, INSERT_BC_VALUES, localX);
386: PetscMalloc4(dim,PetscReal,&coords,dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);
387: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
388: for (c = cStart; c < cEnd; ++c) {
389: PetscScalar *x;
390: PetscReal elemDiff = 0.0;
392: DMPlexComputeCellGeometry(dm, c, v0, J, invJ, &detJ);
393: if (detJ <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);
394: DMPlexVecGetClosure(dm, NULL, localX, c, NULL, &x);
396: for (field = 0, comp = 0, fieldOffset = 0; field < numFields; ++field) {
397: const PetscInt numQuadPoints = quad[field].numQuadPoints;
398: const PetscReal *quadPoints = quad[field].quadPoints;
399: const PetscReal *quadWeights = quad[field].quadWeights;
400: const PetscInt numBasisFuncs = quad[field].numBasisFuncs;
401: const PetscInt numBasisComps = quad[field].numComponents;
402: const PetscReal *basis = quad[field].basis;
403: PetscInt q, d, e, fc, f;
405: if (debug) {
406: char title[1024];
407: PetscSNPrintf(title, 1023, "Solution for Field %d", field);
408: DMPrintCellVector(c, title, numBasisFuncs*numBasisComps, &x[fieldOffset]);
409: }
410: for (q = 0; q < numQuadPoints; ++q) {
411: for (d = 0; d < dim; d++) {
412: coords[d] = v0[d];
413: for (e = 0; e < dim; e++) {
414: coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
415: }
416: }
417: for (fc = 0; fc < numBasisComps; ++fc) {
418: PetscScalar funcVal;
419: PetscScalar interpolant = 0.0;
421: (*funcs[comp+fc])(coords, &funcVal);
422: for (f = 0; f < numBasisFuncs; ++f) {
423: const PetscInt fidx = f*numBasisComps+fc;
424: interpolant += x[fieldOffset+fidx]*basis[q*numBasisFuncs*numBasisComps+fidx];
425: }
426: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d field %d diff %g\n", c, field, PetscSqr(PetscRealPart(interpolant - funcVal))*quadWeights[q]*detJ);}
427: elemDiff += PetscSqr(PetscRealPart(interpolant - funcVal))*quadWeights[q]*detJ;
428: }
429: }
430: comp += numBasisComps;
431: fieldOffset += numBasisFuncs*numBasisComps;
432: }
433: DMPlexVecRestoreClosure(dm, NULL, localX, c, NULL, &x);
434: if (debug) {PetscPrintf(PETSC_COMM_SELF, " elem %d diff %g\n", c, elemDiff);}
435: localDiff += elemDiff;
436: }
437: PetscFree4(coords,v0,J,invJ);
438: DMRestoreLocalVector(dm, &localX);
439: MPI_Allreduce(&localDiff, diff, 1, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)dm));
440: *diff = PetscSqrtReal(*diff);
441: return(0);
442: }
446: /*@
447: DMPlexComputeResidualFEM - Form the local residual F from the local input X using pointwise functions specified by the user
449: Input Parameters:
450: + dm - The mesh
451: . X - Local input vector
452: - user - The user context
454: Output Parameter:
455: . F - Local output vector
457: Note:
458: The second member of the user context must be an FEMContext.
460: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
461: like a GPU, or vectorize on a multicore machine.
463: Level: developer
465: .seealso: DMPlexComputeJacobianActionFEM()
466: @*/
467: PetscErrorCode DMPlexComputeResidualFEM(DM dm, Vec X, Vec F, void *user)
468: {
469: DM_Plex *mesh = (DM_Plex*) dm->data;
470: PetscFEM *fem = (PetscFEM*) &((DM*) user)[1];
471: PetscQuadrature *quad = fem->quad;
472: PetscQuadrature *quadBd = fem->quadBd;
473: PetscSection section;
474: PetscReal *v0, *n, *J, *invJ, *detJ;
475: PetscScalar *elemVec, *u;
476: PetscInt dim, numFields, field, numBatchesTmp = 1, numCells, cStart, cEnd, c;
477: PetscInt cellDof, numComponents;
478: PetscBool has;
479: PetscErrorCode ierr;
482: /* PetscLogEventBegin(ResidualFEMEvent,0,0,0,0); */
483: DMPlexGetDimension(dm, &dim);
484: DMGetDefaultSection(dm, §ion);
485: PetscSectionGetNumFields(section, &numFields);
486: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
487: numCells = cEnd - cStart;
488: for (field = 0, cellDof = 0, numComponents = 0; field < numFields; ++field) {
489: cellDof += quad[field].numBasisFuncs*quad[field].numComponents;
490: numComponents += quad[field].numComponents;
491: }
492: DMPlexProjectFunctionLocal(dm, numComponents, fem->bcFuncs, INSERT_BC_VALUES, X);
493: VecSet(F, 0.0);
494: PetscMalloc6(numCells*cellDof,PetscScalar,&u,numCells*dim,PetscReal,&v0,numCells*dim*dim,PetscReal,&J,numCells*dim*dim,PetscReal,&invJ,numCells,PetscReal,&detJ,numCells*cellDof,PetscScalar,&elemVec);
495: for (c = cStart; c < cEnd; ++c) {
496: PetscScalar *x;
497: PetscInt i;
499: DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);
500: if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
501: DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);
503: for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
504: DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);
505: }
506: for (field = 0; field < numFields; ++field) {
507: const PetscInt numQuadPoints = quad[field].numQuadPoints;
508: const PetscInt numBasisFuncs = quad[field].numBasisFuncs;
509: void (*f0)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f0Funcs[field];
510: void (*f1)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->f1Funcs[field];
511: /* Conforming batches */
512: PetscInt blockSize = numBasisFuncs*numQuadPoints;
513: PetscInt numBlocks = 1;
514: PetscInt batchSize = numBlocks * blockSize;
515: PetscInt numBatches = numBatchesTmp;
516: PetscInt numChunks = numCells / (numBatches*batchSize);
517: /* Remainder */
518: PetscInt numRemainder = numCells % (numBatches * batchSize);
519: PetscInt offset = numCells - numRemainder;
521: (*mesh->integrateResidualFEM)(numChunks*numBatches*batchSize, numFields, field, quad, u, v0, J, invJ, detJ, f0, f1, elemVec);
522: (*mesh->integrateResidualFEM)(numRemainder, numFields, field, quad, &u[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset],
523: f0, f1, &elemVec[offset*cellDof]);
524: }
525: for (c = cStart; c < cEnd; ++c) {
526: if (mesh->printFEM > 1) {DMPrintCellVector(c, "Residual", cellDof, &elemVec[c*cellDof]);}
527: DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);
528: }
529: PetscFree6(u,v0,J,invJ,detJ,elemVec);
530: /* Integration over the boundary:
531: - This can probably be generalized to integration over a set of labels, however
532: the idea here is to do integration where we need the cell normal
533: - We can replace hardcoding with a registration process, and this is how we hook
534: up the system to something like FEniCS
535: */
536: DMPlexHasLabel(dm, "boundary", &has);
537: if (has && quadBd) {
538: DMLabel label;
539: IS pointIS;
540: const PetscInt *points;
541: PetscInt numPoints, p;
543: DMPlexGetLabel(dm, "boundary", &label);
544: DMLabelGetStratumSize(label, 1, &numPoints);
545: DMLabelGetStratumIS(label, 1, &pointIS);
546: ISGetIndices(pointIS, &points);
547: for (field = 0, cellDof = 0, numComponents = 0; field < numFields; ++field) {
548: cellDof += quadBd[field].numBasisFuncs*quadBd[field].numComponents;
549: numComponents += quadBd[field].numComponents;
550: }
551: PetscMalloc7(numPoints*cellDof,PetscScalar,&u,numPoints*dim,PetscReal,&v0,numPoints*dim,PetscReal,&n,numPoints*dim*dim,PetscReal,&J,numPoints*dim*dim,PetscReal,&invJ,numPoints,PetscReal,&detJ,numPoints*cellDof,PetscScalar,&elemVec);
552: for (p = 0; p < numPoints; ++p) {
553: const PetscInt point = points[p];
554: PetscScalar *x;
555: PetscInt i;
557: /* TODO: Add normal determination here */
558: DMPlexComputeCellGeometry(dm, point, &v0[p*dim], &J[p*dim*dim], &invJ[p*dim*dim], &detJ[p]);
559: if (detJ[p] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for face %d", detJ[p], point);
560: DMPlexVecGetClosure(dm, NULL, X, point, NULL, &x);
562: for (i = 0; i < cellDof; ++i) u[p*cellDof+i] = x[i];
563: DMPlexVecRestoreClosure(dm, NULL, X, point, NULL, &x);
564: }
565: for (field = 0; field < numFields; ++field) {
566: const PetscInt numQuadPoints = quadBd[field].numQuadPoints;
567: const PetscInt numBasisFuncs = quadBd[field].numBasisFuncs;
568: void (*f0)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f0BdFuncs[field];
569: void (*f1)(const PetscScalar[], const PetscScalar[], const PetscReal[], const PetscReal[], PetscScalar[]) = fem->f1BdFuncs[field];
570: /* Conforming batches */
571: PetscInt blockSize = numBasisFuncs*numQuadPoints;
572: PetscInt numBlocks = 1;
573: PetscInt batchSize = numBlocks * blockSize;
574: PetscInt numBatches = numBatchesTmp;
575: PetscInt numChunks = numPoints / (numBatches*batchSize);
576: /* Remainder */
577: PetscInt numRemainder = numPoints % (numBatches * batchSize);
578: PetscInt offset = numPoints - numRemainder;
580: (*mesh->integrateBdResidualFEM)(numChunks*numBatches*batchSize, numFields, field, quadBd, u, v0, n, J, invJ, detJ, f0, f1, elemVec);
581: (*mesh->integrateBdResidualFEM)(numRemainder, numFields, field, quadBd, &u[offset*cellDof], &v0[offset*dim], &n[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset],
582: f0, f1, &elemVec[offset*cellDof]);
583: }
584: for (p = 0; p < numPoints; ++p) {
585: const PetscInt point = points[p];
587: if (mesh->printFEM > 1) {DMPrintCellVector(point, "Residual", cellDof, &elemVec[p*cellDof]);}
588: DMPlexVecSetClosure(dm, NULL, F, point, &elemVec[p*cellDof], ADD_VALUES);
589: }
590: ISRestoreIndices(pointIS, &points);
591: ISDestroy(&pointIS);
592: PetscFree7(u,v0,n,J,invJ,detJ,elemVec);
593: }
594: if (mesh->printFEM) {
595: PetscMPIInt rank, numProcs;
596: PetscInt p;
598: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
599: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);
600: PetscPrintf(PetscObjectComm((PetscObject)dm), "Residual:\n");
601: for (p = 0; p < numProcs; ++p) {
602: if (p == rank) {
603: Vec f;
605: VecDuplicate(F, &f);
606: VecCopy(F, f);
607: VecChop(f, 1.0e-10);
608: VecView(f, PETSC_VIEWER_STDOUT_SELF);
609: VecDestroy(&f);
610: PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);
611: }
612: PetscBarrier((PetscObject) dm);
613: }
614: }
615: /* PetscLogEventEnd(ResidualFEMEvent,0,0,0,0); */
616: return(0);
617: }
621: /*@C
622: DMPlexComputeJacobianActionFEM - Form the local action of Jacobian J(u) on the local input X using pointwise functions specified by the user
624: Input Parameters:
625: + dm - The mesh
626: . J - The Jacobian shell matrix
627: . X - Local input vector
628: - user - The user context
630: Output Parameter:
631: . F - Local output vector
633: Note:
634: The second member of the user context must be an FEMContext.
636: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
637: like a GPU, or vectorize on a multicore machine.
639: Level: developer
641: .seealso: DMPlexComputeResidualFEM()
642: @*/
643: PetscErrorCode DMPlexComputeJacobianActionFEM(DM dm, Mat Jac, Vec X, Vec F, void *user)
644: {
645: DM_Plex *mesh = (DM_Plex*) dm->data;
646: PetscFEM *fem = (PetscFEM*) &((DM*) user)[1];
647: PetscQuadrature *quad = fem->quad;
648: PetscSection section;
649: JacActionCtx *jctx;
650: PetscReal *v0, *J, *invJ, *detJ;
651: PetscScalar *elemVec, *u, *a;
652: PetscInt dim, numFields, field, numBatchesTmp = 1, numCells, cStart, cEnd, c;
653: PetscInt cellDof = 0;
654: PetscErrorCode ierr;
657: /* PetscLogEventBegin(JacobianActionFEMEvent,0,0,0,0); */
658: MatShellGetContext(Jac, &jctx);
659: DMPlexGetDimension(dm, &dim);
660: DMGetDefaultSection(dm, §ion);
661: PetscSectionGetNumFields(section, &numFields);
662: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
663: numCells = cEnd - cStart;
664: for (field = 0; field < numFields; ++field) {
665: cellDof += quad[field].numBasisFuncs*quad[field].numComponents;
666: }
667: VecSet(F, 0.0);
668: PetscMalloc7(numCells*cellDof,PetscScalar,&u,numCells*cellDof,PetscScalar,&a,numCells*dim,PetscReal,&v0,numCells*dim*dim,PetscReal,&J,numCells*dim*dim,PetscReal,&invJ,numCells,PetscReal,&detJ,numCells*cellDof,PetscScalar,&elemVec);
669: for (c = cStart; c < cEnd; ++c) {
670: PetscScalar *x;
671: PetscInt i;
673: DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);
674: if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
675: DMPlexVecGetClosure(dm, NULL, jctx->u, c, NULL, &x);
676: for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
677: DMPlexVecRestoreClosure(dm, NULL, jctx->u, c, NULL, &x);
678: DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);
679: for (i = 0; i < cellDof; ++i) a[c*cellDof+i] = x[i];
680: DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);
681: }
682: for (field = 0; field < numFields; ++field) {
683: const PetscInt numQuadPoints = quad[field].numQuadPoints;
684: const PetscInt numBasisFuncs = quad[field].numBasisFuncs;
685: /* Conforming batches */
686: PetscInt blockSize = numBasisFuncs*numQuadPoints;
687: PetscInt numBlocks = 1;
688: PetscInt batchSize = numBlocks * blockSize;
689: PetscInt numBatches = numBatchesTmp;
690: PetscInt numChunks = numCells / (numBatches*batchSize);
691: /* Remainder */
692: PetscInt numRemainder = numCells % (numBatches * batchSize);
693: PetscInt offset = numCells - numRemainder;
695: (*mesh->integrateJacobianActionFEM)(numChunks*numBatches*batchSize, numFields, field, quad, u, a, v0, J, invJ, detJ, fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, elemVec);
696: (*mesh->integrateJacobianActionFEM)(numRemainder, numFields, field, quad, &u[offset*cellDof], &a[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset],
697: fem->g0Funcs, fem->g1Funcs, fem->g2Funcs, fem->g3Funcs, &elemVec[offset*cellDof]);
698: }
699: for (c = cStart; c < cEnd; ++c) {
700: if (mesh->printFEM > 1) {DMPrintCellVector(c, "Jacobian Action", cellDof, &elemVec[c*cellDof]);}
701: DMPlexVecSetClosure(dm, NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);
702: }
703: PetscFree7(u,a,v0,J,invJ,detJ,elemVec);
704: if (mesh->printFEM) {
705: PetscMPIInt rank, numProcs;
706: PetscInt p;
708: MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
709: MPI_Comm_size(PetscObjectComm((PetscObject)dm), &numProcs);
710: PetscPrintf(PetscObjectComm((PetscObject)dm), "Jacobian Action:\n");
711: for (p = 0; p < numProcs; ++p) {
712: if (p == rank) {VecView(F, PETSC_VIEWER_STDOUT_SELF);}
713: PetscBarrier((PetscObject) dm);
714: }
715: }
716: /* PetscLogEventEnd(JacobianActionFEMEvent,0,0,0,0); */
717: return(0);
718: }
722: /*@
723: DMPlexComputeJacobianFEM - Form the local portion of the Jacobian matrix J at the local solution X using pointwise functions specified by the user.
725: Input Parameters:
726: + dm - The mesh
727: . X - Local input vector
728: - user - The user context
730: Output Parameter:
731: . Jac - Jacobian matrix
733: Note:
734: The second member of the user context must be an FEMContext.
736: We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
737: like a GPU, or vectorize on a multicore machine.
739: Level: developer
741: .seealso: FormFunctionLocal()
742: @*/
743: PetscErrorCode DMPlexComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat JacP, MatStructure *str,void *user)
744: {
745: DM_Plex *mesh = (DM_Plex*) dm->data;
746: PetscFEM *fem = (PetscFEM*) &((DM*) user)[1];
747: PetscQuadrature *quad = fem->quad;
748: PetscSection section;
749: PetscReal *v0, *J, *invJ, *detJ;
750: PetscScalar *elemMat, *u;
751: PetscInt dim, numFields, field, fieldI, numBatchesTmp = 1, numCells, cStart, cEnd, c;
752: PetscInt cellDof = 0, numComponents = 0;
753: PetscBool isShell;
754: PetscErrorCode ierr;
757: /* PetscLogEventBegin(JacobianFEMEvent,0,0,0,0); */
758: DMPlexGetDimension(dm, &dim);
759: DMGetDefaultSection(dm, §ion);
760: PetscSectionGetNumFields(section, &numFields);
761: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
762: numCells = cEnd - cStart;
763: for (field = 0; field < numFields; ++field) {
764: cellDof += quad[field].numBasisFuncs*quad[field].numComponents;
765: numComponents += quad[field].numComponents;
766: }
767: DMPlexProjectFunctionLocal(dm, numComponents, fem->bcFuncs, INSERT_BC_VALUES, X);
768: MatZeroEntries(JacP);
769: PetscMalloc6(numCells*cellDof,PetscScalar,&u,numCells*dim,PetscReal,&v0,numCells*dim*dim,PetscReal,&J,numCells*dim*dim,PetscReal,&invJ,numCells,PetscReal,&detJ,numCells*cellDof*cellDof,PetscScalar,&elemMat);
770: for (c = cStart; c < cEnd; ++c) {
771: PetscScalar *x;
772: PetscInt i;
774: DMPlexComputeCellGeometry(dm, c, &v0[c*dim], &J[c*dim*dim], &invJ[c*dim*dim], &detJ[c]);
775: if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
776: DMPlexVecGetClosure(dm, NULL, X, c, NULL, &x);
778: for (i = 0; i < cellDof; ++i) u[c*cellDof+i] = x[i];
779: DMPlexVecRestoreClosure(dm, NULL, X, c, NULL, &x);
780: }
781: PetscMemzero(elemMat, numCells*cellDof*cellDof * sizeof(PetscScalar));
782: for (fieldI = 0; fieldI < numFields; ++fieldI) {
783: const PetscInt numQuadPoints = quad[fieldI].numQuadPoints;
784: const PetscInt numBasisFuncs = quad[fieldI].numBasisFuncs;
785: PetscInt fieldJ;
787: for (fieldJ = 0; fieldJ < numFields; ++fieldJ) {
788: void (*g0)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g0Funcs[fieldI*numFields+fieldJ];
789: void (*g1)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g1Funcs[fieldI*numFields+fieldJ];
790: void (*g2)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g2Funcs[fieldI*numFields+fieldJ];
791: void (*g3)(const PetscScalar[], const PetscScalar[], const PetscReal[], PetscScalar[]) = fem->g3Funcs[fieldI*numFields+fieldJ];
792: /* Conforming batches */
793: PetscInt blockSize = numBasisFuncs*numQuadPoints;
794: PetscInt numBlocks = 1;
795: PetscInt batchSize = numBlocks * blockSize;
796: PetscInt numBatches = numBatchesTmp;
797: PetscInt numChunks = numCells / (numBatches*batchSize);
798: /* Remainder */
799: PetscInt numRemainder = numCells % (numBatches * batchSize);
800: PetscInt offset = numCells - numRemainder;
802: (*mesh->integrateJacobianFEM)(numChunks*numBatches*batchSize, numFields, fieldI, fieldJ, quad, u, v0, J, invJ, detJ, g0, g1, g2, g3, elemMat);
803: (*mesh->integrateJacobianFEM)(numRemainder, numFields, fieldI, fieldJ, quad, &u[offset*cellDof], &v0[offset*dim], &J[offset*dim*dim], &invJ[offset*dim*dim], &detJ[offset],
804: g0, g1, g2, g3, &elemMat[offset*cellDof*cellDof]);
805: }
806: }
807: for (c = cStart; c < cEnd; ++c) {
808: if (mesh->printFEM > 1) {DMPrintCellMatrix(c, "Jacobian", cellDof, cellDof, &elemMat[c*cellDof*cellDof]);}
809: DMPlexMatSetClosure(dm, NULL, NULL, JacP, c, &elemMat[c*cellDof*cellDof], ADD_VALUES);
810: }
811: PetscFree6(u,v0,J,invJ,detJ,elemMat);
813: /* Assemble matrix, using the 2-step process:
814: MatAssemblyBegin(), MatAssemblyEnd(). */
815: MatAssemblyBegin(JacP, MAT_FINAL_ASSEMBLY);
816: MatAssemblyEnd(JacP, MAT_FINAL_ASSEMBLY);
818: if (mesh->printFEM) {
819: PetscPrintf(PETSC_COMM_WORLD, "Jacobian:\n");
820: MatChop(JacP, 1.0e-10);
821: MatView(JacP, PETSC_VIEWER_STDOUT_WORLD);
822: }
823: /* PetscLogEventEnd(JacobianFEMEvent,0,0,0,0); */
824: PetscObjectTypeCompare((PetscObject)Jac, MATSHELL, &isShell);
825: if (isShell) {
826: JacActionCtx *jctx;
828: MatShellGetContext(Jac, &jctx);
829: VecCopy(X, jctx->u);
830: }
831: *str = SAME_NONZERO_PATTERN;
832: return(0);
833: }