Actual source code: ex37.c

  1: /*T
  2:    Concepts: KSP^solving a system of linear equations
  3:    Concepts: KSP^Laplacian, 2d
  4:    Processors: n
  5: T*/

  7: /*
  8: Added at the request of Bob Eisenberg.

 10: Inhomogeneous Laplacian in a 2D ion channel. Modeled by the partial differential equation

 12:    div \epsilon grad u = f,  on \Omega,

 14: with forcing function

 16:    f = \sum_i q_i \delta(x - x_i)

 18: with Dirichlet boundary conditions

 20:    u = 0 for x = -L
 21:    u = V for x =  L

 23: and Neumman boundary conditions

 25:    \hat n \cdot \grad u = 0 for y = -W, W

 27: This uses multigrid to solve the linear system on a 2D radially-symmetric channel boundary:

 29:              28                         29      35         43
 30:               V                          V      V          V 
 31:     2----------------------------------3-----4----12--------------------------------13
 32:     |                                  |     |     |                                 |
 33:     |                                  |     |     |                                 |
 34:     |                               34>|     |     | <36                             |
 35:     |                                  |     |     |                                 |
 36:  27>|                                  |  30>|     |                                 |
 37:     |                                  8     |     11                            42> |
 38:     |                               33>\     |     / <37                             |
 39:     |                                   7 31 | 39 10                                 |
 40:     |                                32> \ V | V / <38                               |
 41:     |                                     6--5--9      41                            |
 42:     |                                        |<40      V                             |
 43:     1----------------------------------------O--------------------------------------14
 44:     |          ^                             |<50                                    |
 45:     |         57                         25-20--19                                   |
 46:     |                               56>  / ^ | ^ \ <48                               |
 47:     |                                   24 51| 49 18                                 |
 48:     |                              55> /     |     \ <47                             |
 49:     | <58                              23    |    17                             44> |
 50:     |                                  | 52> |     |                                 |
 51:     |                                  |     |     |                                 |
 52:     |                              54> |     |(XX) |<46                              |
 53:     |       59                         | 53  | 60  |        45                       |
 54:     |        V                         |  V  |  V  |        V                        |
 55:     26(X)-----------------------------22----21-----16-------------------------------15

 57:     (X) denotes the last vertex, (XX) denotes the last edge
 58: */

 60: static char help[] = "Solves 2D inhomogeneous Laplacian using multigrid in an ion channel.\n\n";

 62:  #include petscmesh.h
 63:  #include petscksp.h
 64:  #include petscmg.h
 65:  #include petscdmmg.h

 67: PetscErrorCode MeshView_Sieve_Newer(ALE::Obj<ALE::Mesh>, PetscViewer);
 68: PetscErrorCode CreateMeshBoundary(ALE::Obj<ALE::Mesh>);
 69: PetscErrorCode updateOperator(Mat, ALE::Obj<ALE::Mesh::field_type>, const ALE::Mesh::point_type&, PetscScalar [], InsertMode);

 76: double refineLimit(const double [], void *);

 78: typedef enum {DIRICHLET, NEUMANN} BCType;

 80: typedef struct {
 81:   ALE::Obj<ALE::Mesh::field_type> epsilon;
 82:   PetscScalar voltage;
 83:   VecScatter  injection;
 84:   PetscReal   refinementLimit;
 85:   PetscReal   refinementExp;

 87:   PetscInt    numCharges;
 88:   PetscReal  *charge;
 89:   PetscReal  *chargeLocation;
 90: } UserContext;

 92: PetscInt debug;

 96: int main(int argc,char **argv)
 97: {
 98:   MPI_Comm       comm;
 99:   DMMG          *dmmg;
100:   UserContext    user;
101:   PetscViewer    viewer;
102:   PetscReal      norm;
103:   PetscInt       dim, l, meshDebug;
104:   PetscTruth     viewEnergy;

107:   PetscInitialize(&argc,&argv,(char *)0,help);
108:   comm = PETSC_COMM_WORLD;
109:   PetscOptionsBegin(comm, "", "Options for the inhomogeneous Poisson equation", "DMMG");
110:     debug = 0;
111:     PetscOptionsInt("-debug", "The debugging flag", "ex37.c", 0, &debug, PETSC_NULL);
112:     meshDebug = 0;
113:     PetscOptionsInt("-mesh_debug", "The mesh debugging flag", "ex37.c", 0, &meshDebug, PETSC_NULL);
114:     dim  = 2;
115:     PetscOptionsInt("-dim", "The mesh dimension", "ex37.c", 2, &dim, PETSC_NULL);
116:     user.refinementLimit = 0.0;
117:     PetscOptionsReal("-refinement_limit", "The area of the largest triangle in the mesh", "ex37.c", 1.0, &user.refinementLimit, PETSC_NULL);
118:     user.refinementExp = 0.0;
119:     PetscOptionsReal("-refinement_exp", "The exponent of the radius for refinement", "ex37.c", 1.0, &user.refinementExp, PETSC_NULL);
120:     user.voltage = 1.0;
121:     PetscOptionsScalar("-voltage", "The voltage of the clamp", "ex37.c", 1.0, &user.voltage, PETSC_NULL);
122:     viewEnergy = PETSC_FALSE;
123:     PetscOptionsTruth("-view_energy", "View the energy density as a field", "ex37.c", PETSC_FALSE, &viewEnergy, PETSC_NULL);
124:   PetscOptionsEnd();

126:   ALE::Obj<ALE::Mesh> meshBoundary = ALE::Mesh(comm, dim-1, meshDebug);
127:   ALE::Obj<ALE::Mesh> mesh;

129:   try {
130:     ALE::LogStage stage = ALE::LogStageRegister("MeshCreation");
131:     ALE::LogStagePush(stage);
132:     PetscPrintf(comm, "Generating mesh\n");
133:     CreateMeshBoundary(meshBoundary);
134:     mesh = ALE::Generator::generate(meshBoundary);
135:     ALE::Obj<ALE::Mesh::sieve_type> topology = mesh->getTopology();
136:     PetscPrintf(comm, "  Read %d elements\n", topology->heightStratum(0)->size());
137:     PetscPrintf(comm, "  Read %d vertices\n", topology->depthStratum(0)->size());
138:     ALE::LogStagePop(stage);
139:     mesh->getBoundary()->view("Initial Mesh Boundary");

141:     stage = ALE::LogStageRegister("MeshDistribution");
142:     ALE::LogStagePush(stage);
143:     PetscPrintf(comm, "Distributing mesh\n");
144:     mesh = mesh->distribute();
145:     ALE::LogStagePop(stage);
146:     mesh->getBoundary()->view("Distributed Mesh Boundary");

148:     if (user.refinementLimit > 0.0) {
149:       stage = ALE::LogStageRegister("MeshRefine");
150:       ALE::LogStagePush(stage);
151:       PetscPrintf(comm, "Refining mesh\n");
152:       if (user.refinementExp == 0.0) {
153:         mesh = ALE::Generator::refine(mesh, user.refinementLimit, true);
154:       } else {
155:         mesh = ALE::Generator::refine(mesh, refineLimit, (void *) &user, true);
156:       }
157:       ALE::LogStagePop(stage);
158:       PetscSynchronizedPrintf(comm, "  [%d]Generated %d local elements\n", mesh->commRank(), mesh->getTopology()->heightStratum(0)->size());
159:       PetscSynchronizedPrintf(comm, "  [%d]Generated %d local vertices\n", mesh->commRank(), mesh->getTopology()->depthStratum(0)->size());
160:       PetscSynchronizedFlush(comm);
161:     }
162:     topology = mesh->getTopology();

164:     stage = ALE::LogStageRegister("BndValues");
165:     ALE::LogStagePush(stage);
166:     PetscPrintf(comm, "Calculating boundary values\n");
167:     ALE::Obj<ALE::Mesh::field_type> boundary = mesh->getBoundary();
168:     ALE::Obj<ALE::Mesh::sieve_type::traits::depthSequence> vertices = topology->depthStratum(0);
169:     ALE::Mesh::field_type::patch_type groundPatch(0, 1);
170:     ALE::Mesh::field_type::patch_type voltagePatch(0, 3);
171:     ALE::Mesh::field_type::patch_type patch;

173:     for(ALE::Mesh::sieve_type::traits::depthSequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
174:       if (boundary->getIndex(groundPatch, *v_iter).index > 0) {
175:         double values[1] = {0.0};

177:         boundary->update(groundPatch, *v_iter, values);
178:       } else if (boundary->getIndex(voltagePatch, *v_iter).index > 0) {
179:         double values[1] = {user.voltage};

181:         boundary->update(voltagePatch, *v_iter, values);
182:       }
183:     }
184:     if (debug) {boundary->view("Mesh Boundary");}
185:     ALE::LogStagePop(stage);

187:     user.numCharges = 1;
188:     PetscMalloc2(user.numCharges,PetscReal,&user.charge,user.numCharges*2,PetscReal,&user.chargeLocation);
189:     user.charge[0] = 1.0;
190:     user.chargeLocation[0] = 0.0;
191:     user.chargeLocation[1] = 0.0;

193:     ALE::Obj<ALE::Mesh::field_type> u = mesh->getField("u");
194:     ALE::Obj<ALE::Mesh::field_type> b = mesh->getField("b");
195:     u->setPatch(topology->leaves(), ALE::Mesh::field_type::patch_type());
196:     u->setFiberDimensionByDepth(patch, 0, 1);
197:     u->orderPatches();
198:     if (debug) {u->view("u");}
199:     u->createGlobalOrder();
200:     b->setPatch(topology->leaves(), ALE::Mesh::field_type::patch_type());
201:     b->setFiberDimensionByDepth(patch, 0, 1);
202:     b->orderPatches();
203:     if (debug) {b->view("b");}
204:     b->createGlobalOrder();
205:     ALE::Obj<ALE::Mesh::sieve_type::traits::heightSequence> elements = topology->heightStratum(0);
206:     ALE::Obj<ALE::Mesh::bundle_type> vertexBundle = mesh->getBundle(0);
207:     std::string orderName("element");

209:     for(ALE::Mesh::sieve_type::traits::heightSequence::iterator e_iter = elements->begin(); e_iter != elements->end(); e_iter++) {
210:       // setFiberDimensionByDepth() does not work here since we only want it to apply to the patch cone
211:       //   What we really need is the depthStratum relative to the patch
212:       ALE::Obj<ALE::Mesh::bundle_type::order_type::coneSequence> cone = vertexBundle->getPatch(orderName, *e_iter);

214:       u->setPatch(orderName, cone, *e_iter);
215:       b->setPatch(orderName, cone, *e_iter);
216:       for(ALE::Mesh::bundle_type::order_type::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
217:         u->setFiberDimension(orderName, *e_iter, *c_iter, 1);
218:         b->setFiberDimension(orderName, *e_iter, *c_iter, 1);
219:       }
220:     }
221:     u->orderPatches(orderName);
222:     b->orderPatches(orderName);
223:     CheckElementGeometry(mesh);

225:     CreateDielectricField(mesh, mesh->getField("epsilon"));
226:     user.epsilon = mesh->getField("epsilon");

228:     Mesh petscMesh;
229:     MeshCreate(comm, &petscMesh);
230:     MeshSetMesh(petscMesh, mesh);
231:     DMMGCreate(comm,3,PETSC_NULL,&dmmg);
232:     DMMGSetDM(dmmg, (DM) petscMesh);
233:     MeshDestroy(petscMesh);
234:     for (l = 0; l < DMMGGetLevels(dmmg); l++) {
235:       DMMGSetUser(dmmg,l,&user);
236:     }

238:     DMMGSetKSP(dmmg,ComputeRHS,ComputeMatrix);
239:     MeshGetGlobalScatter(mesh, "u", DMMGGetx(dmmg), &user.injection);

241:     DMMGSolve(dmmg);

243:     if (debug) {
244:       PetscPrintf(mesh->comm(), "Solution vector:");
245:       VecView(DMMGGetx(dmmg), PETSC_VIEWER_STDOUT_WORLD);
246:     }
247:     //ComputeError(dmmg[DMMGGetLevels(dmmg)-1], DMMGGetx(dmmg), &error);
248:     //PetscPrintf(comm,"Error norm %g\n",error);

250:     MatMult(DMMGGetJ(dmmg),DMMGGetx(dmmg),DMMGGetr(dmmg));
251:     VecAXPY(DMMGGetr(dmmg),-1.0,DMMGGetRHS(dmmg));
252:     VecNorm(DMMGGetr(dmmg),NORM_2,&norm);
253:     PetscPrintf(comm,"Residual norm %g\n",norm);
254:     VecAssemblyBegin(DMMGGetx(dmmg));
255:     VecAssemblyEnd(DMMGGetx(dmmg));

257:     stage = ALE::LogStageRegister("MeshOutput");
258:     ALE::LogStagePush(stage);
259:     PetscPrintf(comm, "Creating VTK mesh file\n");
260:     PetscViewerCreate(comm, &viewer);
261:     PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);
262:     PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_VTK);
263:     PetscViewerFileSetName(viewer, "channel.vtk");
264:     MeshView_Sieve_Newer(mesh, viewer);
265:     //VecView(DMMGGetRHS(dmmg), viewer);
266:     if (viewEnergy) {
267:       ALE::Obj<ALE::Mesh::field_type> u = mesh->getField("u");
268:       Vec energy, locU;

270:       VecCreateSeqWithArray(PETSC_COMM_SELF, u->getSize(patch), u->restrict(patch), &locU);
271:       VecScatterBegin(user.injection, DMMGGetx(dmmg), locU, INSERT_VALUES, SCATTER_REVERSE);
272:       VecScatterEnd(user.injection, DMMGGetx(dmmg), locU, INSERT_VALUES, SCATTER_REVERSE);
273:       VecDestroy(locU);

275:       PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK_CELL);
276:       CreateEnergyDensity(mesh, u, mesh->getField("epsilon"), &energy);
277:       VecView(energy, viewer);
278:       PetscViewerPopFormat(viewer);
279:     } else {
280:       VecView(DMMGGetx(dmmg), viewer);
281:     }
282:     PetscViewerDestroy(viewer);
283:     ALE::LogStagePop(stage);

285:     PetscFree2(user.charge,user.chargeLocation);
286:     DMMGDestroy(dmmg);
287:   } catch (ALE::Exception e) {
288:     std::cout << e << std::endl;
289:   }
290:   PetscFinalize();
291:   return 0;
292: }

294: double refineLimit(const double centroid[], void *ctx) {
295:   UserContext *user = (UserContext *) ctx;
296:   double       r2   = centroid[0]*centroid[0] + centroid[1]*centroid[1];

298:   return user->refinementLimit*pow(r2, user->refinementExp*0.5);
299: }

303: /*
304:   2D radially-symmetric channel boundary:
305:                     29                   30    31                 32
306:                      V                    V     V                  V
307:     2----------------------------------3-----4----12--------------------------------13
308:     |                                  |     |     |                                 |
309:     |                                  |     |     |                                 |
310:     |                               39>|     |     |<46                              |
311:     |                                  |     |     |                                 |
312:  28>|                                  |  59>|     |                                 |<33
313:     |                                  8     |     11                                |
314:     |                               40>\     |     /<45                              |
315:     |                                   7 42 |43 10                                  |
316:     |                                 41>\ V | V /<44                                |
317:     |               55                    6--5--9                 57                 |
318:     |                V                       |<56                  V                 |
319:     1----------------------------------------O--------------------------------------14
320:     |                                        |<58                                    |
321:     |                                    25-20--19                                   |
322:     |                                 49>/ ^ | ^ \<52                                |
323:     |                                   24 50| 51 18                                 |
324:     |                               48>/     |     \<53                              |
325:  27>|                                  23    |    17                                 |<34
326:     |                                  |  60>|     |                                 |
327:     |                                  |     |     |                                 |
328:     |                               47>|     |     |<54                              |
329:     |               38                 | 37  | 36  |              35                 |
330:     |                V                 |  V  |  V  |               V                 |
331:     26(X)-----------------------------22----21-----16-------------------------------15

333:     (X) denotes the last vertex, (XX) denotes the last edge
334: */
335: PetscErrorCode CreateMeshBoundary(ALE::Obj<ALE::Mesh> mesh)
336: {
337:   ALE::Obj<ALE::Mesh::sieve_type> topology = mesh->getTopology();
338:   PetscScalar       coords[54] =  {/*O*/      0.0,      0.0,
339:                                    /*1*/   -112.5,      0.0,
340:                                    /*2*/   -112.5,     50.0,
341:                                    /*3*/    -12.5,     50.0,
342:                                    /*4*/      0.0,     50.0,
343:                                    /*5*/      0.0,      3.0,
344:                                    /*6*/     -2.5,      3.0,
345:                                    /*7*/   -35.0/6.0,  10.0,
346:                                    /*8*/    -12.5,     15.0,
347:                                    /*9*/      2.5,      3.0,
348:                                    /*10*/   35.0/6.0,  10.0,
349:                                    /*11*/    12.5,     15.0,
350:                                    /*12*/    12.5,     50.0,
351:                                    /*13*/   112.5,     50.0,
352:                                    /*14*/   112.5,      0.0,
353:                                    /*15*/   112.5,    -50.0,
354:                                    /*16*/    12.5,    -50.0,
355:                                    /*17*/    12.5,    -15.0,
356:                                    /*18*/   35.0/6.0, -10.0,
357:                                    /*19*/     2.5,     -3.0,
358:                                    /*20*/     0.0,     -3.0,
359:                                    /*21*/     0.0,    -50.0,
360:                                    /*22*/   -12.5,    -50.0,
361:                                    /*23*/   -12.5,    -15.0,
362:                                    /*24*/  -35.0/6.0, -10.0,
363:                                    /*25*/    -2.5,     -3.0,
364:                                    /*26*/  -112.5,    -50.0};
365:   PetscInt    connectivity[68] = {26, 1, /* 1: phi = 0 */
366:                                   1, 2,  /* 1: phi = 0 */
367:                                   2, 3,  /* 2: grad phi = 0 */
368:                                   3, 4,  /* 2: grad phi = 0 */
369:                                   4, 12, /* 2: grad phi = 0 */
370:                                   12,13, /* 2: grad phi = 0 */
371:                                   13,14, /* 3: phi = V */
372:                                   14,15, /* 3: phi = V */
373:                                   15,16, /* 4: grad phi = 0 */
374:                                   16,21, /* 4: grad phi = 0 */
375:                                   21,22, /* 4: grad phi = 0 */
376:                                   22,26, /* 4: grad phi = 0 */
377:                                   3, 8,  /* 5: top lipid boundary */
378:                                   8, 7,  /* 5: top lipid boundary */
379:                                   7, 6,  /* 5: top lipid boundary */
380:                                   6, 5,  /* 5: top lipid boundary */
381:                                   5,  9, /* 5: top lipid boundary */
382:                                   9, 10, /* 5: top lipid boundary */
383:                                   10,11, /* 5: top lipid boundary */
384:                                   11,12, /* 5: top lipid boundary */
385:                                   22,23, /* 6: bottom lipid boundary */
386:                                   23,24, /* 6: bottom lipid boundary */
387:                                   24,25, /* 6: bottom lipid boundary */
388:                                   25,20, /* 6: bottom lipid boundary */
389:                                   20,19, /* 6: bottom lipid boundary */
390:                                   19,18, /* 6: bottom lipid boundary */
391:                                   18,17, /* 6: bottom lipid boundary */
392:                                   17,16, /* 6: bottom lipid boundary */
393:                                   0, 1,  /* 7: symmetry preservation */
394:                                   0, 5,  /* 7: symmetry preservation */
395:                                   0, 14, /* 7: symmetry preservation */
396:                                   0, 20, /* 7: symmetry preservation */
397:                                   4, 5,  /* 7: symmetry preservation */
398:                                   21,20  /* 7: symmetry preservation */
399:                                   };
400:   ALE::Mesh::point_type vertices[27];

403:   PetscInt order = 0;
404:   if (mesh->commRank() == 0) {
405:     ALE::Mesh::point_type edge;

407:     /* Create topology and ordering */
408:     for(int v = 0; v < 27; v++) {
409:       vertices[v] = ALE::Mesh::point_type(0, v);
410:     }
411:     for(int e = 27; e < 61; e++) {
412:       int ee = e - 27;
413:       edge = ALE::Mesh::point_type(0, e);
414:       topology->addArrow(vertices[connectivity[2*ee]],   edge, order++);
415:       topology->addArrow(vertices[connectivity[2*ee+1]], edge, order++);
416:     }
417:   }
418:   topology->stratify();
419:   mesh->createVertexBundle(34, connectivity, 27);
420:   mesh->createSerialCoordinates(2, 0, coords);
421:   /* Create boundary conditions */
422:   if (mesh->commRank() == 0) {
423:     for(int e = 47; e < 55; e++) {
424:       int ee = e - 27;

426:       topology->setMarker(ALE::Mesh::point_type(0, e), 6);
427:       topology->setMarker(ALE::Mesh::point_type(0, connectivity[2*ee+0]), 6);
428:       topology->setMarker(ALE::Mesh::point_type(0, connectivity[2*ee+1]), 6);
429:     }
430:     for(int e = 39; e < 47; e++) {
431:       int ee = e - 27;

433:       topology->setMarker(ALE::Mesh::point_type(0, e), 5);
434:       topology->setMarker(ALE::Mesh::point_type(0, connectivity[2*ee+0]), 5);
435:       topology->setMarker(ALE::Mesh::point_type(0, connectivity[2*ee+1]), 5);
436:     }
437:     for(int e = 35; e < 39; e++) {
438:       int ee = e - 27;

440:       topology->setMarker(ALE::Mesh::point_type(0, e), 4);
441:       topology->setMarker(ALE::Mesh::point_type(0, connectivity[2*ee+0]), 4);
442:       topology->setMarker(ALE::Mesh::point_type(0, connectivity[2*ee+1]), 4);
443:     }
444:     for(int e = 29; e < 33; e++) {
445:       int ee = e - 27;

447:       topology->setMarker(ALE::Mesh::point_type(0, e), 2);
448:       topology->setMarker(ALE::Mesh::point_type(0, connectivity[2*ee+0]), 2);
449:       topology->setMarker(ALE::Mesh::point_type(0, connectivity[2*ee+1]), 2);
450:     }
451:     for(int e = 33; e < 35; e++) {
452:       int ee = e - 27;

454:       topology->setMarker(ALE::Mesh::point_type(0, e), 3);
455:       topology->setMarker(ALE::Mesh::point_type(0, connectivity[2*ee+0]), 3);
456:       topology->setMarker(ALE::Mesh::point_type(0, connectivity[2*ee+1]), 3);
457:     }
458:     for(int e = 27; e < 29; e++) {
459:       int ee = e - 27;

461:       topology->setMarker(ALE::Mesh::point_type(0, e), 1);
462:       topology->setMarker(ALE::Mesh::point_type(0, connectivity[2*ee+0]), 1);
463:       topology->setMarker(ALE::Mesh::point_type(0, connectivity[2*ee+1]), 1);
464:     }
465:   }
466:   return(0);
467: }

469: #define NUM_QUADRATURE_POINTS 9

471: /* Quadrature points */
472: static double points[18] = {
473:   -0.794564690381,
474:   -0.822824080975,
475:   -0.866891864322,
476:   -0.181066271119,
477:   -0.952137735426,
478:   0.575318923522,
479:   -0.0885879595127,
480:   -0.822824080975,
481:   -0.409466864441,
482:   -0.181066271119,
483:   -0.787659461761,
484:   0.575318923522,
485:   0.617388771355,
486:   -0.822824080975,
487:   0.0479581354402,
488:   -0.181066271119,
489:   -0.623181188096,
490:   0.575318923522};

492: /* Quadrature weights */
493: static double weights[9] = {
494:   0.223257681932,
495:   0.2547123404,
496:   0.0775855332238,
497:   0.357212291091,
498:   0.407539744639,
499:   0.124136853158,
500:   0.223257681932,
501:   0.2547123404,
502:   0.0775855332238};

504: #define NUM_BASIS_FUNCTIONS 3

506: /* Nodal basis function evaluations */
507: static double Basis[27] = {
508:   0.808694385678,
509:   0.10271765481,
510:   0.0885879595127,
511:   0.52397906772,
512:   0.0665540678392,
513:   0.409466864441,
514:   0.188409405952,
515:   0.0239311322871,
516:   0.787659461761,
517:   0.455706020244,
518:   0.455706020244,
519:   0.0885879595127,
520:   0.29526656778,
521:   0.29526656778,
522:   0.409466864441,
523:   0.10617026912,
524:   0.10617026912,
525:   0.787659461761,
526:   0.10271765481,
527:   0.808694385678,
528:   0.0885879595127,
529:   0.0665540678392,
530:   0.52397906772,
531:   0.409466864441,
532:   0.0239311322871,
533:   0.188409405952,
534:   0.787659461761};

536: /* Nodal basis function derivative evaluations */
537: static double BasisDerivatives[54] = {
538:   -0.5,
539:   -0.5,
540:   0.5,
541:   4.74937635818e-17,
542:   0.0,
543:   0.5,
544:   -0.5,
545:   -0.5,
546:   0.5,
547:   4.74937635818e-17,
548:   0.0,
549:   0.5,
550:   -0.5,
551:   -0.5,
552:   0.5,
553:   4.74937635818e-17,
554:   0.0,
555:   0.5,
556:   -0.5,
557:   -0.5,
558:   0.5,
559:   4.74937635818e-17,
560:   0.0,
561:   0.5,
562:   -0.5,
563:   -0.5,
564:   0.5,
565:   4.74937635818e-17,
566:   0.0,
567:   0.5,
568:   -0.5,
569:   -0.5,
570:   0.5,
571:   4.74937635818e-17,
572:   0.0,
573:   0.5,
574:   -0.5,
575:   -0.5,
576:   0.5,
577:   4.74937635818e-17,
578:   0.0,
579:   0.5,
580:   -0.5,
581:   -0.5,
582:   0.5,
583:   4.74937635818e-17,
584:   0.0,
585:   0.5,
586:   -0.5,
587:   -0.5,
588:   0.5,
589:   4.74937635818e-17,
590:   0.0,
591:   0.5};

595: PetscErrorCode ElementGeometry(ALE::Obj<ALE::Mesh> mesh, const ALE::Mesh::point_type& e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
596: {
597:   const double  *coords = mesh->getCoordinates()->restrict(std::string("element"), e);
598:   int            dim = mesh->getDimension();
599:   PetscReal      det, invDet;

602:   if (debug) {
603:     MPI_Comm comm = mesh->comm();
604:     int      rank = mesh->commRank();

606:     PetscSynchronizedPrintf(comm, "[%d]Element (%d, %d)\n", rank, e.prefix, e.index);
607:     PetscSynchronizedPrintf(comm, "[%d]Coordinates:\n[%d]  ", rank, rank);
608:     for(int f = 0; f <= dim; f++) {
609:       PetscSynchronizedPrintf(comm, " (");
610:       for(int d = 0; d < dim; d++) {
611:         if (d > 0) PetscSynchronizedPrintf(comm, ", ");
612:         PetscSynchronizedPrintf(comm, "%g", coords[f*dim+d]);
613:       }
614:       PetscSynchronizedPrintf(comm, ")");
615:     }
616:     PetscSynchronizedPrintf(comm, "\n");
617:   }
618:   if (v0) {
619:     for(int d = 0; d < dim; d++) {
620:       v0[d] = coords[d];
621:     }
622:   }
623:   if (J) {
624:     for(int d = 0; d < dim; d++) {
625:       for(int f = 0; f < dim; f++) {
626:         J[d*dim+f] = 0.5*(coords[(f+1)*dim+d] - coords[0*dim+d]);
627:       }
628:     }
629:     if (debug) {
630:       MPI_Comm comm = mesh->comm();
631:       int      rank = mesh->commRank();

633:       for(int d = 0; d < dim; d++) {
634:         if (d == 0) {
635:           PetscSynchronizedPrintf(comm, "[%d]J = /", rank);
636:         } else if (d == dim-1) {
637:           PetscSynchronizedPrintf(comm, "[%d]    \\", rank);
638:         } else {
639:           PetscSynchronizedPrintf(comm, "[%d]    |", rank);
640:         }
641:         for(int e = 0; e < dim; e++) {
642:           PetscSynchronizedPrintf(comm, " %g", J[d*dim+e]);
643:         }
644:         if (d == 0) {
645:           PetscSynchronizedPrintf(comm, " \\\n");
646:         } else if (d == dim-1) {
647:           PetscSynchronizedPrintf(comm, " /\n");
648:         } else {
649:           PetscSynchronizedPrintf(comm, " |\n");
650:         }
651:       }
652:     }
653:     if (dim == 2) {
654:       det = J[0]*J[3] - J[1]*J[2];
655:     } else if (dim == 3) {
656:       det = J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
657:             J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
658:             J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
659:     }
660:     invDet = 1.0/det;
661:     if (detJ) {
662:       if (det < 0) {SETERRQ(PETSC_ERR_ARG_WRONG, "Negative Matrix determinant");}
663:       *detJ = det;
664:     }
665:     if (invJ) {
666:       if (dim == 2) {
667:         invJ[0] =  invDet*J[3];
668:         invJ[1] = -invDet*J[1];
669:         invJ[2] = -invDet*J[2];
670:         invJ[3] =  invDet*J[0];
671:       } else if (dim == 3) {
672:         // FIX: This may be wrong
673:         invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
674:         invJ[0*3+1] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
675:         invJ[0*3+2] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
676:         invJ[1*3+0] = invDet*(J[0*3+1]*J[2*3+2] - J[0*3+2]*J[2*3+1]);
677:         invJ[1*3+1] = invDet*(J[0*3+2]*J[2*3+0] - J[0*3+0]*J[2*3+2]);
678:         invJ[1*3+2] = invDet*(J[0*3+0]*J[2*3+1] - J[0*3+1]*J[2*3+0]);
679:         invJ[2*3+0] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
680:         invJ[2*3+1] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
681:         invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
682:       }
683:       if (debug) {
684:         MPI_Comm comm = mesh->comm();
685:         int      rank = mesh->commRank();

687:         for(int d = 0; d < dim; d++) {
688:           if (d == 0) {
689:             PetscSynchronizedPrintf(comm, "[%d]Jinv = /", rank);
690:           } else if (d == dim-1) {
691:             PetscSynchronizedPrintf(comm, "[%d]       \\", rank);
692:           } else {
693:             PetscSynchronizedPrintf(comm, "[%d]       |", rank);
694:           }
695:           for(int e = 0; e < dim; e++) {
696:             PetscSynchronizedPrintf(comm, " %g", invJ[d*dim+e]);
697:           }
698:           if (d == 0) {
699:             PetscSynchronizedPrintf(comm, " \\\n");
700:           } else if (d == dim-1) {
701:             PetscSynchronizedPrintf(comm, " /\n");
702:           } else {
703:             PetscSynchronizedPrintf(comm, " |\n");
704:           }
705:         }
706:       }
707:     }
708:   }
709:   return(0);
710: }

714: PetscErrorCode ElementContains(ALE::Obj<ALE::Mesh> mesh, const ALE::Mesh::point_type& e, double point[], PetscTruth *contains)
715: {
716:   const double  *coords = mesh->getCoordinates()->restrict(std::string("element"), e);
717:   int            dim = mesh->getDimension();

720:   *contains = PETSC_TRUE;
721:   for(int p = 0; p < 3; p++) {
722:     double side[2] = {coords[((p+1)%3)*dim+0] - coords[p*dim+0], coords[((p+1)%3)*dim+1] - coords[p*dim+1]};
723:     double ray[2]  = {point[0]                - coords[p*dim+0], point[1]                - coords[p*dim+1]};
724:     double cross   = side[0]*ray[1] - side[1]*ray[0];

726:     if (cross <= 0.0) {
727:       *contains = PETSC_FALSE;
728:       break;
729:     }
730:   }
731:   return(0);
732: }

736: PetscErrorCode CheckElementGeometry(ALE::Obj<ALE::Mesh> mesh)
737: {
738:   ALE::Obj<ALE::Mesh::sieve_type::traits::heightSequence> elements = mesh->getTopology()->heightStratum(0);
739:   PetscInt       dim = mesh->getDimension();
740:   PetscReal     *v0, *Jac;
741:   PetscReal      detJ;

745:   PetscMalloc2(dim,PetscReal,&v0,dim*dim,PetscReal,&Jac);
746:   for(ALE::Mesh::sieve_type::traits::heightSequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
747:     ElementGeometry(mesh, *e_iter, v0, Jac, PETSC_NULL, &detJ);
748:   }
749:   PetscSynchronizedFlush(mesh->comm());
750:   PetscFree2(v0,Jac);
751:   return(0);
752: }

756: PetscErrorCode ComputeChargeDensity(ALE::Obj<ALE::Mesh> mesh, double points[], double jac[], double v0[], double detJ, double basis[], ALE::Obj<ALE::Mesh::field_type> field, UserContext *user)
757: {
758:   ALE::Obj<ALE::Mesh::sieve_type::traits::heightSequence> elements = mesh->getTopology()->heightStratum(0);
759:   PetscReal      elementVec[NUM_BASIS_FUNCTIONS];
760:   PetscTruth     contains;

764:   for(int q = 0; q < user->numCharges; q++) {
765:     double x = user->chargeLocation[q*2+0];
766:     double y = user->chargeLocation[q*2+1];

768:     for(ALE::Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); e_itor++) {
769:       ElementContains(mesh, *e_itor, user->chargeLocation, &contains);
770:       if (contains) {
771:         double     minDist = 1.0e30;
772:         int        minP    = -1;
773:         PetscReal *v0, *Jac, detJ;

775:         ElementGeometry(mesh, *e_itor, v0, Jac, PETSC_NULL, &detJ);
776:         for(int p = 0; p < NUM_QUADRATURE_POINTS; p++) {
777:           double xi, eta, x_p, y_p, dist;

779:           xi = points[p*2+0] + 1.0;
780:           eta = points[p*2+1] + 1.0;
781:           x_p = jac[0]*xi + jac[1]*eta + v0[0];
782:           y_p = jac[2]*xi + jac[3]*eta + v0[1];
783:           dist = (x - x_p)*(x - x_p) + (y - y_p)*(y - y_p);
784:           if (dist < minDist) {
785:             minDist = dist;
786:             minP    = p;
787:           }
788:         }
789:         for(int f = 0; f < NUM_BASIS_FUNCTIONS; f++) {
790:           elementVec[f] += basis[minP*NUM_BASIS_FUNCTIONS+f]*user->charge[q]*detJ;
791:         }
792:         /* Assembly */
793:         field->updateAdd("element", *e_itor, elementVec);
794:         break;
795:       }
796:     }
797:   }
798:   return(0);
799: }

803: PetscErrorCode ComputeRHS(DMMG dmmg, Vec b)
804: {
805:   ALE::Obj<ALE::Mesh> m;
806:   Mesh                mesh = (Mesh) dmmg->dm;
807:   UserContext        *user = (UserContext *) dmmg->user;
808:   MPI_Comm            comm;
809:   PetscReal           elementVec[NUM_BASIS_FUNCTIONS];
810:   PetscReal           *v0, *Jac;
811:   PetscReal           xi, eta, x_q, y_q, detJ, funcValue;
812:   PetscInt            dim;
813:   PetscInt            f, q;
814:   PetscErrorCode      ierr;

817:   PetscObjectGetComm((PetscObject) mesh, &comm);
818:   MeshGetMesh(mesh, &m);
819:   dim  = m->getDimension();
820:   PetscMalloc(dim * sizeof(PetscReal), &v0);
821:   PetscMalloc(dim*dim * sizeof(PetscReal), &Jac);
822:   ALE::Obj<ALE::Mesh::field_type> field = m->getField("b");
823:   ALE::Obj<ALE::Mesh::sieve_type::traits::heightSequence> elements = m->getTopology()->heightStratum(0);
824:   ALE::Mesh::field_type::patch_type patch;
825:   for(ALE::Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); e_itor++) {
826:     ElementGeometry(m, *e_itor, v0, Jac, PETSC_NULL, &detJ);
827:     /* Element integral */
828:     PetscMemzero(elementVec, NUM_BASIS_FUNCTIONS*sizeof(PetscScalar));
829:     for(q = 0; q < NUM_QUADRATURE_POINTS; q++) {
830:       xi = points[q*2+0] + 1.0;
831:       eta = points[q*2+1] + 1.0;
832:       x_q = Jac[0]*xi + Jac[1]*eta + v0[0];
833:       y_q = Jac[2]*xi + Jac[3]*eta + v0[1];
834:       funcValue = 0.0;
835:       for(f = 0; f < NUM_BASIS_FUNCTIONS; f++) {
836:         elementVec[f] += Basis[q*NUM_BASIS_FUNCTIONS+f]*funcValue*weights[q]*detJ;
837:       }
838:     }
839:     if (debug) {PetscSynchronizedPrintf(comm, "elementVec = [%g %g %g]\n", elementVec[0], elementVec[1], elementVec[2]);}
840:     /* Assembly */
841:     field->updateAdd("element", *e_itor, elementVec);
842:     if (debug) {PetscSynchronizedFlush(comm);}
843:   }
844:   ComputeChargeDensity(m, points, v0, Jac, detJ, Basis, field, user);
845:   PetscFree(v0);
846:   PetscFree(Jac);

848:   Vec locB;
849:   VecCreateSeqWithArray(PETSC_COMM_SELF, field->getSize(patch), field->restrict(patch), &locB);
850:   VecScatterBegin(user.injection, locB, b, ADD_VALUES, SCATTER_FORWARD);
851:   VecScatterEnd(user.injection, locB, b, ADD_VALUES, SCATTER_FORWARD);
852:   VecDestroy(locB);

854:   {
855:     /* Zero out BC rows */
856:     ALE::Mesh::field_type::patch_type patch;
857:     ALE::Mesh::field_type::patch_type groundPatch(0, 1);
858:     ALE::Mesh::field_type::patch_type voltagePatch(0, 3);
859:     ALE::Obj<ALE::Mesh::field_type> boundary = m->getBoundary();
860:     ALE::Obj<ALE::Mesh::field_type::order_type::coneSequence> groundCone = boundary->getPatch(groundPatch);
861:     ALE::Obj<ALE::Mesh::field_type::order_type::coneSequence> voltageCone = boundary->getPatch(voltagePatch);
862:     PetscScalar *boundaryValues;
863:     PetscInt    *boundaryIndices;
864:     PetscInt     numBoundaryIndices = 0;
865:     PetscInt     k = 0;

867:     for(ALE::Mesh::field_type::order_type::coneSequence::iterator p = groundCone->begin(); p != groundCone->end(); ++p) {
868:       numBoundaryIndices += field->getGlobalOrder()->getIndex(patch, *p).index;
869:     }
870:     for(ALE::Mesh::field_type::order_type::coneSequence::iterator p = voltageCone->begin(); p != voltageCone->end(); ++p) {
871:       numBoundaryIndices += field->getGlobalOrder()->getIndex(patch, *p).index;
872:     }
873:     PetscMalloc2(numBoundaryIndices,PetscInt,&boundaryIndices,numBoundaryIndices,PetscScalar,&boundaryValues);
874:     for(ALE::Mesh::field_type::order_type::coneSequence::iterator p = groundCone->begin(); p != groundCone->end(); ++p) {
875:       const ALE::Mesh::field_type::index_type& idx = field->getGlobalOrder()->getIndex(patch, *p);
876:       const double *data = boundary->restrict(groundPatch, *p);

878:       for(int i = 0; i < idx.index; i++) {
879:         boundaryIndices[k] = idx.prefix + i;
880:         boundaryValues[k] = data[i];
881:         k++;
882:       }
883:     }
884:     for(ALE::Mesh::field_type::order_type::coneSequence::iterator p = voltageCone->begin(); p != voltageCone->end(); ++p) {
885:       const ALE::Mesh::field_type::index_type& idx = field->getGlobalOrder()->getIndex(patch, *p);
886:       const double *data = boundary->restrict(voltagePatch, *p);

888:       for(int i = 0; i < idx.index; i++) {
889:         boundaryIndices[k] = idx.prefix + i;
890:         boundaryValues[k] = data[i];
891:         k++;
892:       }
893:     }
894:     if (debug) {
895:       boundary->view("Boundary for rhs conditions");
896:       for(int i = 0; i < numBoundaryIndices; i++) {
897:         PetscSynchronizedPrintf(comm, "[%d]boundaryIndices[%d] = %d\n", m->commRank(), i, boundaryIndices[i]);
898:       }
899:     }
900:     PetscSynchronizedFlush(comm);
901:     VecSetValues(b, numBoundaryIndices, boundaryIndices, boundaryValues, INSERT_VALUES);
902:     PetscFree2(boundaryIndices, boundaryValues);
903:   }
904:   if (debug) {
905:     PetscPrintf(comm, "Rhs vector:");
906:     VecView(b, PETSC_VIEWER_STDOUT_WORLD);
907:   }
908:   return(0);
909: }

913: PetscErrorCode ComputeMatrix(DMMG dmmg, Mat J, Mat jac)
914: {
915:   ALE::Obj<ALE::Mesh> m;
916:   Mesh              mesh = (Mesh) dmmg->dm;
917:   UserContext      *user = (UserContext *) dmmg->user;
918:   MPI_Comm          comm;
919:   PetscReal         elementMat[NUM_BASIS_FUNCTIONS*NUM_BASIS_FUNCTIONS];
920:   PetscReal        *v0, *Jac, *Jinv, *t_der, *b_der;
921:   PetscReal         xi, eta, x_q, y_q, detJ;
922:   PetscInt          dim;
923:   PetscInt          f, g, q;
924:   PetscMPIInt       rank;
925:   PetscErrorCode    ierr;

928:   PetscObjectGetComm((PetscObject) mesh, &comm);
929:   MPI_Comm_rank(comm, &rank);
930:   MeshGetMesh(mesh, &m);
931:   dim  = m->getDimension();
932:   PetscMalloc(dim * sizeof(PetscReal), &v0);
933:   PetscMalloc(dim * sizeof(PetscReal), &t_der);
934:   PetscMalloc(dim * sizeof(PetscReal), &b_der);
935:   PetscMalloc(dim*dim * sizeof(PetscReal), &Jac);
936:   PetscMalloc(dim*dim * sizeof(PetscReal), &Jinv);
937:   ALE::Obj<ALE::Mesh::field_type> field = m->getField("u");
938:   ALE::Obj<ALE::Mesh::sieve_type::traits::heightSequence> elements = m->getTopology()->heightStratum(0);
939:   for(ALE::Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); e_itor++) {
940:     ALE::Mesh::field_type::patch_type patch;
941:     double eps = user->epsilon->restrict(patch, *e_itor)[0];

943:     CHKMEMQ;
944:     ElementGeometry(m, *e_itor, v0, Jac, Jinv, &detJ);
945:     /* Element integral */
946:     PetscMemzero(elementMat, NUM_BASIS_FUNCTIONS*NUM_BASIS_FUNCTIONS*sizeof(PetscScalar));
947:     for(q = 0; q < NUM_QUADRATURE_POINTS; q++) {
948:       xi = points[q*2+0] + 1.0;
949:       eta = points[q*2+1] + 1.0;
950:       x_q = Jac[0]*xi + Jac[1]*eta + v0[0];
951:       y_q = Jac[2]*xi + Jac[3]*eta + v0[1];
952:       for(f = 0; f < NUM_BASIS_FUNCTIONS; f++) {
953:         t_der[0] = Jinv[0]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+0] + Jinv[2]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+1];
954:         t_der[1] = Jinv[1]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+0] + Jinv[3]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+1];
955:         for(g = 0; g < NUM_BASIS_FUNCTIONS; g++) {
956:           b_der[0] = Jinv[0]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+0] + Jinv[2]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+1];
957:           b_der[1] = Jinv[1]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+0] + Jinv[3]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+1];
958:           elementMat[f*NUM_BASIS_FUNCTIONS+g] += eps*(t_der[0]*b_der[0] + t_der[1]*b_der[1])*weights[q]*detJ;
959:         }
960:       }
961:     }
962:     if (debug) {
963:       PetscSynchronizedPrintf(comm, "[%d]elementMat = [%g %g %g]\n                [%g %g %g]\n                [%g %g %g]\n",
964:                                      rank, elementMat[0], elementMat[1], elementMat[2], elementMat[3], elementMat[4],
965:                                      elementMat[5], elementMat[6], elementMat[7], elementMat[8]);
966:     }
967:     /* Assembly */
968:     updateOperator(jac, field, *e_itor, elementMat, ADD_VALUES);
969:     if (debug) {PetscSynchronizedFlush(comm);}
970:   }
971:   PetscFree(v0);
972:   PetscFree(t_der);
973:   PetscFree(b_der);
974:   PetscFree(Jac);
975:   PetscFree(Jinv);
976:   MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
977:   MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);

979:   {
980:     /* Zero out BC rows */
981:     ALE::Mesh::field_type::patch_type patch;
982:     ALE::Mesh::field_type::patch_type groundPatch(0, 1);
983:     ALE::Mesh::field_type::patch_type voltagePatch(0, 3);
984:     ALE::Obj<ALE::Mesh::field_type> boundary = m->getBoundary();
985:     ALE::Obj<ALE::Mesh::field_type::order_type::coneSequence> groundCone = boundary->getPatch(groundPatch);
986:     ALE::Obj<ALE::Mesh::field_type::order_type::coneSequence> voltageCone = boundary->getPatch(voltagePatch);
987:     PetscInt *boundaryIndices;
988:     PetscInt  numBoundaryIndices = 0;
989:     PetscInt  k = 0;

991:     for(ALE::Mesh::field_type::order_type::coneSequence::iterator p = groundCone->begin(); p != groundCone->end(); ++p) {
992:       numBoundaryIndices += field->getGlobalOrder()->getIndex(patch, *p).index;
993:     }
994:     for(ALE::Mesh::field_type::order_type::coneSequence::iterator p = voltageCone->begin(); p != voltageCone->end(); ++p) {
995:       numBoundaryIndices += field->getGlobalOrder()->getIndex(patch, *p).index;
996:     }
997:     PetscMalloc(numBoundaryIndices * sizeof(PetscInt), &boundaryIndices);
998:     for(ALE::Mesh::field_type::order_type::coneSequence::iterator p = groundCone->begin(); p != groundCone->end(); ++p) {
999:       const ALE::Mesh::field_type::index_type& idx = field->getGlobalOrder()->getIndex(patch, *p);

1001:       for(int i = 0; i < idx.index; i++) {
1002:         boundaryIndices[k++] = idx.prefix + i;
1003:       }
1004:     }
1005:     for(ALE::Mesh::field_type::order_type::coneSequence::iterator p = voltageCone->begin(); p != voltageCone->end(); ++p) {
1006:       const ALE::Mesh::field_type::index_type& idx = field->getGlobalOrder()->getIndex(patch, *p);

1008:       for(int i = 0; i < idx.index; i++) {
1009:         boundaryIndices[k++] = idx.prefix + i;
1010:       }
1011:     }
1012:     if (debug) {
1013:       for(int i = 0; i < numBoundaryIndices; i++) {
1014:         PetscSynchronizedPrintf(comm, "[%d]boundaryIndices[%d] = %d\n", rank, i, boundaryIndices[i]);
1015:       }
1016:     }
1017:     PetscSynchronizedFlush(comm);
1018:     MatZeroRows(jac, numBoundaryIndices, boundaryIndices, 1.0);
1019:     PetscFree(boundaryIndices);
1020:   }
1021:   MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
1022:   MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
1023:   return(0);
1024: }

1028: PetscErrorCode ComputeDielectric(double x, double y, double *epsilon) {
1029:   double water   = 80.0;
1030:   double lipid   = 2.0;
1031:   double channel = 10.0;

1034:   *epsilon = -80;
1035:   if ((x >= -112.5) && (x <= -12.5)) {
1036:     // Left water bath
1037:     *epsilon = water;
1038:   } else if ((x >= 12.5) && (x <= 112.5)) {
1039:     // Right water bath
1040:     *epsilon = water;
1041:   } else {
1042:     if ((y >= 15.0) && (y <= 50.0)) {
1043:       // Top lipid
1044:       *epsilon = lipid;
1045:     } else if ((y <= -15.0) && (y >= -50.0)) {
1046:       // Bottom lipid
1047:       *epsilon = lipid;
1048:     } else {
1049:       if ((x >= -12.5) && (x <= -2.5)) {
1050:         // Left lipid or water
1051:         if (x <= -35.0/6.0) {
1052:           // Left parallelogram
1053:           if (y >= 0.0) {
1054:             // Top half
1055:             double slope = (15.0 - 10.0)/(-12.5 + 35.0/6.0);

1057:             if (y <= 15.0 + slope*(x + 12.5)) {
1058:               // Middle water
1059:               *epsilon = water;
1060:             } else {
1061:               // Middle lipid
1062:               *epsilon = lipid;
1063:             }
1064:           } else {
1065:             // Bottom half
1066:             double slope = (-15.0 + 10.0)/(-12.5 + 35.0/6.0);

1068:             if (y >= -15.0 + slope*(x + 12.5)) {
1069:               // Middle water
1070:               *epsilon = water;
1071:             } else {
1072:               // Middle lipid
1073:               *epsilon = lipid;
1074:             }
1075:           }
1076:         } else {
1077:           // Right parallelogram
1078:           if (y >= 0.0) {
1079:             // Top half
1080:             double slope = (10.0 - 3.0)/(-35.0/6.0 + 2.5);

1082:             if (y <= 10.0 + slope*(x + 35.0/6.0)) {
1083:               // Middle water
1084:               *epsilon = water;
1085:             } else {
1086:               // Middle lipid
1087:               *epsilon = lipid;
1088:             }
1089:           } else {
1090:             // Bottom half
1091:             double slope = (-10.0 + 3.0)/(-35.0/6.0 + 2.5);

1093:             if (y >= -10.0 + slope*(x + 35.0/6.0)) {
1094:               // Middle water
1095:               *epsilon = water;
1096:             } else {
1097:               // Middle lipid
1098:               *epsilon = lipid;
1099:             }
1100:           }
1101:         }
1102:       } else if ((x >= 2.5) && (x <= 12.5)) {
1103:         // Right lipid or water
1104:         if (x >= 35.0/6.0) {
1105:           // Right parallelogram
1106:           if (y >= 0.0) {
1107:             // Top half
1108:             double slope = (15.0 - 10.0)/(12.5 - 35.0/6.0);

1110:             if (y <= 15.0 + slope*(x - 12.5)) {
1111:               // Middle water
1112:               *epsilon = water;
1113:             } else {
1114:               // Middle lipid
1115:               *epsilon = lipid;
1116:             }
1117:           } else {
1118:             // Bottom half
1119:             double slope = (-15.0 + 10.0)/(12.5 - 35.0/6.0);

1121:             if (y >= -15.0 + slope*(x - 12.5)) {
1122:               // Middle water
1123:               *epsilon = water;
1124:             } else {
1125:               // Middle lipid
1126:               *epsilon = lipid;
1127:             }
1128:           }
1129:         } else {
1130:           // Left parallelogram
1131:           if (y >= 0.0) {
1132:             // Top half
1133:             double slope = (10.0 - 3.0)/(35.0/6.0 - 2.5);

1135:             if (y <= 10.0 + slope*(x - 35.0/6.0)) {
1136:               // Middle water
1137:               *epsilon = water;
1138:             } else {
1139:               // Middle lipid
1140:               *epsilon = lipid;
1141:             }
1142:           } else {
1143:             // Bottom half
1144:             double slope = (-10.0 + 3.0)/(35.0/6.0 - 2.5);

1146:             if (y >= -10.0 + slope*(x - 35.0/6.0)) {
1147:               // Middle water
1148:               *epsilon = water;
1149:             } else {
1150:               // Middle lipid
1151:               *epsilon = lipid;
1152:             }
1153:           }
1154:         }
1155:       } else {
1156:         if ((y <= 3.0) && (y >= -3.0)) {
1157:           // Channel
1158:           *epsilon = channel;
1159:         } else {
1160:           // Central lipid
1161:           *epsilon = lipid;
1162:         }
1163:       }
1164:     }
1165:   }
1166:   return(0);
1167: }

1171: /*
1172:   Creates a vector whose value is the dielectric constant on each element
1173: */
1174: PetscErrorCode CreateDielectricField(ALE::Obj<ALE::Mesh> mesh, ALE::Obj<ALE::Mesh::field_type> epsilon)
1175: {
1176:   ALE::Obj<ALE::Mesh::sieve_type> topology = mesh->getTopology();
1177:   ALE::Obj<ALE::Mesh::sieve_type::traits::heightSequence> elements = topology->heightStratum(0);
1178:   ALE::Obj<ALE::Mesh::field_type> coordinates = mesh->getCoordinates();
1179:   ALE::Mesh::field_type::patch_type patch;
1180:   std::string orderName("element");

1184:   ALE_LOG_EVENT_BEGIN;
1185:   epsilon->setPatch(topology->leaves(), patch);
1186:   epsilon->setFiberDimensionByHeight(patch, 0, 1);
1187:   epsilon->orderPatches();
1188:   epsilon->createGlobalOrder();

1190:   for(ALE::Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); ++e_itor) {
1191:     const double *coords = coordinates->restrict(orderName, *e_itor);
1192:     double centroidX = (coords[0]+coords[2]+coords[4])/3.0;
1193:     double centroidY = (coords[1]+coords[3]+coords[5])/3.0;

1195:     double eps;
1196:     ComputeDielectric(centroidX, centroidY, &eps);
1197:     epsilon->update(patch, *e_itor, &eps);
1198:   }
1199:   ALE_LOG_EVENT_END;
1200:   return(0);
1201: }

1205: /*
1206:   Creates a vector whose value is the energy on each element
1207: */
1208: PetscErrorCode CreateEnergyDensity(ALE::Obj<ALE::Mesh> mesh, ALE::Obj<ALE::Mesh::field_type> field, ALE::Obj<ALE::Mesh::field_type> epsilon, Vec *energy)
1209: {
1210:   ALE::Obj<ALE::Mesh::sieve_type> topology = mesh->getTopology();
1211:   ALE::Obj<ALE::Mesh::sieve_type::traits::heightSequence> elements = topology->heightStratum(0);
1212:   ALE::Obj<ALE::Mesh::field_type> coordinates = mesh->getCoordinates();
1213:   ALE::Obj<ALE::Mesh::field_type> e2 = mesh->getField("energy");
1214:   ALE::Mesh::field_type::patch_type patch;
1215:   std::string orderName("element");
1216:   PetscInt       dim = mesh->getDimension();
1217:   PetscReal      elementMat[NUM_BASIS_FUNCTIONS*NUM_BASIS_FUNCTIONS];
1218:   PetscReal     *v0, *Jac, *Jinv, *t_der, *b_der;
1219:   PetscReal      xi, eta, x_q, y_q, detJ;

1223:   ALE_LOG_EVENT_BEGIN;
1224:   e2->setPatch(mesh->getTopology()->leaves(), patch);
1225:   e2->setFiberDimensionByHeight(patch, 0, 1);
1226:   e2->orderPatches();
1227:   e2->createGlobalOrder();

1229:   PetscMalloc(dim * sizeof(PetscReal), &v0);
1230:   PetscMalloc(dim * sizeof(PetscReal), &t_der);
1231:   PetscMalloc(dim * sizeof(PetscReal), &b_der);
1232:   PetscMalloc(dim*dim * sizeof(PetscReal), &Jac);
1233:   PetscMalloc(dim*dim * sizeof(PetscReal), &Jinv);
1234:   for(ALE::Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); ++e_itor) {
1235:     const double *phi = field->restrict(orderName, *e_itor);
1236:     double eps = epsilon->restrict(patch, *e_itor)[0];
1237:     double elementEnergy = 0.0;

1239:     ElementGeometry(mesh, *e_itor, v0, Jac, Jinv, &detJ);
1240:     /* Element integral */
1241:     PetscMemzero(elementMat, NUM_BASIS_FUNCTIONS*NUM_BASIS_FUNCTIONS*sizeof(PetscScalar));
1242:     for(int q = 0; q < NUM_QUADRATURE_POINTS; q++) {
1243:       xi = points[q*2+0] + 1.0;
1244:       eta = points[q*2+1] + 1.0;
1245:       x_q = Jac[0]*xi + Jac[1]*eta + v0[0];
1246:       y_q = Jac[2]*xi + Jac[3]*eta + v0[1];
1247:       for(int f = 0; f < NUM_BASIS_FUNCTIONS; f++) {
1248:         t_der[0] = Jinv[0]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+0] + Jinv[2]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+1];
1249:         t_der[1] = Jinv[1]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+0] + Jinv[3]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+1];
1250:         for(int g = 0; g < NUM_BASIS_FUNCTIONS; g++) {
1251:           b_der[0] = Jinv[0]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+0] + Jinv[2]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+1];
1252:           b_der[1] = Jinv[1]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+0] + Jinv[3]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+1];
1253:           elementMat[f*NUM_BASIS_FUNCTIONS+g] += eps*(t_der[0]*b_der[0] + t_der[1]*b_der[1])*weights[q]*detJ;
1254:         }
1255:       }
1256:     }
1257:     for(int f = 0; f < NUM_BASIS_FUNCTIONS; f++) {
1258:       for(int g = 0; g < NUM_BASIS_FUNCTIONS; g++) {
1259:         elementEnergy += phi[f]*elementMat[f*NUM_BASIS_FUNCTIONS+g]*phi[g];
1260:       }
1261:     }
1262:     e2->update(patch, *e_itor, &elementEnergy);
1263:   }
1264:   PetscFree(v0);
1265:   PetscFree(t_der);
1266:   PetscFree(b_der);
1267:   PetscFree(Jac);
1268:   PetscFree(Jinv);

1270:   VecScatter injection;
1271:   MeshCreateVector(mesh, mesh->getBundle(mesh->getDimension()), energy);
1272:   MeshGetGlobalScatter(mesh, "energy", *energy, &injection);

1274:   Vec locEnergy;
1275:   VecCreateSeqWithArray(PETSC_COMM_SELF, e2->getSize(patch), e2->restrict(patch), &locEnergy);
1276:   VecScatterBegin(injection, locEnergy, *energy, INSERT_VALUES, SCATTER_FORWARD);
1277:   VecScatterEnd(injection, locEnergy, *energy, INSERT_VALUES, SCATTER_FORWARD);
1278:   VecDestroy(locEnergy);
1279:   ALE_LOG_EVENT_END;
1280:   return(0);
1281: }