Actual source code: ex12.c

  1: static char help[] = "Bratu nonlinear PDE in 2d.\n\
  2: We solve the  Bratu (SFI - solid fuel ignition) problem in a 2D rectangular\n\
  3: domain, using a parallel unstructured mesh (DMMESH) to discretize it.\n\
  4: The command line options include:\n\
  5:   -lambda <parameter>, where <parameter> indicates the problem's nonlinearity\n\
  6:      problem SFI:  <parameter> = Bratu parameter (0 <= lambda <= 6.81)\n\n";

  8: /*T
  9:    Concepts: SNES^parallel Bratu example
 10:    Concepts: DMMESH^using unstructured grids;
 11:    Processors: n
 12: T*/

 14: /* ------------------------------------------------------------------------

 16:     Solid Fuel Ignition (SFI) problem.  This problem is modeled by
 17:     the partial differential equation

 19:             -Laplacian u - lambda*exp(u) = f(x,y),  0 < x,y < 1,

 21:     with boundary conditions

 23:              u = 0  for  x = 0, x = 1, y = 0, y = 1.

 25:     A finite element approximation with the usual P_1 linear basis
 26:     is used to discretize the boundary value problem to obtain a nonlinear
 27:     system of equations.

 29:     Program usage:  mpiexec -n <procs> ex5 [-help] [all PETSc options]
 30:      e.g.,
 31:       ./ex5 -draw_pause -1
 32:       mpiexec -n 2 ./ex5 -log_summary

 34:   ------------------------------------------------------------------------- */

 36: /*
 37:    Include "petscdmmesh.h" so that we can use unstructured meshes (DMMESHs).
 38:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 39:    file automatically includes:
 40:      petscsys.h    - base PETSc routines
 41:      petscvec.h    - vectors               petscmat.h - matrices
 42:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 43:      petscviewer.h - viewers               petscpc.h  - preconditioners
 44: */
 45: #include <petscdmmesh.h>
 46: #include <petscsnes.h>

 48: typedef enum {RUN_FULL, RUN_TEST, RUN_MESH} RunType;
 49: typedef enum {NEUMANN, DIRICHLET} BCType;

 51: /*------------------------------------------------------------------------------
 52:   This code can be generated using config/PETSc/FEM.py

 54:     import PETSc.FEM
 55:     from FIAT.reference_element import default_simplex
 56:     from FIAT.lagrange import Lagrange

 58:     generator = PETSc.FEM.QuadratureGenerator()
 59:     generator.setup()
 60:     dim      = 2
 61:     order    = 1
 62:     elements = [Lagrange(default_simplex(dim), order))]
 63:     generator.run(elements, filename)
 64:  -----------------------------------------------------------------------------*/
 65: #include <stdlib.h>

 67: #define NUM_QUADRATURE_POINTS_0 1

 69: /* Quadrature points
 70:    - (x1,y1,x2,y2,...) */
 71: static PetscReal points_0[1] = {0.0};

 73: /* Quadrature weights
 74:    - (v1,v2,...) */
 75: static PetscReal weights_0[1] = {2.0};

 77: #define NUM_BASIS_FUNCTIONS_0 2

 79: /* Nodal basis function evaluations
 80:     - basis function is fastest varying, then point */
 81: static PetscReal Basis_0[2] = {
 82:   0.5,
 83:   0.5};

 85: /* Nodal basis function derivative evaluations,
 86:     - derivative direction fastest varying, then basis function, then point */
 87: static PetscReal BasisDerivatives_0[2] = {
 88:   -0.5,
 89:   0.5};

 91: #define NUM_DUAL_POINTS_0 2

 93: /* Dual points
 94:    - (x1,y1,x2,y2,...) */
 95: static PetscReal dualPoints_0[2] = {
 96:   -1.0,
 97:   1.0};



103: PetscReal IntegrateDualBasis_gen_0(const PetscReal *v0, const PetscReal *J, const int dualIndex, PetscReal (*func)(const PetscReal *coords))
104: {
105:   PetscReal refCoords[1];
106:   PetscReal coords[1];

108:   switch(dualIndex) {
109:     case 0:
110:       refCoords[0] = -1.0;
111:       break;
112:     case 1:
113:       refCoords[0] = 1.0;
114:       break;
115:     default:
116:       printf("dualIndex: %d\n", dualIndex);
117:       throw ALE::Exception("Bad dual index");
118:   }
119:   for(int d = 0; d < 1; d++)
120:   {
121:     coords[d] = v0[d];
122:     for(int e = 0; e < 1; e++)
123:     {
124:       coords[d] += J[d * 1 + e] * (refCoords[e] + 1.0);
125:     }
126:   }
127:   return (*func)(coords);
128: }



134: PetscReal IntegrateBdDualBasis_gen_0(const PetscReal *v0, const PetscReal *J, const int dualIndex, PetscReal (*func)(const PetscReal *coords))
135: {
136:   PetscReal refCoords[1];
137:   PetscReal coords[2];

139:   switch(dualIndex) {
140:     case 0:
141:       refCoords[0] = -1.0;
142:       break;
143:     case 1:
144:       refCoords[0] = 1.0;
145:       break;
146:     default:
147:       printf("dualIndex: %d\n", dualIndex);
148:       throw ALE::Exception("Bad dual index");
149:   }
150:   for(int d = 0; d < 2; d++)
151:   {
152:     coords[d] = v0[d];
153:     for(int e = 0; e < 1; e++)
154:     {
155:       coords[d] += J[d * 2 + e] * (refCoords[e] + 1.0);
156:     }
157:   }
158:   return (*func)(coords);
159: }

161: #define NUM_QUADRATURE_POINTS_1 1

163: /* Quadrature points
164:    - (x1,y1,x2,y2,...) */
165: static PetscReal points_1[2] = {
166:   -0.333333333333,
167:   -0.333333333333};

169: /* Quadrature weights
170:    - (v1,v2,...) */
171: static PetscReal weights_1[1] = {2.0};

173: #define NUM_BASIS_FUNCTIONS_1 3

175: /* Nodal basis function evaluations
176:     - basis function is fastest varying, then point */
177: static PetscReal Basis_1[3] = {
178:   0.333333333333,
179:   0.333333333333,
180:   0.333333333333};

182: /* Nodal basis function derivative evaluations,
183:     - derivative direction fastest varying, then basis function, then point */
184: static PetscReal BasisDerivatives_1[6] = {
185:   -0.5,
186:   -0.5,
187:   0.5,
188:   0.0,
189:   0.0,
190:   0.5};

192: #define NUM_DUAL_POINTS_1 3

194: /* Dual points
195:    - (x1,y1,x2,y2,...) */
196: static PetscReal dualPoints_1[6] = {
197:   -1.0,
198:   -1.0,
199:   1.0,
200:   -1.0,
201:   -1.0,
202:   1.0};



208: PetscReal IntegrateDualBasis_gen_1(const PetscReal *v0, const PetscReal *J, const int dualIndex, PetscReal (*func)(const PetscReal *coords))
209: {
210:   PetscReal refCoords[2];
211:   PetscReal coords[2];

213:   switch(dualIndex) {
214:     case 0:
215:       refCoords[0] = -1.0;
216:       refCoords[1] = -1.0;
217:       break;
218:     case 1:
219:       refCoords[0] = 1.0;
220:       refCoords[1] = -1.0;
221:       break;
222:     case 2:
223:       refCoords[0] = -1.0;
224:       refCoords[1] = 1.0;
225:       break;
226:     default:
227:       printf("dualIndex: %d\n", dualIndex);
228:       throw ALE::Exception("Bad dual index");
229:   }
230:   for(int d = 0; d < 2; d++)
231:   {
232:     coords[d] = v0[d];
233:     for(int e = 0; e < 2; e++)
234:     {
235:       coords[d] += J[d * 2 + e] * (refCoords[e] + 1.0);
236:     }
237:   }
238:   return (*func)(coords);
239: }



245: PetscReal IntegrateBdDualBasis_gen_1(const PetscReal *v0, const PetscReal *J, const int dualIndex, PetscReal (*func)(const PetscReal *coords))
246: {
247:   PetscReal refCoords[2];
248:   PetscReal coords[3];

250:   switch(dualIndex) {
251:     case 0:
252:       refCoords[0] = -1.0;
253:       refCoords[1] = -1.0;
254:       break;
255:     case 1:
256:       refCoords[0] = 1.0;
257:       refCoords[1] = -1.0;
258:       break;
259:     case 2:
260:       refCoords[0] = -1.0;
261:       refCoords[1] = 1.0;
262:       break;
263:     default:
264:       printf("dualIndex: %d\n", dualIndex);
265:       throw ALE::Exception("Bad dual index");
266:   }
267:   for(int d = 0; d < 3; d++)
268:   {
269:     coords[d] = v0[d];
270:     for(int e = 0; e < 2; e++)
271:     {
272:       coords[d] += J[d * 3 + e] * (refCoords[e] + 1.0);
273:     }
274:   }
275:   return (*func)(coords);
276: }


279: #define NUM_QUADRATURE_POINTS_2 1

281: /* Quadrature points
282:    - (x1,y1,x2,y2,...) */
283: static PetscReal points_2[3] = {
284:   -0.5,
285:   -0.5,
286:   -0.5};

288: /* Quadrature weights
289:    - (v1,v2,...) */
290: static PetscReal weights_2[1] = {1.33333333333};

292: #define NUM_BASIS_FUNCTIONS_2 4

294: /* Nodal basis function evaluations
295:     - basis function is fastest varying, then point */
296: static PetscReal Basis_2[4] = {
297:   0.25,
298:   0.25,
299:   0.25,
300:   0.25};

302: /* Nodal basis function derivative evaluations,
303:     - derivative direction fastest varying, then basis function, then point */
304: static PetscReal BasisDerivatives_2[12] = {
305:   -0.5,
306:   -0.5,
307:   -0.5,
308:   0.5,
309:   0.0,
310:   0.0,
311:   0.0,
312:   0.5,
313:   0.0,
314:   0.0,
315:   0.0,
316:   0.5};

318: #define NUM_DUAL_POINTS_2 4

320: /* Dual points
321:    - (x1,y1,x2,y2,...) */
322: static PetscReal dualPoints_2[12] = {
323:   -1.0,
324:   -1.0,
325:   -1.0,
326:   1.0,
327:   -1.0,
328:   -1.0,
329:   -1.0,
330:   1.0,
331:   -1.0,
332:   -1.0,
333:   -1.0,
334:   1.0};



340: PetscReal IntegrateDualBasis_gen_2(const PetscReal *v0, const PetscReal *J, const int dualIndex, PetscReal (*func)(const PetscReal *coords))
341: {
342:   PetscReal refCoords[3];
343:   PetscReal coords[3];

345:   switch(dualIndex) {
346:     case 0:
347:       refCoords[0] = -1.0;
348:       refCoords[1] = -1.0;
349:       refCoords[2] = -1.0;
350:       break;
351:     case 1:
352:       refCoords[0] = 1.0;
353:       refCoords[1] = -1.0;
354:       refCoords[2] = -1.0;
355:       break;
356:     case 2:
357:       refCoords[0] = -1.0;
358:       refCoords[1] = 1.0;
359:       refCoords[2] = -1.0;
360:       break;
361:     case 3:
362:       refCoords[0] = -1.0;
363:       refCoords[1] = -1.0;
364:       refCoords[2] = 1.0;
365:       break;
366:     default:
367:       printf("dualIndex: %d\n", dualIndex);
368:       throw ALE::Exception("Bad dual index");
369:   }
370:   for(int d = 0; d < 3; d++)
371:   {
372:     coords[d] = v0[d];
373:     for(int e = 0; e < 3; e++)
374:     {
375:       coords[d] += J[d * 3 + e] * (refCoords[e] + 1.0);
376:     }
377:   }
378:   return (*func)(coords);
379: }



385: PetscReal IntegrateBdDualBasis_gen_2(const PetscReal *v0, const PetscReal *J, const int dualIndex, PetscReal (*func)(const PetscReal *coords))
386: {
387:   PetscReal refCoords[3];
388:   PetscReal coords[4];

390:   switch(dualIndex) {
391:     case 0:
392:       refCoords[0] = -1.0;
393:       refCoords[1] = -1.0;
394:       refCoords[2] = -1.0;
395:       break;
396:     case 1:
397:       refCoords[0] = 1.0;
398:       refCoords[1] = -1.0;
399:       refCoords[2] = -1.0;
400:       break;
401:     case 2:
402:       refCoords[0] = -1.0;
403:       refCoords[1] = 1.0;
404:       refCoords[2] = -1.0;
405:       break;
406:     case 3:
407:       refCoords[0] = -1.0;
408:       refCoords[1] = -1.0;
409:       refCoords[2] = 1.0;
410:       break;
411:     default:
412:       printf("dualIndex: %d\n", dualIndex);
413:       throw ALE::Exception("Bad dual index");
414:   }
415:   for(int d = 0; d < 4; d++)
416:   {
417:     coords[d] = v0[d];
418:     for(int e = 0; e < 3; e++)
419:     {
420:       coords[d] += J[d * 4 + e] * (refCoords[e] + 1.0);
421:     }
422:   }
423:   return (*func)(coords);
424: }

426: /*------------------------------------------------------------------------------
427:   end of generated code
428:  -----------------------------------------------------------------------------*/

430: /*
431:    User-defined application context - contains data needed by the
432:    application-provided call-back routines, FormJacobianLocal() and
433:    FormFunctionLocal().
434: */
435: typedef struct {
436:   DM            dm;                /* The unstructured mesh data structure */
437:   PetscInt      debug;             /* The debugging level */
438:   PetscMPIInt   rank;              /* The process rank */
439:   PetscMPIInt   numProcs;          /* The number of processes */
440:   RunType       run;               /* The run type */
441:   PetscInt      dim;               /* The topological mesh dimension */
442:   PetscBool     interpolate;       /* Generate intermediate mesh elements */
443:   PetscReal     refinementLimit;   /* The largest allowable cell volume */
444:   char          partitioner[2048]; /* The graph partitioner */
445:   /* Element quadrature */
446:   PetscQuadrature q;
447:   /* Problem specific parameters */
448:   BCType        bcType;            /* The type of boundary conditions */
449:   PetscReal     lambda;            /* The Bratu problem parameter */
450:   PetscScalar (*rhsFunc)(const PetscReal []);   /* The rhs function f(x,y,z) */
451:   PetscScalar (*exactFunc)(const PetscReal []); /* The exact solution function u(x,y,z) */
452:   Vec           exactSol;          /* The discrete exact solution */
453:   Vec           error;             /* The discrete cell-wise error */
454: } AppCtx;

456: /*
457:    User-defined routines
458: */

463: PetscReal lambda = 0.0;
464: PetscScalar guess(const PetscReal coords[]) {
465:   PetscScalar scale = lambda/(lambda+1.0);
466:   return scale*(0.5 - fabs(coords[0]-0.5))*(0.5 - fabs(coords[1]-0.5));
467: }

469: PetscScalar zero(const PetscReal coords[]) {
470:   return 0.0;
471: }

473: PetscScalar constant(const PetscReal x[]) {
474:   return -4.0;
475: };

477: PetscScalar nonlinear_2d(const PetscReal x[]) {
478:   return -4.0 - lambda*PetscExpScalar(x[0]*x[0] + x[1]*x[1]);
479: };

481: PetscScalar linear_2d(const PetscReal x[]) {
482:   return -6.0*(x[0] - 0.5) - 6.0*(x[1] - 0.5);
483: };

485: PetscScalar quadratic_2d(const PetscReal x[]) {
486:   return x[0]*x[0] + x[1]*x[1];
487: };

489: PetscScalar cubic_2d(const PetscReal x[]) {
490:   return x[0]*x[0]*x[0] - 1.5*x[0]*x[0] + x[1]*x[1]*x[1] - 1.5*x[1]*x[1] + 0.5;
491: };

493: PetscScalar nonlinear_3d(const double x[]) {
494:   return -4.0 - lambda*PetscExpScalar((2.0/3.0)*(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]));
495: };

497: PetscScalar linear_3d(const double x[]) {
498:   return -6.0*(x[0] - 0.5) - 6.0*(x[1] - 0.5) - 6.0*(x[2] - 0.5);
499: };

501: PetscScalar quadratic_3d(const double x[]) {
502:   return (2.0/3.0)*(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
503: };

505: PetscScalar cubic_3d(const double x[]) {
506:   return x[0]*x[0]*x[0] - 1.5*x[0]*x[0] + x[1]*x[1]*x[1] - 1.5*x[1]*x[1] + x[2]*x[2]*x[2] - 1.5*x[2]*x[2] + 0.75;
507: };

511: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) {
512:   const char    *runTypes[3] = {"full", "test", "mesh"};
513:   const char    *bcTypes[2]  = {"neumann", "dirichlet"};
514:   PetscReal      bratu_lambda_max = 6.81, bratu_lambda_min = 0.0;
515:   PetscInt       run, bc;

519:   options->debug           = 0;
520:   options->run             = RUN_FULL;
521:   options->dim             = 2;
522:   options->interpolate     = PETSC_FALSE;
523:   options->refinementLimit = 0.0;
524:   options->bcType          = DIRICHLET;
525:   options->lambda          = 6.0;
526:   options->rhsFunc         = zero;

528:   MPI_Comm_size(comm, &options->numProcs);
529:   MPI_Comm_rank(comm, &options->rank);
530:   PetscOptionsBegin(comm, "", "Bratu Problem Options", "DMMESH");
531:   PetscOptionsInt("-debug", "The debugging level", "ex12.c", options->debug, &options->debug, PETSC_NULL);
532:   run = options->run;
533:   PetscOptionsEList("-run_type", "The run type", "ex12.c", runTypes, 3, runTypes[options->run], &run, PETSC_NULL);
534:   options->run = (RunType) run;
535:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex12.c", options->dim, &options->dim, PETSC_NULL);
536:   PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex12.c", options->interpolate, &options->interpolate, PETSC_NULL);
537:   PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex12.c", options->refinementLimit, &options->refinementLimit, PETSC_NULL);
538:   PetscStrcpy(options->partitioner, "chaco");
539:   PetscOptionsString("-partitioner", "The graph partitioner", "pflotran.cxx", options->partitioner, options->partitioner, 2048, PETSC_NULL);
540:   bc = options->bcType;
541:   PetscOptionsEList("-bc_type","Type of boundary condition","ex12.c",bcTypes,2,bcTypes[options->bcType],&bc,PETSC_NULL);
542:   options->bcType = (BCType) bc;
543:   PetscOptionsReal("-lambda", "The parameter controlling nonlinearity", "ex12.c", options->lambda, &options->lambda, PETSC_NULL);
544:   if (options->lambda >= bratu_lambda_max || options->lambda < bratu_lambda_min) {
545:     SETERRQ3(PETSC_COMM_WORLD, 1, "Lambda, %g, is out of range, [%g, %g)", options->lambda, bratu_lambda_min, bratu_lambda_max);
546:   }
547:   PetscOptionsEnd();
548:   lambda = options->lambda;
549:   return(0);
550: };

554: PetscErrorCode SetupQuadrature(AppCtx *user) {
556:   switch(user->dim) {
557:   case 1:
558:     user->q.numQuadPoints = NUM_QUADRATURE_POINTS_0;
559:     user->q.quadPoints    = points_0;
560:     user->q.quadWeights   = weights_0;
561:     user->q.numBasisFuncs = NUM_BASIS_FUNCTIONS_0;
562:     user->q.basis         = Basis_0;
563:     user->q.basisDer      = BasisDerivatives_0;
564:     break;
565:   case 2:
566:     user->q.numQuadPoints = NUM_QUADRATURE_POINTS_1;
567:     user->q.quadPoints    = points_1;
568:     user->q.quadWeights   = weights_1;
569:     user->q.numBasisFuncs = NUM_BASIS_FUNCTIONS_1;
570:     user->q.basis         = Basis_1;
571:     user->q.basisDer      = BasisDerivatives_1;
572:     break;
573:   case 3:
574:     user->q.numQuadPoints = NUM_QUADRATURE_POINTS_2;
575:     user->q.quadPoints    = points_2;
576:     user->q.quadWeights   = weights_2;
577:     user->q.numBasisFuncs = NUM_BASIS_FUNCTIONS_2;
578:     user->q.basis         = Basis_2;
579:     user->q.basisDer      = BasisDerivatives_2;
580:     break;
581:   default:
582:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
583:   }
584:   return(0);
585: }

589: PetscErrorCode SetupSection(AppCtx *user) {
590:   PetscSection   section;
591:   /* These can be generated using config/PETSc/FEM.py */
592:   PetscInt       numDof_0[2] = {1, 0};
593:   PetscInt       numDof_1[3] = {1, 0, 0};
594:   PetscInt       numDof_2[4] = {1, 0, 0, 0};
595:   PetscInt      *numDof;
596:   const char    *bcLabel = PETSC_NULL;

600:   switch(user->dim) {
601:   case 1:
602:     numDof = numDof_0;
603:     break;
604:   case 2:
605:     numDof = numDof_1;
606:     break;
607:   case 3:
608:     numDof = numDof_2;
609:     break;
610:   default:
611:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid spatial dimension %d", user->dim);
612:   }
613:   if (user->bcType == DIRICHLET) {
614:     bcLabel = "marker";
615:   }
616:   DMMeshCreateSection(user->dm, user->dim, numDof, bcLabel, 1, &section);
617:   DMMeshSetSection(user->dm, "default", section);
618:   return(0);
619: }

623: PetscErrorCode SetupExactSolution(AppCtx *user) {
625:   switch(user->dim) {
626:   case 2:
627:     if (user->bcType == DIRICHLET) {
628:       if (user->lambda > 0.0) {
629:         user->rhsFunc   = nonlinear_2d;
630:         user->exactFunc = quadratic_2d;
631:       } else {
632:         user->rhsFunc   = constant;
633:         user->exactFunc = quadratic_2d;
634:       }
635:     } else {
636:       user->rhsFunc   = linear_2d;
637:       user->exactFunc = cubic_2d;
638:     }
639:     break;
640:   case 3:
641:     if (user->bcType == DIRICHLET) {
642:       if (user->lambda > 0.0) {
643:         user->rhsFunc   = nonlinear_3d;
644:         user->exactFunc = quadratic_3d;
645:       } else {
646:         user->rhsFunc   = constant;
647:         user->exactFunc = quadratic_3d;
648:       }
649:     } else {
650:       user->rhsFunc   = linear_3d;
651:       user->exactFunc = cubic_3d;
652:     }
653:     break;
654:   default:
655:     SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d", user->dim);
656:   }
657:   return(0);
658: }

662: PetscErrorCode ComputeError(Vec X, PetscReal *error, AppCtx *user) {
663:   PetscScalar    (*exactFunc)(const PetscReal []) = user->exactFunc;
664:   const PetscInt   debug         = user->debug;
665:   const PetscInt   dim           = user->dim;
666:   const PetscInt   numQuadPoints = user->q.numQuadPoints;
667:   const PetscReal *quadPoints    = user->q.quadPoints;
668:   const PetscReal *quadWeights   = user->q.quadWeights;
669:   const PetscInt   numBasisFuncs = user->q.numBasisFuncs;
670:   const PetscReal *basis         = user->q.basis;
671:   Vec              localX;
672:   PetscReal       *coords, *v0, *J, *invJ, detJ;
673:   PetscReal        localError;
674:   PetscInt         cStart, cEnd;
675:   PetscErrorCode   ierr;

678:   DMGetLocalVector(user->dm, &localX);
679:   DMGlobalToLocalBegin(user->dm, X, INSERT_VALUES, localX);
680:   DMGlobalToLocalEnd(user->dm, X, INSERT_VALUES, localX);
681:   PetscMalloc4(dim,PetscReal,&coords,dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);
682:   DMMeshGetHeightStratum(user->dm, 0, &cStart, &cEnd);
683:   for(PetscInt c = cStart; c < cEnd; ++c) {
684:     const PetscScalar *x;
685:     PetscReal          elemError = 0.0;

687:     DMMeshComputeCellGeometry(user->dm, c, v0, J, invJ, &detJ);
688:     if (detJ <= 0.0) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);}
689:     DMMeshVecGetClosure(user->dm, localX, c, &x);
690:     if (debug) {DMMeshPrintCellVector(c, "Solution", numBasisFuncs, x);}
691:     for(int q = 0; q < numQuadPoints; ++q) {
692:       for(int d = 0; d < dim; d++) {
693:         coords[d] = v0[d];
694:         for(int e = 0; e < dim; e++) {
695:           coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
696:         }
697:       }
698:       const PetscScalar funcVal     = (*exactFunc)(coords);
699:       PetscReal         interpolant = 0.0;
700:       for(int f = 0; f < numBasisFuncs; ++f) {
701:         interpolant += x[f]*basis[q*numBasisFuncs+f];
702:       }
703:       elemError += PetscSqr(interpolant - funcVal)*quadWeights[q]*detJ;
704:     }
705:     if (debug) {PetscPrintf(PETSC_COMM_SELF, "  elem %d error %g\n", c, elemError);}
706:     localError += elemError;
707:   }
708:   PetscFree4(coords,v0,J,invJ);
709:   DMRestoreLocalVector(user->dm, &localX);
710:   MPI_Allreduce(&localError, error, 1, MPIU_REAL, MPI_SUM, PETSC_COMM_WORLD);
711:   *error = PetscSqrtReal(*error);
712:   return(0);
713: }

717: int main(int argc, char **argv)
718: {
719:   SNES           snes;                 /* nonlinear solver */
720:   Vec            u,r;                  /* solution, residual vectors */
721:   Mat            A,J;                  /* Jacobian matrix */
722:   AppCtx         user;                 /* user-defined work context */
723:   PetscInt       its;                  /* iterations for convergence */
724:   PetscReal      error;                /* L_2 error in the solution */

727:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
728:      Initialize program
729:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
730:   PetscInitialize(&argc, &argv, PETSC_NULL, help);

732:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
733:      Initialize problem parameters
734:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
735:   ProcessOptions(PETSC_COMM_WORLD, &user);

737:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
738:      Create nonlinear solver context
739:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
740:   SNESCreate(PETSC_COMM_WORLD, &snes);

742:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
743:      Create unstructured mesh (DMMESH) to manage parallel grid and vectors
744:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
745:   DMMeshCreateBoxMesh(PETSC_COMM_WORLD, user.dim, user.interpolate, &user.dm);
746:   {
747:     DM refinedMesh     = PETSC_NULL;
748:     DM distributedMesh = PETSC_NULL;

750:     /* Refine mesh using a volume constraint */
751:     DMMeshRefine(user.dm, user.refinementLimit, user.interpolate, &refinedMesh);
752:     if (refinedMesh) {
753:       DMDestroy(&user.dm);
754:       user.dm = refinedMesh;
755:     }
756: #if 0
757:     {
758:       ALE::Obj<PETSC_MESH_TYPE> oldMesh;
759:       DMMeshGetMesh(user.dm, oldMesh);
760:       oldMesh->setDebug(1);
761:     }
762: #endif
763:     /* Distribute mesh over processes */
764:     DMMeshDistribute(user.dm, user.partitioner, &distributedMesh);
765:     if (distributedMesh) {
766:       DMDestroy(&user.dm);
767:       user.dm = distributedMesh;
768:     }
769:     /* Mark boundary cells for higher order element calculations */
770:     if (user.bcType == DIRICHLET) {
771:       DMMeshMarkBoundaryCells(user.dm, "marker", 1, 2);
772:     }
773:   }
774:   DMSetFromOptions(user.dm);
775:   SNESSetDM(snes, user.dm);

777:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
778:      Setup dof layout.

780:      For a DMDA, this is automatic given the number of dof at each vertex.
781:      However, for a DMMesh, we need to specify this.
782:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
783:   SetupExactSolution(&user);
784:   SetupQuadrature(&user);
785:   SetupSection(&user);
786:   if (user.bcType == NEUMANN) {
787:     /* With Neumann conditions, we tell DMMG that constants are in the null space of the operator
788:          Should have a nice one like DMMG that sets it for all MG PCs */
789:     KSP          ksp;
790:     MatNullSpace nullsp;

792:     SNESGetKSP(snes, &ksp);
793:     MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, PETSC_NULL, &nullsp);
794:     KSPSetNullSpace(ksp, nullsp);
795:     MatNullSpaceDestroy(&nullsp);
796:   }

798:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
799:      Extract global vectors from DM; then duplicate for remaining
800:      vectors that are the same types
801:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
802:   DMCreateGlobalVector(user.dm, &u);
803:   VecDuplicate(u, &r);

805:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
806:      Create matrix data structure; set Jacobian evaluation routine

808:      Set Jacobian matrix data structure and default Jacobian evaluation
809:      routine. User can override with:
810:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
811:                 (unless user explicitly sets preconditioner)
812:      -snes_mf_operator : form preconditioning matrix as set by the user,
813:                          but use matrix-free approx for Jacobian-vector
814:                          products within Newton-Krylov method
815:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
816:   /* J can be type of MATAIJ, MATBAIJ or MATSBAIJ */
817:   DMGetMatrix(user.dm, MATAIJ, &J);
818:   A    = J;
819:   SNESSetJacobian(snes, A, J, SNESMeshFormJacobian, &user);

821:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
822:      Set local function evaluation routine
823:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
824:   DMMeshSetLocalFunction(user.dm, (DMMeshLocalFunction1) FormFunctionLocal);
825:   DMMeshSetLocalJacobian(user.dm, (DMMeshLocalJacobian1) FormJacobianLocal);
826:   SNESSetFunction(snes, r, SNESMeshFormFunction, &user);

828:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
829:      Customize nonlinear solver; set runtime options
830:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
831:   SNESSetFromOptions(snes);

833:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
834:      Setup boundary conditions
835:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
836:   FormInitialGuess(u, user.exactFunc, INSERT_ALL_VALUES, &user);
837:   if (user.run == RUN_FULL) {
838:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
839:      Evaluate initial guess
840:      Note: The user should initialize the vector u, with the initial guess
841:      for the nonlinear solver prior to calling SNESSolve().  In particular,
842:      to employ an initial guess of zero, the user should explicitly set
843:      this vector to zero by calling VecSet().
844:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
845:     FormInitialGuess(u, guess, INSERT_VALUES, &user);
846:     if (user.debug) {
847:       PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
848:       VecView(u, PETSC_VIEWER_STDOUT_WORLD);
849:     }
850:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
851:      Solve nonlinear system
852:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
853:     SNESSolve(snes, PETSC_NULL, u);
854:     SNESGetIterationNumber(snes, &its);
855:     PetscPrintf(PETSC_COMM_WORLD, "Number of Newton iterations = %D\n", its);
856:     ComputeError(u, &error, &user);
857:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);
858:     PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
859:     VecView(u, PETSC_VIEWER_STDOUT_WORLD);
860:   } else {
861:     PetscReal res;

863:     /* Check discretization error */
864:     PetscPrintf(PETSC_COMM_WORLD, "Initial guess\n");
865:     VecView(u, PETSC_VIEWER_STDOUT_WORLD);
866:     ComputeError(u, &error, &user);
867:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Error: %g\n", error);
868:     /* Check residual */
869:     SNESMeshFormFunction(snes, u, r, &user);
870:     VecNorm(r, NORM_2, &res);
871:     PetscPrintf(PETSC_COMM_WORLD, "L_2 Residual: %g\n", res);
872:   }

874:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
875:      Output results
876:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
877:   if (0) {
878:     PetscViewer viewer;

880:     PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
881:     /*PetscViewerSetType(viewer, PETSCVIEWERDRAW);
882:       PetscViewerDrawSetInfo(viewer, PETSC_NULL, "Solution", PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE); */
883:     PetscViewerSetType(viewer, PETSCVIEWERASCII);
884:     PetscViewerFileSetName(viewer, "ex12_sol.vtk");
885:     PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_VTK);
886:     DMView(user.dm, viewer);
887:     VecView(u, viewer);
888:     PetscViewerDestroy(&viewer);
889:   }

891:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
892:      Free work space.  All PETSc objects should be destroyed when they
893:      are no longer needed.
894:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
895:   if (A != J) {
896:     MatDestroy(&A);
897:   }
898:   MatDestroy(&J);
899:   VecDestroy(&u);
900:   VecDestroy(&r);
901:   SNESDestroy(&snes);
902:   DMDestroy(&user.dm);
903:   PetscFinalize();
904:   return 0;
905: }

909: /*
910:   FormInitialGuess - Forms initial approximation.

912:   Input Parameters:
913: + user - user-defined application context
914: - guessFunc - The coordinate function to use for the guess

916:   Output Parameter:
917: . X - vector
918: */
919: PetscErrorCode FormInitialGuess(Vec X, PetscScalar (*guessFunc)(const PetscReal []), InsertMode mode, AppCtx *user)
920: {
921:   Vec            localX, coordinates;
922:   PetscSection   section, cSection;
923:   PetscInt       vStart, vEnd;

927:   DMGetLocalVector(user->dm, &localX);
928:   DMMeshGetDepthStratum(user->dm, 0, &vStart, &vEnd);
929:   DMMeshGetDefaultSection(user->dm, &section);
930:   DMMeshGetCoordinateSection(user->dm, &cSection);
931:   DMMeshGetCoordinateVec(user->dm, &coordinates);
932:   for(PetscInt v = vStart; v < vEnd; ++v) {
933:     PetscScalar  values[1];
934:     PetscScalar *coords;

936:     VecGetValuesSection(coordinates, cSection, v, &coords);
937:     values[0] = (*guessFunc)(coords);
938:     VecSetValuesSection(localX, section, v, values, mode);
939:   }
940:   VecDestroy(&coordinates);
941:   PetscSectionDestroy(&section);
942:   PetscSectionDestroy(&cSection);

944:   DMLocalToGlobalBegin(user->dm, localX, INSERT_VALUES, X);
945:   DMLocalToGlobalEnd(user->dm, localX, INSERT_VALUES, X);
946:   DMRestoreLocalVector(user->dm, &localX);
947: #if 0
948:   /* This is necessary for higher order elements */
949:   MeshGetSectionReal(mesh, "exactSolution", &this->_options.exactSol.section);
950:   const Obj<PETSC_MESH_TYPE::real_section_type>& s = this->_mesh->getRealSection("exactSolution");
951:   this->_mesh->setupField(s);
952:   const Obj<PETSC_MESH_TYPE::label_sequence>&     cells       = this->_mesh->heightStratum(0);
953:   const Obj<PETSC_MESH_TYPE::real_section_type>&  coordinates = this->_mesh->getRealSection("coordinates");
954:   const int                                       localDof    = this->_mesh->sizeWithBC(s, *cells->begin());
955:   PETSC_MESH_TYPE::real_section_type::value_type *values      = new PETSC_MESH_TYPE::real_section_type::value_type[localDof];
956:   PetscReal                                      *v0          = new PetscReal[dim()];
957:   PetscReal                                      *J           = new PetscReal[dim()*dim()];
958:   PetscReal                                       detJ;
959:   ALE::ISieveVisitor::PointRetriever<PETSC_MESH_TYPE::sieve_type> pV((int) pow(this->_mesh->getSieve()->getMaxConeSize(), this->_mesh->depth())+1, true);

961:   for(PETSC_MESH_TYPE::label_sequence::iterator c_iter = cells->begin(); c_iter != cells->end(); ++c_iter) {
962:     ALE::ISieveTraversal<PETSC_MESH_TYPE::sieve_type>::orientedClosure(*this->_mesh->getSieve(), *c_iter, pV);
963:     const PETSC_MESH_TYPE::point_type *oPoints = pV.getPoints();
964:     const int                          oSize   = pV.getSize();
965:     int                                v       = 0;

967:     this->_mesh->computeElementGeometry(coordinates, *c_iter, v0, J, PETSC_NULL, detJ);
968:     for(int cl = 0; cl < oSize; ++cl) {
969:       const int pointDim = s->getFiberDimension(oPoints[cl]);

971:       if (pointDim) {
972:         for(int d = 0; d < pointDim; ++d, ++v) {
973:           values[v] = (*this->_options.integrate)(v0, J, v, initFunc);
974:         }
975:       }
976:     }
977:     this->_mesh->update(s, *c_iter, values);
978:     pV.clear();
979:   }
980:   delete [] values;
981:   delete [] v0;
982:   delete [] J;
983: #endif
984:   return(0);
985: }

989: /*
990:    FormFunctionLocal - Evaluates nonlinear function, F(x).
991: */
992: PetscErrorCode FormFunctionLocal(DM dm, Vec X, Vec F, AppCtx *user)
993: {
994:   PetscScalar    (*rhsFunc)(const PetscReal []) = user->rhsFunc;
995:   const PetscInt   debug         = user->debug;
996:   const PetscInt   dim           = user->dim;
997:   const PetscInt   numQuadPoints = user->q.numQuadPoints;
998:   const PetscReal *quadPoints    = user->q.quadPoints;
999:   const PetscReal *quadWeights   = user->q.quadWeights;
1000:   const PetscInt   numBasisFuncs = user->q.numBasisFuncs;
1001:   const PetscReal *basis         = user->q.basis;
1002:   const PetscReal *basisDer      = user->q.basisDer;
1003:   const PetscReal  lambda        = user->lambda;
1004:   PetscReal       *coords, *v0, *J, *invJ, detJ;
1005:   PetscScalar     *realSpaceDer, *fieldGrad, *elemVec;
1006:   PetscInt         cStart, cEnd;
1007:   PetscErrorCode   ierr;

1010:   VecSet(F, 0.0);
1011:   PetscMalloc3(dim,PetscScalar,&realSpaceDer,dim,PetscScalar,&fieldGrad,numBasisFuncs,PetscScalar,&elemVec);
1012:   PetscMalloc4(dim,PetscReal,&coords,dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);
1013:   DMMeshGetHeightStratum(user->dm, 0, &cStart, &cEnd);
1014:   for(PetscInt c = cStart; c < cEnd; ++c) {
1015:     const PetscScalar *x;

1017:     PetscMemzero(elemVec, numBasisFuncs * sizeof(PetscScalar));
1018:     DMMeshComputeCellGeometry(user->dm, c, v0, J, invJ, &detJ);
1019:     if (detJ <= 0.0) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);}
1020:     DMMeshVecGetClosure(user->dm, X, c, &x);
1021:     if (debug) {DMMeshPrintCellVector(c, "Solution", numBasisFuncs, x);}

1023:     for(int q = 0; q < numQuadPoints; ++q) {
1024:       PetscScalar fieldVal = 0.0;

1026:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
1027:       for(int d = 0; d < dim; ++d) {
1028:         fieldGrad[d] = 0.0;
1029:         coords[d] = v0[d];
1030:         for(int e = 0; e < dim; ++e) {
1031:           coords[d] += J[d*dim+e]*(quadPoints[q*dim+e] + 1.0);
1032:         }
1033:         if (debug) {PetscPrintf(PETSC_COMM_SELF, "    coords[%d] %g\n", d, coords[d]);}
1034:       }
1035:       for(int f = 0; f < numBasisFuncs; ++f) {
1036:         fieldVal += x[f]*basis[q*numBasisFuncs+f];

1038:         for(int d = 0; d < dim; ++d) {
1039:           realSpaceDer[d] = 0.0;
1040:           for(int e = 0; e < dim; ++e) {
1041:             realSpaceDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+f)*dim+e];
1042:           }
1043:           fieldGrad[d] += realSpaceDer[d]*x[f];
1044:         }
1045:       }
1046:       if (debug) {
1047:         for(int d = 0; d < dim; ++d) {
1048:           PetscPrintf(PETSC_COMM_SELF, "    fieldGrad[%d] %g\n", d, fieldGrad[d]);
1049:         }
1050:       }
1051:       const PetscScalar funcVal = (*rhsFunc)(coords);
1052:       for(int f = 0; f < numBasisFuncs; ++f) {
1053:         /* Constant term: -f(x) */
1054:         elemVec[f] -= basis[q*numBasisFuncs+f]*funcVal*quadWeights[q]*detJ;
1055:         /* Linear term: -\Delta u */
1056:         PetscScalar product = 0.0;
1057:         for(int d = 0; d < dim; ++d) {
1058:           realSpaceDer[d] = 0.0;
1059:           for(int e = 0; e < dim; ++e) {
1060:             realSpaceDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+f)*dim+e];
1061:           }
1062:           product += realSpaceDer[d]*fieldGrad[d];
1063:         }
1064:         elemVec[f] += product*quadWeights[q]*detJ;
1065:         /* Nonlinear term: -\lambda e^{u} */
1066:         elemVec[f] -= basis[q*numBasisFuncs+f]*lambda*PetscExpScalar(fieldVal)*quadWeights[q]*detJ;
1067:       }
1068:     }
1069:     if (debug) {DMMeshPrintCellVector(c, "Residual", numBasisFuncs, elemVec);}
1070:     DMMeshVecSetClosure(user->dm, F, c, elemVec, ADD_VALUES);
1071:   }
1072:   PetscLogFlops((cEnd-cStart)*numQuadPoints*numBasisFuncs*(dim*(dim*5+4)+14));
1073:   PetscFree3(realSpaceDer,fieldGrad,elemVec);
1074:   PetscFree4(coords,v0,J,invJ);

1076:   PetscPrintf(PETSC_COMM_WORLD, "Residual:\n");
1077:   for(int p = 0; p < user->numProcs; ++p) {
1078:     if (p == user->rank) {VecView(F, PETSC_VIEWER_STDOUT_SELF);}
1079:     PetscBarrier((PetscObject) user->dm);
1080:   }
1081:   return(0);
1082: }

1086: /*
1087:    FormJacobianLocal - Evaluates Jacobian matrix.
1088: */
1089: PetscErrorCode FormJacobianLocal(DM dm, Vec X, Mat Jac, AppCtx *user)
1090: {
1091:   const PetscInt   debug         = user->debug;
1092:   const PetscInt   dim           = user->dim;
1093:   const PetscInt   numQuadPoints = user->q.numQuadPoints;
1094:   const PetscReal *quadWeights   = user->q.quadWeights;
1095:   const PetscInt   numBasisFuncs = user->q.numBasisFuncs;
1096:   const PetscReal *basis         = user->q.basis;
1097:   const PetscReal *basisDer      = user->q.basisDer;
1098:   const PetscReal  lambda        = user->lambda;
1099:   PetscReal       *v0, *J, *invJ, detJ;
1100:   PetscScalar     *realSpaceTestDer, *realSpaceBasisDer, *elemMat;
1101:   PetscInt         cStart, cEnd;
1102:   PetscErrorCode   ierr;

1105:   MatZeroEntries(Jac);
1106:   PetscMalloc3(dim,PetscScalar,&realSpaceTestDer,dim,PetscScalar,&realSpaceBasisDer,numBasisFuncs*numBasisFuncs,PetscScalar,&elemMat);
1107:   PetscMalloc3(dim,PetscReal,&v0,dim*dim,PetscReal,&J,dim*dim,PetscReal,&invJ);
1108:   DMMeshGetHeightStratum(user->dm, 0, &cStart, &cEnd);
1109:   for(PetscInt c = cStart; c < cEnd; ++c) {
1110:     const PetscScalar *x;

1112:     PetscMemzero(elemMat, numBasisFuncs*numBasisFuncs * sizeof(PetscScalar));
1113:     DMMeshComputeCellGeometry(user->dm, c, v0, J, invJ, &detJ);
1114:     if (detJ <= 0.0) {SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ, c);}
1115:     DMMeshVecGetClosure(user->dm, X, c, &x);

1117:     for(int q = 0; q < numQuadPoints; ++q) {
1118:       PetscScalar fieldVal = 0.0;

1120:       for(int f = 0; f < numBasisFuncs; ++f) {
1121:         fieldVal += x[f]*basis[q*numBasisFuncs+f];
1122:       }
1123:       for(int f = 0; f < numBasisFuncs; ++f) {
1124:         for(int d = 0; d < dim; ++d) {
1125:           realSpaceTestDer[d] = 0.0;
1126:           for(int e = 0; e < dim; ++e) {
1127:             realSpaceTestDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+f)*dim+e];
1128:           }
1129:         }
1130:         for(int g = 0; g < numBasisFuncs; ++g) {
1131:           for(int d = 0; d < dim; ++d) {
1132:             realSpaceBasisDer[d] = 0.0;
1133:             for(int e = 0; e < dim; ++e) {
1134:               realSpaceBasisDer[d] += invJ[e*dim+d]*basisDer[(q*numBasisFuncs+g)*dim+e];
1135:             }
1136:           }
1137:           /* Linear term: -\Delta u */
1138:           PetscScalar product = 0.0;
1139:           for(int d = 0; d < dim; ++d) product += realSpaceTestDer[d]*realSpaceBasisDer[d];
1140:           elemMat[f*numBasisFuncs+g] += product*quadWeights[q]*detJ;
1141:           /* Nonlinear term: -\lambda e^{u} */
1142:           elemMat[f*numBasisFuncs+g] -= basis[q*numBasisFuncs+f]*basis[q*numBasisFuncs+g]*lambda*PetscExpScalar(fieldVal)*quadWeights[q]*detJ;
1143:         }
1144:       }
1145:     }
1146:     if (debug) {DMMeshPrintCellMatrix(c, "Jacobian", numBasisFuncs, numBasisFuncs, elemMat);}
1147:     DMMeshMatSetClosure(user->dm, Jac, c, elemMat, ADD_VALUES);
1148:   }
1149:   PetscLogFlops((cEnd-cStart)*numQuadPoints*numBasisFuncs*(dim*(dim*5+4)+14));
1150:   PetscFree3(realSpaceTestDer,realSpaceBasisDer,elemMat);
1151:   PetscFree3(v0,J,invJ);

1153:   /* Assemble matrix, using the 2-step process:
1154:        MatAssemblyBegin(), MatAssemblyEnd(). */
1155:   MatAssemblyBegin(Jac, MAT_FINAL_ASSEMBLY);
1156:   MatAssemblyEnd(Jac, MAT_FINAL_ASSEMBLY);
1157:   /* Tell the matrix we will never add a new nonzero location to the
1158:      matrix. If we do, it will generate an error. */
1159:   MatSetOption(Jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
1160:   return(0);
1161: }