Actual source code: ex52.c
1: static char help[] = "Testbed for FEM operations on the GPU\n\n";
3: #include<petscdmmesh.h>
4: #include<petscsnes.h>
6: #ifndef PETSC_HAVE_SIEVE
7: #error This example requires Sieve. Reconfigure using --with-sieve
8: #endif
10: typedef struct {
11: DM dm; /* REQUIRED in order to use SNES evaluation functions */
12: PetscInt debug; /* The debugging level */
13: PetscMPIInt rank; /* The process rank */
14: PetscMPIInt numProcs; /* The number of processes */
15: PetscInt dim; /* The topological mesh dimension */
16: PetscBool interpolate; /* Generate intermediate mesh elements */
17: PetscReal refinementLimit; /* The largest allowable cell volume */
18: char partitioner[2048]; /* The graph partitioner */
19: PetscBool computeFunction; /* The flag for computing a residual */
20: PetscBool computeJacobian; /* The flag for computing a Jacobian */
21: PetscBool batch; /* The flag for batch assembly */
22: /* Element quadrature */
23: PetscQuadrature q;
24: } AppCtx;
26: /*------------------------------------------------------------------------------
27: This code can be generated using 'bin/pythonscripts/PetscGenerateFEMQuadrature.py 2 1 src/snes/examples/tutorials/ex52.h'
28: -----------------------------------------------------------------------------*/
29: #include "ex52.h"
33: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) {
37: options->debug = 0;
38: options->dim = 2;
39: options->interpolate = PETSC_FALSE;
40: options->refinementLimit = 0.0;
41: options->computeFunction = PETSC_FALSE;
42: options->computeJacobian = PETSC_FALSE;
43: options->batch = PETSC_FALSE;
45: MPI_Comm_size(comm, &options->numProcs);
46: MPI_Comm_rank(comm, &options->rank);
47: PetscOptionsBegin(comm, "", "Bratu Problem Options", "DMMESH");
48: PetscOptionsInt("-debug", "The debugging level", "ex52.c", options->debug, &options->debug, PETSC_NULL);
49: PetscOptionsInt("-dim", "The topological mesh dimension", "ex52.c", options->dim, &options->dim, PETSC_NULL);
50: PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex52.c", options->interpolate, &options->interpolate, PETSC_NULL);
51: PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex52.c", options->refinementLimit, &options->refinementLimit, PETSC_NULL);
52: PetscStrcpy(options->partitioner, "chaco");
53: PetscOptionsString("-partitioner", "The graph partitioner", "pflotran.cxx", options->partitioner, options->partitioner, 2048, PETSC_NULL);
54: PetscOptionsBool("-compute_function", "Compute the residual", "ex52.c", options->computeFunction, &options->computeFunction, PETSC_NULL);
55: PetscOptionsBool("-compute_jacobian", "Compute the Jacobian", "ex52.c", options->computeJacobian, &options->computeJacobian, PETSC_NULL);
56: PetscOptionsBool("-batch", "Use the batch assembly method", "ex52.c", options->batch, &options->batch, PETSC_NULL);
57: PetscOptionsEnd();
58: return(0);
59: };
63: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
64: {
65: PetscInt dim = user->dim;
66: PetscBool interpolate = user->interpolate;
67: PetscReal refinementLimit = user->refinementLimit;
68: const char *partitioner = user->partitioner;
72: DMMeshCreateBoxMesh(comm, dim, interpolate, dm);
73: {
74: DM refinedMesh = PETSC_NULL;
75: DM distributedMesh = PETSC_NULL;
77: /* Refine mesh using a volume constraint */
78: DMMeshRefine(*dm, refinementLimit, interpolate, &refinedMesh);
79: if (refinedMesh) {
80: DMDestroy(dm);
81: *dm = refinedMesh;
82: }
83: /* Distribute mesh over processes */
84: DMMeshDistribute(*dm, partitioner, &distributedMesh);
85: if (distributedMesh) {
86: DMDestroy(dm);
87: *dm = distributedMesh;
88: }
89: }
90: DMSetFromOptions(*dm);
91: user->dm = *dm;
92: return(0);
93: }
97: PetscErrorCode SetupQuadrature(AppCtx *user) {
99: switch(user->dim) {
100: case 2:
101: user->q.numQuadPoints = NUM_QUADRATURE_POINTS_0;
102: user->q.quadPoints = points_0;
103: user->q.quadWeights = weights_0;
104: user->q.numBasisFuncs = NUM_BASIS_FUNCTIONS_0;
105: user->q.basis = Basis_0;
106: user->q.basisDer = BasisDerivatives_0;
107: break;
108: default:
109: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
110: }
111: return(0);
112: }
116: PetscErrorCode SetupSection(DM dm, AppCtx *user) {
117: PetscSection section;
118: /* These can be generated using config/PETSc/FEM.py */
119: PetscInt numDof_0[2] = {1, 0};
120: PetscInt numDof_1[3] = {1, 0, 0};
121: PetscInt numDof_2[4] = {1, 0, 0, 0};
122: PetscInt dim = user->dim;
123: PetscInt *numDof;
124: const char *bcLabel = PETSC_NULL;
128: switch(user->dim) {
129: case 1:
130: numDof = numDof_0;
131: break;
132: case 2:
133: numDof = numDof_1;
134: break;
135: case 3:
136: numDof = numDof_2;
137: break;
138: default:
139: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid spatial dimension %d", user->dim);
140: }
141: //if (user->bcType == DIRICHLET) {
142: // bcLabel = "marker";
143: //}
144: DMMeshCreateSection(dm, dim, numDof, bcLabel, 1, §ion);
145: DMMeshSetSection(dm, "default", section);
146: return(0);
147: }
151: /*
152: dm - The mesh
153: X - Local intput vector
154: F - Local output vector
155: */
156: PetscErrorCode FormFunctionLocal(DM dm, Vec X, Vec F, AppCtx *user)
157: {
158: //PetscScalar (*rhsFunc)(const PetscReal []) = user->rhsFunc;
159: const PetscInt debug = user->debug;
160: const PetscInt dim = user->dim;
161: const PetscInt numQuadPoints = user->q.numQuadPoints;
162: const PetscReal *quadPoints = user->q.quadPoints;
163: const PetscReal *quadWeights = user->q.quadWeights;
164: const PetscInt numBasisFuncs = user->q.numBasisFuncs;
165: const PetscReal *basis = user->q.basis;
166: const PetscReal *basisDer = user->q.basisDer;
167: PetscReal *coords, *v0, *J, *invJ, detJ;
168: PetscScalar *realSpaceDer, *fieldGrad, *elemVec;
169: PetscInt cStart, cEnd;
170: PetscErrorCode ierr;
173: VecSet(F, 0.0);
174: PetscMalloc3(dim,PetscScalar,&realSpaceDer,dim,PetscScalar,&fieldGrad,numBasisFuncs,PetscScalar,&elemVec);
175: PetscMalloc4(dim,PetscReal,&coords,dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);
176: DMMeshGetHeightStratum(dm, 0, &cStart, &cEnd);
177: for(PetscInt c = cStart; c < cEnd; ++c) {
178: const PetscScalar *x;
180: PetscMemzero(elemVec, numBasisFuncs * sizeof(PetscScalar));
181: DMMeshComputeCellGeometry(dm, c, v0, J, invJ, &detJ);
182: if (detJ <= 0.0) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);}
183: DMMeshVecGetClosure(dm, X, c, &x);
184: if (debug) {DMMeshPrintCellVector(c, "Solution", numBasisFuncs, x);}
186: for(int q = 0; q < numQuadPoints; ++q) {
187: PetscScalar fieldVal = 0.0;
189: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
190: for(int d = 0; d < dim; ++d) {
191: fieldGrad[d] = 0.0;
192: coords[d] = v0[d];
193: for(int e = 0; e < dim; ++e) {
194: coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
195: }
196: if (debug) {PetscPrintf(PETSC_COMM_SELF, " coords[%d] %g\n", d, coords[d]);}
197: }
198: for(int f = 0; f < numBasisFuncs; ++f) {
199: fieldVal += x[f]*basis[q*numBasisFuncs+f];
201: for(int d = 0; d < dim; ++d) {
202: realSpaceDer[d] = 0.0;
203: for(int e = 0; e < dim; ++e) {
204: realSpaceDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+f)*dim+e];
205: }
206: fieldGrad[d] += realSpaceDer[d]*x[f];
207: }
208: }
209: if (debug) {
210: for(int d = 0; d < dim; ++d) {
211: PetscPrintf(PETSC_COMM_SELF, " fieldGrad[%d] %g\n", d, fieldGrad[d]);
212: }
213: }
214: const PetscScalar funcVal = 0.0; //(*rhsFunc)(coords);
215: for(int f = 0; f < numBasisFuncs; ++f) {
216: /* Constant term: -f(x) */
217: elemVec[f] -= basis[q*numBasisFuncs+f]*funcVal*quadWeights[q]*detJ;
218: /* Linear term: -\Delta u */
219: PetscScalar product = 0.0;
220: for(int d = 0; d < dim; ++d) {
221: realSpaceDer[d] = 0.0;
222: for(int e = 0; e < dim; ++e) {
223: realSpaceDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+f)*dim+e];
224: }
225: product += realSpaceDer[d]*fieldGrad[d];
226: }
227: elemVec[f] += product*quadWeights[q]*detJ;
228: /* Nonlinear term: -\lambda e^{u} */
229: elemVec[f] -= basis[q*numBasisFuncs+f]*0.0*PetscExpScalar(fieldVal)*quadWeights[q]*detJ;
230: }
231: }
232: if (debug) {DMMeshPrintCellVector(c, "Residual", numBasisFuncs, elemVec);}
233: DMMeshVecSetClosure(dm, F, c, elemVec, ADD_VALUES);
234: }
235: PetscLogFlops((cEnd-cStart)*numQuadPoints*numBasisFuncs*(dim*(dim*5+4)+14));
236: PetscFree3(realSpaceDer,fieldGrad,elemVec);
237: PetscFree4(coords,v0,J,invJ);
239: PetscPrintf(PETSC_COMM_WORLD, "Residual:\n");
240: for(int p = 0; p < user->numProcs; ++p) {
241: if (p == user->rank) {VecView(F, PETSC_VIEWER_STDOUT_SELF);}
242: PetscBarrier((PetscObject) dm);
243: }
244: return(0);
245: }
251: /*
252: dm - The mesh
253: X - Local intput vector
254: F - Local output vector
255: */
256: PetscErrorCode FormFunctionLocalBatch(DM dm, Vec X, Vec F, AppCtx *user)
257: {
258: //PetscScalar (*rhsFunc)(const PetscReal []) = user->rhsFunc;
259: const PetscInt debug = user->debug;
260: const PetscInt dim = user->dim;
261: const PetscInt numQuadPoints = user->q.numQuadPoints;
262: const PetscReal *quadPoints = user->q.quadPoints;
263: const PetscReal *quadWeights = user->q.quadWeights;
264: const PetscInt numBasisFuncs = user->q.numBasisFuncs;
265: const PetscReal *basis = user->q.basis;
266: const PetscReal *basisDer = user->q.basisDer;
267: PetscReal *coords, *v0, *J, *invJ, *detJ;
268: PetscScalar *realSpaceDer, *fieldGrad, *elemVec;
269: PetscInt cStart, cEnd;
270: PetscErrorCode ierr;
273: VecSet(F, 0.0);
274: //PetscMalloc3(dim,PetscScalar,&realSpaceDer,dim,PetscScalar,&fieldGrad,numBasisFuncs,PetscScalar,&elemVec);
275: PetscMalloc4(dim,PetscReal,&coords,dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);
276: DMMeshGetHeightStratum(dm, 0, &cStart, &cEnd);
277: const PetscInt numCells = cEnd - cStart;
278: PetscScalar *u;
280: PetscMalloc3(numCells*numBasisFuncs,PetscScalar,&u,numCells,PetscReal,&detJ,numCells*numBasisFuncs,PetscScalar,&elemVec);
281: for(PetscInt c = cStart; c < cEnd; ++c) {
282: const PetscScalar *x;
284: DMMeshComputeCellGeometry(dm, c, v0, J, invJ, &detJ[c]);
285: DMMeshVecGetClosure(dm, X, c, &x);
287: for(int f = 0; f < numBasisFuncs; ++f) {
288: u[c*numBasisFuncs+f] = x[f];
289: }
290: }
291: IntegrateElementBatch(numCells, numBasisFuncs, u, detJ, numQuadPoints, quadPoints, quadWeights, basis, basisDer, elemVec);
292: for(PetscInt c = cStart; c < cEnd; ++c) {
293: DMMeshVecSetClosure(dm, F, c, &elemVec[c*numBasisFuncs], ADD_VALUES);
294: }
295: PetscFree3(u,detJ,elemVec);
297: PetscLogFlops(0);
298: //PetscFree3(realSpaceDer,fieldGrad,elemVec);
299: PetscFree4(coords,v0,J,invJ);
301: PetscPrintf(PETSC_COMM_WORLD, "Residual:\n");
302: for(int p = 0; p < user->numProcs; ++p) {
303: if (p == user->rank) {VecView(F, PETSC_VIEWER_STDOUT_SELF);}
304: PetscBarrier((PetscObject) dm);
305: }
306: return(0);
307: }
311: /*
312: dm - The mesh
313: X - The local input vector
314: Jac - The output matrix
315: */
316: PetscErrorCode FormJacobianLocal(DM dm, Vec X, Mat Jac, AppCtx *user)
317: {
318: const PetscInt debug = user->debug;
319: const PetscInt dim = user->dim;
320: const PetscInt numQuadPoints = user->q.numQuadPoints;
321: const PetscReal *quadWeights = user->q.quadWeights;
322: const PetscInt numBasisFuncs = user->q.numBasisFuncs;
323: const PetscReal *basis = user->q.basis;
324: const PetscReal *basisDer = user->q.basisDer;
325: PetscReal *v0, *J, *invJ, detJ;
326: PetscScalar *realSpaceTestDer, *realSpaceBasisDer, *elemMat;
327: PetscInt cStart, cEnd;
328: PetscErrorCode ierr;
331: MatZeroEntries(Jac);
332: PetscMalloc3(dim,PetscScalar,&realSpaceTestDer,dim,PetscScalar,&realSpaceBasisDer,numBasisFuncs*numBasisFuncs,PetscScalar,&elemMat);
333: PetscMalloc3(dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);
334: DMMeshGetHeightStratum(dm, 0, &cStart, &cEnd);
335: for(PetscInt c = cStart; c < cEnd; ++c) {
336: const PetscScalar *x;
338: PetscMemzero(elemMat, numBasisFuncs*numBasisFuncs * sizeof(PetscScalar));
339: DMMeshComputeCellGeometry(dm, c, v0, J, invJ, &detJ);
340: if (detJ <= 0.0) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);}
341: DMMeshVecGetClosure(dm, X, c, &x);
343: for(int q = 0; q < numQuadPoints; ++q) {
344: PetscScalar fieldVal = 0.0;
346: for(int f = 0; f < numBasisFuncs; ++f) {
347: fieldVal += x[f]*basis[q*numBasisFuncs+f];
348: }
349: for(int f = 0; f < numBasisFuncs; ++f) {
350: for(int d = 0; d < dim; ++d) {
351: realSpaceTestDer[d] = 0.0;
352: for(int e = 0; e < dim; ++e) {
353: realSpaceTestDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+f)*dim+e];
354: }
355: }
356: for(int g = 0; g < numBasisFuncs; ++g) {
357: for(int d = 0; d < dim; ++d) {
358: realSpaceBasisDer[d] = 0.0;
359: for(int e = 0; e < dim; ++e) {
360: realSpaceBasisDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+g)*dim+e];
361: }
362: }
363: /* Linear term: -\Delta u */
364: PetscScalar product = 0.0;
365: for(int d = 0; d < dim; ++d) product += realSpaceTestDer[d]*realSpaceBasisDer[d];
366: elemMat[f*numBasisFuncs+g] += product*quadWeights[q]*detJ;
367: /* Nonlinear term: -\lambda e^{u} */
368: elemMat[f*numBasisFuncs+g] -= basis[q*numBasisFuncs+f]*basis[q*numBasisFuncs+g]*0.0*PetscExpScalar(fieldVal)*quadWeights[q]*detJ;
369: }
370: }
371: }
372: if (debug) {DMMeshPrintCellMatrix(c, "Jacobian", numBasisFuncs, numBasisFuncs, elemMat);}
373: DMMeshMatSetClosure(dm, Jac, c, elemMat, ADD_VALUES);
374: }
375: PetscLogFlops((cEnd-cStart)*numQuadPoints*numBasisFuncs*(dim*(dim*5+4)+14));
376: PetscFree3(realSpaceTestDer,realSpaceBasisDer,elemMat);
377: PetscFree3(v0,J,invJ);
379: /* Assemble matrix, using the 2-step process:
380: MatAssemblyBegin(), MatAssemblyEnd(). */
381: MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);
382: MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);
383: /* Tell the matrix we will never add a new nonzero location to the
384: matrix. If we do, it will generate an error. */
385: MatSetOption(Jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
386: return(0);
387: }
391: /*
392: dm - The mesh
393: X - The local input vector
394: Jac - The output matrix
395: */
396: PetscErrorCode FormJacobianLocalBatch(DM dm, Vec X, Mat Jac, AppCtx *user)
397: {
398: const PetscInt debug = user->debug;
399: const PetscInt dim = user->dim;
400: const PetscInt numQuadPoints = user->q.numQuadPoints;
401: const PetscReal *quadWeights = user->q.quadWeights;
402: const PetscInt numBasisFuncs = user->q.numBasisFuncs;
403: const PetscReal *basis = user->q.basis;
404: const PetscReal *basisDer = user->q.basisDer;
405: PetscReal *v0, *J, *invJ, detJ;
406: PetscScalar *realSpaceTestDer, *realSpaceBasisDer, *elemMat;
407: PetscInt cStart, cEnd;
408: PetscErrorCode ierr;
411: MatZeroEntries(Jac);
412: PetscMalloc3(dim,PetscScalar,&realSpaceTestDer,dim,PetscScalar,&realSpaceBasisDer,numBasisFuncs*numBasisFuncs,PetscScalar,&elemMat);
413: PetscMalloc3(dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);
414: DMMeshGetHeightStratum(dm, 0, &cStart, &cEnd);
415: for(PetscInt c = cStart; c < cEnd; ++c) {
416: const PetscScalar *x;
418: PetscMemzero(elemMat, numBasisFuncs*numBasisFuncs * sizeof(PetscScalar));
419: DMMeshComputeCellGeometry(dm, c, v0, J, invJ, &detJ);
420: if (detJ <= 0.0) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);}
421: DMMeshVecGetClosure(dm, X, c, &x);
423: for(int q = 0; q < numQuadPoints; ++q) {
424: PetscScalar fieldVal = 0.0;
426: for(int f = 0; f < numBasisFuncs; ++f) {
427: fieldVal += x[f]*basis[q*numBasisFuncs+f];
428: }
429: for(int f = 0; f < numBasisFuncs; ++f) {
430: for(int d = 0; d < dim; ++d) {
431: realSpaceTestDer[d] = 0.0;
432: for(int e = 0; e < dim; ++e) {
433: realSpaceTestDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+f)*dim+e];
434: }
435: }
436: for(int g = 0; g < numBasisFuncs; ++g) {
437: for(int d = 0; d < dim; ++d) {
438: realSpaceBasisDer[d] = 0.0;
439: for(int e = 0; e < dim; ++e) {
440: realSpaceBasisDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+g)*dim+e];
441: }
442: }
443: /* Linear term: -\Delta u */
444: PetscScalar product = 0.0;
445: for(int d = 0; d < dim; ++d) product += realSpaceTestDer[d]*realSpaceBasisDer[d];
446: elemMat[f*numBasisFuncs+g] += product*quadWeights[q]*detJ;
447: /* Nonlinear term: -\lambda e^{u} */
448: elemMat[f*numBasisFuncs+g] -= basis[q*numBasisFuncs+f]*basis[q*numBasisFuncs+g]*0.0*PetscExpScalar(fieldVal)*quadWeights[q]*detJ;
449: }
450: }
451: }
452: if (debug) {DMMeshPrintCellMatrix(c, "Jacobian", numBasisFuncs, numBasisFuncs, elemMat);}
453: DMMeshMatSetClosure(dm, Jac, c, elemMat, ADD_VALUES);
454: }
455: PetscLogFlops((cEnd-cStart)*numQuadPoints*numBasisFuncs*(dim*(dim*5+4)+14));
456: PetscFree3(realSpaceTestDer,realSpaceBasisDer,elemMat);
457: PetscFree3(v0,J,invJ);
459: /* Assemble matrix, using the 2-step process:
460: MatAssemblyBegin(), MatAssemblyEnd(). */
461: MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);
462: MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);
463: /* Tell the matrix we will never add a new nonzero location to the
464: matrix. If we do, it will generate an error. */
465: MatSetOption(Jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
466: return(0);
467: }
471: int main(int argc, char **argv)
472: {
473: DM dm;
474: SNES snes;
475: AppCtx user;
478: PetscInitialize(&argc, &argv, PETSC_NULL, help);
479: #ifndef PETSC_HAVE_CUDA
480: SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "This example requires CUDA support.");
481: #endif
482: ProcessOptions(PETSC_COMM_WORLD, &user);
483: CreateMesh(PETSC_COMM_WORLD, &user, &dm);
484: SetupQuadrature(&user);
485: SetupSection(dm, &user);
487: SNESCreate(PETSC_COMM_WORLD, &snes);
488: if (user.computeFunction) {
489: Vec X, F;
491: DMGetGlobalVector(dm, &X);
492: DMGetGlobalVector(dm, &F);
493: if (user.batch) {
494: DMMeshSetLocalFunction(dm, (PetscErrorCode (*)(DM, Vec, Vec, void*)) FormFunctionLocalBatch);
495: } else {
496: DMMeshSetLocalFunction(dm, (PetscErrorCode (*)(DM, Vec, Vec, void*)) FormFunctionLocal);
497: }
498: SNESMeshFormFunction(snes, X, F, &user);
499: DMRestoreGlobalVector(dm, &X);
500: DMRestoreGlobalVector(dm, &F);
501: }
502: if (user.computeJacobian) {
503: Vec X;
504: Mat J;
505: MatStructure flag;
507: DMGetGlobalVector(dm, &X);
508: DMGetMatrix(dm, MATAIJ, &J);
509: if (user.batch) {
510: DMMeshSetLocalJacobian(dm, (PetscErrorCode (*)(DM, Vec, Mat, void*)) FormJacobianLocalBatch);
511: } else {
512: DMMeshSetLocalJacobian(dm, (PetscErrorCode (*)(DM, Vec, Mat, void*)) FormJacobianLocal);
513: }
514: SNESMeshFormJacobian(snes, X, &J, &J, &flag, &user);
515: MatDestroy(&J);
516: DMRestoreGlobalVector(dm, &X);
517: }
518: SNESDestroy(&snes);
519: DMDestroy(&dm);
520: PetscFinalize();
521: return 0;
522: }