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: }