Actual source code: ex33.c
1: /*T
2: Concepts: KSP^solving a system of linear equations
3: Concepts: KSP^Laplacian, 2d
4: Processors: n
5: T*/
7: /*
8: Laplacian in 2D. Modeled by the partial differential equation
10: -div grad u = f, 0 < x,y < 2,
12: with forcing function
14: f = 2
15: f = 2y (1 - y) + 2x (1 - x)
17: with Dirichlet boundary conditions
19: u = 0 for x = 0, x = 2 and u = x (2 - x) for y = 0, y = 2
20: u = 0 for x = 0, x = 2, y = 0, y = 2
22: or pure Neumman boundary conditions.
24: This has exact solution
26: u = x (2 - x)
27: u = x (2 - x) y (2 - y)
29: This uses multigrid to solve the linear system
31: The 2D test mesh
33: 13
34: 14--29----31---12
35: |\ |\ |
36: 2 2 5 2 3 7 3
37: 7 6 8 0 2
38: | 4 \ | 6 \ |
39: | \| \|
40: 15--20-16-24---11
41: |\ |\ |
42: 1 1 1 2 2 3 2
43: 8 7 1 2 5
44: | 0 \ | 2 \ |
45: | \| \|
46: 8--19----23---10
47: 9
49: */
51: static char help[] = "Solves 2D inhomogeneous Laplacian using multigrid.\n\n";
53: #include petscmesh.h
54: #include petscksp.h
55: #include petscmg.h
56: #include petscdmmg.h
58: PetscErrorCode MeshView_Sieve_Newer(ALE::Obj<ALE::Two::Mesh>, PetscViewer);
59: PetscErrorCode CreateMeshBoundary(ALE::Obj<ALE::Two::Mesh>);
60: PetscErrorCode updateOperator(Mat, ALE::Obj<ALE::Two::Mesh::field_type>, const ALE::Two::Mesh::point_type&, PetscScalar [], InsertMode);
67: typedef enum {DIRICHLET, NEUMANN} BCType;
69: typedef struct {
70: BCType bcType;
71: VecScatter injection;
72: } UserContext;
74: PetscInt debug;
78: int main(int argc,char **argv)
79: {
80: MPI_Comm comm;
81: DMMG *dmmg;
82: UserContext user;
83: PetscViewer viewer;
84: const char *bcTypes[2] = {"dirichlet", "neumann"};
85: PetscReal refinementLimit, norm, error;
86: PetscInt dim, bc, l, meshDebug;
89: PetscInitialize(&argc,&argv,(char *)0,help);
90: comm = PETSC_COMM_WORLD;
91: PetscOptionsBegin(comm, "", "Options for the inhomogeneous Poisson equation", "DMMG");
92: debug = 0;
93: PetscOptionsInt("-debug", "The debugging flag", "ex33.c", 0, &debug, PETSC_NULL);
94: meshDebug = 0;
95: PetscOptionsInt("-mesh_debug", "The mesh debugging flag", "ex33.c", 0, &meshDebug, PETSC_NULL);
96: dim = 2;
97: PetscOptionsInt("-dim", "The mesh dimension", "ex33.c", 2, &dim, PETSC_NULL);
98: refinementLimit = 0.0;
99: PetscOptionsReal("-refinement_limit", "The area of the largest triangle in the mesh", "ex33.c", 1.0, &refinementLimit, PETSC_NULL);
100: bc = (PetscInt)DIRICHLET;
101: PetscOptionsEList("-bc_type","Type of boundary condition","ex33.c",bcTypes,2,bcTypes[0],&bc,PETSC_NULL);
102: user.bcType = (BCType) bc;
103: PetscOptionsEnd();
105: ALE::Obj<ALE::Two::Mesh> meshBoundary = ALE::Two::Mesh(comm, dim-1, meshDebug);
106: ALE::Obj<ALE::Two::Mesh> mesh;
108: try {
109: ALE::LogStage stage = ALE::LogStageRegister("MeshCreation");
110: ALE::LogStagePush(stage);
111: PetscPrintf(comm, "Generating mesh\n");
112: CreateMeshBoundary(meshBoundary);
113: mesh = ALE::Two::Generator::generate(meshBoundary);
114: ALE::Obj<ALE::Two::Mesh::sieve_type> topology = mesh->getTopology();
115: PetscPrintf(comm, " Generated %d elements\n", topology->heightStratum(0)->size());
116: PetscPrintf(comm, " Generated %d vertices\n", topology->depthStratum(0)->size());
117: ALE::LogStagePop(stage);
119: stage = ALE::LogStageRegister("MeshDistribution");
120: ALE::LogStagePush(stage);
121: PetscPrintf(comm, "Distributing mesh\n");
122: mesh = mesh->distribute();
123: ALE::LogStagePop(stage);
125: if (refinementLimit > 0.0) {
126: stage = ALE::LogStageRegister("MeshRefine");
127: ALE::LogStagePush(stage);
128: PetscPrintf(comm, "Refining mesh\n");
129: mesh = ALE::Two::Generator::refine(mesh, refinementLimit);
130: ALE::LogStagePop(stage);
131: }
132: topology = mesh->getTopology();
134: stage = ALE::LogStageRegister("BndValues");
135: ALE::LogStagePush(stage);
136: PetscPrintf(comm, "Calculating boundary values\n");
137: ALE::Obj<ALE::Two::Mesh::field_type> boundary = mesh->getBoundary();
138: ALE::Obj<ALE::Two::Mesh::sieve_type::traits::depthSequence> vertices = topology->depthStratum(0);
139: ALE::Two::Mesh::field_type::patch_type bdPatch(0, 1);
140: ALE::Two::Mesh::field_type::patch_type patch;
142: for(ALE::Two::Mesh::sieve_type::traits::depthSequence::iterator v_iter = vertices->begin(); v_iter != vertices->end(); ++v_iter) {
143: if (boundary->getIndex(bdPatch, *v_iter).index > 0) {
144: const double *coords = mesh->getCoordinates()->restrict(patch, *v_iter);
145: double values[1];
147: values[0] = coords[0]*(2.0 - coords[0]);
148: boundary->update(bdPatch, *v_iter, values);
149: }
150: }
151: if (debug) {boundary->view("Mesh Boundary");}
152: ALE::LogStagePop(stage);
154: ALE::Obj<ALE::Two::Mesh::field_type> u = mesh->getField("u");
155: ALE::Obj<ALE::Two::Mesh::field_type> b = mesh->getField("b");
156: u->setPatch(topology->leaves(), ALE::Two::Mesh::field_type::patch_type());
157: u->setFiberDimensionByDepth(patch, 0, 1);
158: u->orderPatches();
159: if (debug) {u->view("u");}
160: u->createGlobalOrder();
161: b->setPatch(topology->leaves(), ALE::Two::Mesh::field_type::patch_type());
162: b->setFiberDimensionByDepth(patch, 0, 1);
163: b->orderPatches();
164: if (debug) {b->view("b");}
165: b->createGlobalOrder();
166: ALE::Obj<ALE::Two::Mesh::sieve_type::traits::heightSequence> elements = topology->heightStratum(0);
167: ALE::Obj<ALE::Two::Mesh::bundle_type> vertexBundle = mesh->getBundle(0);
168: std::string orderName("element");
170: for(ALE::Two::Mesh::sieve_type::traits::heightSequence::iterator e_iter = elements->begin(); e_iter != elements->end(); e_iter++) {
171: // setFiberDimensionByDepth() does not work here since we only want it to apply to the patch cone
172: // What we really need is the depthStratum relative to the patch
173: ALE::Obj<ALE::Two::Mesh::bundle_type::order_type::coneSequence> cone = vertexBundle->getPatch(orderName, *e_iter);
175: u->setPatch(orderName, cone, *e_iter);
176: b->setPatch(orderName, cone, *e_iter);
177: for(ALE::Two::Mesh::bundle_type::order_type::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
178: u->setFiberDimension(orderName, *e_iter, *c_iter, 1);
179: b->setFiberDimension(orderName, *e_iter, *c_iter, 1);
180: }
181: }
182: u->orderPatches(orderName);
183: b->orderPatches(orderName);
184: CheckElementGeometry(mesh);
186: Mesh petscMesh;
187: MeshCreate(comm, &petscMesh);
188: MeshSetMesh(petscMesh, mesh);
189: DMMGCreate(comm,3,PETSC_NULL,&dmmg);
190: DMMGSetDM(dmmg, (DM) petscMesh);
191: MeshDestroy(petscMesh);
192: for (l = 0; l < DMMGGetLevels(dmmg); l++) {
193: DMMGSetUser(dmmg,l,&user);
194: }
196: DMMGSetKSP(dmmg,ComputeRHS,ComputeMatrix);
197: if (user.bcType == NEUMANN) {
198: DMMGSetNullSpace(dmmg,PETSC_TRUE,0,PETSC_NULL);
199: }
200: MeshGetGlobalScatter(mesh, "u", DMMGGetx(dmmg), &user.injection);
202: DMMGSolve(dmmg);
204: if (debug) {
205: PetscPrintf(mesh->comm(), "Solution vector:");
206: VecView(DMMGGetx(dmmg), PETSC_VIEWER_STDOUT_WORLD);
207: }
208: ComputeError(dmmg[DMMGGetLevels(dmmg)-1], DMMGGetx(dmmg), &error);
209: PetscPrintf(comm,"Error norm %g\n",error);
211: MatMult(DMMGGetJ(dmmg),DMMGGetx(dmmg),DMMGGetr(dmmg));
212: VecAXPY(DMMGGetr(dmmg),-1.0,DMMGGetRHS(dmmg));
213: VecNorm(DMMGGetr(dmmg),NORM_2,&norm);
214: PetscPrintf(comm,"Residual norm %g\n",norm);
215: VecAssemblyBegin(DMMGGetx(dmmg));
216: VecAssemblyEnd(DMMGGetx(dmmg));
218: stage = ALE::LogStageRegister("MeshOutput");
219: ALE::LogStagePush(stage);
220: PetscPrintf(comm, "Creating VTK mesh file\n");
221: PetscViewerCreate(comm, &viewer);
222: PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);
223: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_VTK);
224: PetscViewerFileSetName(viewer, "poisson.vtk");
225: MeshView_Sieve_Newer(mesh, viewer);
226: //VecView(DMMGGetRHS(dmmg), viewer);
227: VecView(DMMGGetx(dmmg), viewer);
228: PetscViewerDestroy(viewer);
229: ALE::LogStagePop(stage);
231: DMMGDestroy(dmmg);
232: } catch (ALE::Exception e) {
233: std::cout << e << std::endl;
234: }
235: PetscFinalize();
236: return 0;
237: }
241: /*
242: Simple square boundary:
244: 6--14-5--13-4
245: | | |
246: 15 19 12
247: | | |
248: 7--20-8--18-3
249: | | |
250: 16 17 11
251: | | |
252: 0--9--1--10-2
253: */
254: PetscErrorCode CreateSquareBoundary(ALE::Obj<ALE::Two::Mesh> mesh)
255: {
256: MPI_Comm comm = mesh->comm();
257: ALE::Obj<ALE::Two::Mesh::sieve_type> topology = mesh->getTopology();
258: PetscScalar coords[18] = {0.0, 0.0,
259: 1.0, 0.0,
260: 2.0, 0.0,
261: 2.0, 1.0,
262: 2.0, 2.0,
263: 1.0, 2.0,
264: 0.0, 2.0,
265: 0.0, 1.0,
266: 1.0, 1.0};
267: PetscInt connectivity[40] = {0, 1,
268: 1, 2,
269: 2, 3,
270: 3, 4,
271: 4, 5,
272: 5, 6,
273: 6, 7,
274: 7, 0,
275: 1, 8,
276: 3, 8,
277: 5, 8,
278: 7, 8};
279: ALE::Two::Mesh::point_type vertices[9];
280: PetscInt order = 0;
281: PetscMPIInt rank;
282: PetscErrorCode ierr;
285: MPI_Comm_rank(comm, &rank);
286: if (rank == 0) {
287: ALE::Two::Mesh::point_type edge;
289: /* Create topology and ordering */
290: for(int v = 0; v < 9; v++) {
291: vertices[v] = ALE::Two::Mesh::point_type(0, v);
292: }
293: for(int e = 9; e < 17; e++) {
294: edge = ALE::Two::Mesh::point_type(0, e);
295: topology->addArrow(vertices[e-9], edge, order++);
296: topology->addArrow(vertices[(e-8)%8], edge, order++);
297: }
298: edge = ALE::Two::Mesh::point_type(0, 17);
299: topology->addArrow(vertices[1], edge, order++);
300: topology->addArrow(vertices[8], edge, order++);
301: edge = ALE::Two::Mesh::point_type(0, 18);
302: topology->addArrow(vertices[3], edge, order++);
303: topology->addArrow(vertices[8], edge, order++);
304: edge = ALE::Two::Mesh::point_type(0, 19);
305: topology->addArrow(vertices[5], edge, order++);
306: topology->addArrow(vertices[8], edge, order++);
307: edge = ALE::Two::Mesh::point_type(0, 20);
308: topology->addArrow(vertices[7], edge, order++);
309: topology->addArrow(vertices[8], edge, order++);
310: }
311: topology->stratify();
312: mesh->createVertexBundle(20, connectivity);
313: mesh->createSerialCoordinates(2, 0, coords);
314: /* Create boundary conditions */
315: if (rank == 0) {
316: for(int v = 0; v < 8; v++) {
317: topology->setMarker(ALE::Two::Mesh::point_type(0, v), 1);
318: }
319: for(int e = 9; e < 17; e++) {
320: topology->setMarker(ALE::Two::Mesh::point_type(0, e), 1);
321: }
322: }
323: return(0);
324: }
328: /*
329: Simple cube boundary:
331: 7-----6
332: /| /|
333: 3-----2 |
334: | | | |
335: | 4---|-5
336: |/ |/
337: 0-----1
338: */
339: PetscErrorCode CreateCubeBoundary(ALE::Obj<ALE::Two::Mesh> mesh)
340: {
341: MPI_Comm comm = mesh->comm();
342: ALE::Obj<ALE::Two::Mesh::sieve_type> topology = mesh->getTopology();
343: PetscScalar coords[24] = {0.0, 0.0, 0.0,
344: 1.0, 0.0, 0.0,
345: 1.0, 1.0, 0.0,
346: 0.0, 1.0, 0.0,
347: 0.0, 0.0, 1.0,
348: 1.0, 0.0, 1.0,
349: 1.0, 1.0, 1.0,
350: 0.0, 1.0, 1.0};
351: PetscInt connectivity[24] = {0, 1, 2, 3,
352: 7, 6, 5, 4,
353: 0, 4, 5, 1,
354: 1, 5, 6, 2,
355: 2, 6, 7, 3,
356: 3, 7, 4, 0};
357: ALE::Obj<std::set<ALE::Two::Mesh::point_type> > cone = std::set<ALE::Two::Mesh::point_type>();
358: ALE::Two::Mesh::point_type vertices[8];
359: ALE::Two::Mesh::point_type edges[12];
360: ALE::Two::Mesh::point_type edge;
361: PetscInt embedDim = 3;
362: PetscInt order = 0;
363: PetscMPIInt rank;
364: PetscErrorCode ierr;
367: MPI_Comm_rank(comm, &rank);
368: if (rank == 0) {
369: ALE::Two::Mesh::point_type face;
371: /* Create topology and ordering */
372: /* Vertices: 0 .. 3 on the bottom of the cube, 4 .. 7 on the top */
373: for(int v = 0; v < 8; v++) {
374: vertices[v] = ALE::Two::Mesh::point_type(0, v);
375: }
377: /* Edges on the bottom: Sieve element numbers e = 8 .. 11, edge numbers e - 8 = 0 .. 3 */
378: for(int e = 8; e < 12; e++) {
379: edge = ALE::Two::Mesh::point_type(0, e);
380: edges[e-8] = edge;
381: topology->addArrow(vertices[e-8], edge, order++);
382: topology->addArrow(vertices[(e-7)%4], edge, order++);
383: }
384: /* Edges on the top: Sieve element numbers e = 12 .. 15, edge numbers e - 8 = 4 .. 7 */
385: for(int e = 12; e < 16; e++) {
386: edge = ALE::Two::Mesh::point_type(0, e);
387: edges[e-8] = edge;
388: topology->addArrow(vertices[e-8], edge, order++);
389: topology->addArrow(vertices[(e-11)%4+4], edge, order++);
390: }
391: /* Edges from bottom to top: Sieve element numbers e = 16 .. 19, edge numbers e - 8 = 8 .. 11 */
392: for(int e = 16; e < 20; e++) {
393: edge = ALE::Two::Mesh::point_type(0, e);
394: edges[e-8] = edge;
395: topology->addArrow(vertices[e-16], edge, order++);
396: topology->addArrow(vertices[e-16+4], edge, order++);
397: }
399: /* Bottom face */
400: face = ALE::Two::Mesh::point_type(0, 20);
401: topology->addArrow(edges[0], face, order++);
402: topology->addArrow(edges[1], face, order++);
403: topology->addArrow(edges[2], face, order++);
404: topology->addArrow(edges[3], face, order++);
405: /* Top face */
406: face = ALE::Two::Mesh::point_type(0, 21);
407: topology->addArrow(edges[4], face, order++);
408: topology->addArrow(edges[5], face, order++);
409: topology->addArrow(edges[6], face, order++);
410: topology->addArrow(edges[7], face, order++);
411: /* Side faces: f = 22 .. 25 */
412: for(int f = 22; f < 26; f++) {
413: face = ALE::Two::Mesh::point_type(0, f);
414: int v = f - 22;
415: /* Covered by edges f - 22, f - 22 + 4, f - 22 + 8, (f - 21)%4 + 8 */
416: topology->addArrow(edges[v], face, order++);
417: topology->addArrow(edges[(v+1)%4+8], face, order++);
418: topology->addArrow(edges[v+4], face, order++);
419: topology->addArrow(edges[v+8], face, order++);
420: }
421: }/* if(rank == 0) */
422: topology->stratify();
423: if (rank == 0) {
424: ALE::Obj<ALE::Two::Mesh::bundle_type> vertexBundle = mesh->getBundle(0);
425: ALE::Obj<ALE::Two::Mesh::bundle_type::PointArray> points = ALE::Two::Mesh::bundle_type::PointArray();
426: const std::string orderName("element");
427: /* Bottom face */
428: ALE::Two::Mesh::point_type face = ALE::Two::Mesh::point_type(0, 20);
429: points->clear();
430: points->push_back(vertices[0]);
431: points->push_back(vertices[1]);
432: points->push_back(vertices[2]);
433: points->push_back(vertices[3]);
434: vertexBundle->setPatch(orderName, points, face);
435: /* Top face */
436: face = ALE::Two::Mesh::point_type(0, 21);
437: points->clear();
438: points->push_back(vertices[4]);
439: points->push_back(vertices[5]);
440: points->push_back(vertices[6]);
441: points->push_back(vertices[7]);
442: vertexBundle->setPatch(orderName, points, face);
443: /* Side faces: f = 22 .. 25 */
444: for(int f = 22; f < 26; f++) {
445: face = ALE::Two::Mesh::point_type(0, f);
446: int v = f - 22;
447: /* Covered by edges f - 22, f - 22 + 4, f - 22 + 8, (f - 21)%4 + 8 */
448: points->clear();
449: points->push_back(vertices[v]);
450: points->push_back(vertices[(v+1)%4]);
451: points->push_back(vertices[(v+1)%4+4]);
452: points->push_back(vertices[v+4]);
453: vertexBundle->setPatch(orderName, points, face);
454: }
455: }
456: mesh->createVertexBundle(6, connectivity);
457: mesh->createSerialCoordinates(embedDim, 0, coords);
459: /* Create boundary conditions: set marker 1 to all of the sieve elements,
460: since everything is on the boundary (no internal faces, edges or vertices) */
461: if (rank == 0) {
462: /* set marker to the base of the topology sieve -- the faces and the edges */
463: topology->setMarker(topology->base(), 1);
464: /* set marker to the vertices -- the 0-depth stratum */
465: topology->setMarker(topology->depthStratum(0), 1);
466: }
468: return(0);
469: }
473: /*
474: Simple square boundary:
476: 6--14-5--13-4
477: | | |
478: 15 19 12
479: | | |
480: 7--20-8--18-3
481: | | |
482: 16 17 11
483: | | |
484: 0--9--1--10-2
485: */
486: PetscErrorCode CreateMeshBoundary(ALE::Obj<ALE::Two::Mesh> mesh)
487: {
488: int dim = mesh->getDimension();
492: if (dim == 1) {
493: CreateSquareBoundary(mesh);
494: } else if (dim == 2) {
495: CreateCubeBoundary(mesh);
496: } else {
497: SETERRQ1(PETSC_ERR_SUP, "Cannot construct a boundary of dimension %d", dim);
498: }
499: return(0);
500: }
502: #ifndef MESH_3D
504: #define NUM_QUADRATURE_POINTS 9
506: /* Quadrature points */
507: static double points[18] = {
508: -0.794564690381,
509: -0.822824080975,
510: -0.866891864322,
511: -0.181066271119,
512: -0.952137735426,
513: 0.575318923522,
514: -0.0885879595127,
515: -0.822824080975,
516: -0.409466864441,
517: -0.181066271119,
518: -0.787659461761,
519: 0.575318923522,
520: 0.617388771355,
521: -0.822824080975,
522: 0.0479581354402,
523: -0.181066271119,
524: -0.623181188096,
525: 0.575318923522};
527: /* Quadrature weights */
528: static double weights[9] = {
529: 0.223257681932,
530: 0.2547123404,
531: 0.0775855332238,
532: 0.357212291091,
533: 0.407539744639,
534: 0.124136853158,
535: 0.223257681932,
536: 0.2547123404,
537: 0.0775855332238};
539: #define NUM_BASIS_FUNCTIONS 3
541: /* Nodal basis function evaluations */
542: static double Basis[27] = {
543: 0.808694385678,
544: 0.10271765481,
545: 0.0885879595127,
546: 0.52397906772,
547: 0.0665540678392,
548: 0.409466864441,
549: 0.188409405952,
550: 0.0239311322871,
551: 0.787659461761,
552: 0.455706020244,
553: 0.455706020244,
554: 0.0885879595127,
555: 0.29526656778,
556: 0.29526656778,
557: 0.409466864441,
558: 0.10617026912,
559: 0.10617026912,
560: 0.787659461761,
561: 0.10271765481,
562: 0.808694385678,
563: 0.0885879595127,
564: 0.0665540678392,
565: 0.52397906772,
566: 0.409466864441,
567: 0.0239311322871,
568: 0.188409405952,
569: 0.787659461761};
571: /* Nodal basis function derivative evaluations */
572: static double BasisDerivatives[54] = {
573: -0.5,
574: -0.5,
575: 0.5,
576: 4.74937635818e-17,
577: 0.0,
578: 0.5,
579: -0.5,
580: -0.5,
581: 0.5,
582: 4.74937635818e-17,
583: 0.0,
584: 0.5,
585: -0.5,
586: -0.5,
587: 0.5,
588: 4.74937635818e-17,
589: 0.0,
590: 0.5,
591: -0.5,
592: -0.5,
593: 0.5,
594: 4.74937635818e-17,
595: 0.0,
596: 0.5,
597: -0.5,
598: -0.5,
599: 0.5,
600: 4.74937635818e-17,
601: 0.0,
602: 0.5,
603: -0.5,
604: -0.5,
605: 0.5,
606: 4.74937635818e-17,
607: 0.0,
608: 0.5,
609: -0.5,
610: -0.5,
611: 0.5,
612: 4.74937635818e-17,
613: 0.0,
614: 0.5,
615: -0.5,
616: -0.5,
617: 0.5,
618: 4.74937635818e-17,
619: 0.0,
620: 0.5,
621: -0.5,
622: -0.5,
623: 0.5,
624: 4.74937635818e-17,
625: 0.0,
626: 0.5};
628: #else
630: #define NUM_QUADRATURE_POINTS 27
632: /* Quadrature points */
633: static double points[81] = {
634: -0.809560240317,
635: -0.835756864273,
636: -0.854011951854,
637: -0.865851516496,
638: -0.884304792128,
639: -0.305992467923,
640: -0.939397037651,
641: -0.947733495427,
642: 0.410004419777,
643: -0.876607962782,
644: -0.240843539439,
645: -0.854011951854,
646: -0.913080888692,
647: -0.465239359176,
648: -0.305992467923,
649: -0.960733394129,
650: -0.758416359732,
651: 0.410004419777,
652: -0.955631394718,
653: 0.460330056095,
654: -0.854011951854,
655: -0.968746121484,
656: 0.0286773243482,
657: -0.305992467923,
658: -0.985880737721,
659: -0.53528439884,
660: 0.410004419777,
661: -0.155115591937,
662: -0.835756864273,
663: -0.854011951854,
664: -0.404851369974,
665: -0.884304792128,
666: -0.305992467923,
667: -0.731135462175,
668: -0.947733495427,
669: 0.410004419777,
670: -0.452572254354,
671: -0.240843539439,
672: -0.854011951854,
673: -0.61438408645,
674: -0.465239359176,
675: -0.305992467923,
676: -0.825794030022,
677: -0.758416359732,
678: 0.410004419777,
679: -0.803159052121,
680: 0.460330056095,
681: -0.854011951854,
682: -0.861342428212,
683: 0.0286773243482,
684: -0.305992467923,
685: -0.937360010468,
686: -0.53528439884,
687: 0.410004419777,
688: 0.499329056443,
689: -0.835756864273,
690: -0.854011951854,
691: 0.0561487765469,
692: -0.884304792128,
693: -0.305992467923,
694: -0.522873886699,
695: -0.947733495427,
696: 0.410004419777,
697: -0.0285365459258,
698: -0.240843539439,
699: -0.854011951854,
700: -0.315687284208,
701: -0.465239359176,
702: -0.305992467923,
703: -0.690854665916,
704: -0.758416359732,
705: 0.410004419777,
706: -0.650686709523,
707: 0.460330056095,
708: -0.854011951854,
709: -0.753938734941,
710: 0.0286773243482,
711: -0.305992467923,
712: -0.888839283216,
713: -0.53528439884,
714: 0.410004419777};
716: /* Quadrature weights */
717: static double weights[27] = {
718: 0.0701637994372,
719: 0.0653012061324,
720: 0.0133734490519,
721: 0.0800491405774,
722: 0.0745014590358,
723: 0.0152576273199,
724: 0.0243830167241,
725: 0.022693189565,
726: 0.0046474825267,
727: 0.1122620791,
728: 0.104481929812,
729: 0.021397518483,
730: 0.128078624924,
731: 0.119202334457,
732: 0.0244122037118,
733: 0.0390128267586,
734: 0.0363091033041,
735: 0.00743597204272,
736: 0.0701637994372,
737: 0.0653012061324,
738: 0.0133734490519,
739: 0.0800491405774,
740: 0.0745014590358,
741: 0.0152576273199,
742: 0.0243830167241,
743: 0.022693189565,
744: 0.0046474825267};
746: #define NUM_BASIS_FUNCTIONS 4
748: /* Nodal basis function evaluations */
749: static double Basis[108] = {
750: 0.749664528222,
751: 0.0952198798417,
752: 0.0821215678634,
753: 0.0729940240731,
754: 0.528074388273,
755: 0.0670742417521,
756: 0.0578476039361,
757: 0.347003766038,
758: 0.23856305665,
759: 0.0303014811743,
760: 0.0261332522867,
761: 0.705002209888,
762: 0.485731727037,
763: 0.0616960186091,
764: 0.379578230281,
765: 0.0729940240731,
766: 0.342156357896,
767: 0.0434595556538,
768: 0.267380320412,
769: 0.347003766038,
770: 0.154572667042,
771: 0.0196333029355,
772: 0.120791820134,
773: 0.705002209888,
774: 0.174656645238,
775: 0.0221843026408,
776: 0.730165028048,
777: 0.0729940240731,
778: 0.12303063253,
779: 0.0156269392579,
780: 0.514338662174,
781: 0.347003766038,
782: 0.0555803583921,
783: 0.00705963113955,
784: 0.23235780058,
785: 0.705002209888,
786: 0.422442204032,
787: 0.422442204032,
788: 0.0821215678634,
789: 0.0729940240731,
790: 0.297574315013,
791: 0.297574315013,
792: 0.0578476039361,
793: 0.347003766038,
794: 0.134432268912,
795: 0.134432268912,
796: 0.0261332522867,
797: 0.705002209888,
798: 0.273713872823,
799: 0.273713872823,
800: 0.379578230281,
801: 0.0729940240731,
802: 0.192807956775,
803: 0.192807956775,
804: 0.267380320412,
805: 0.347003766038,
806: 0.0871029849888,
807: 0.0871029849888,
808: 0.120791820134,
809: 0.705002209888,
810: 0.0984204739396,
811: 0.0984204739396,
812: 0.730165028048,
813: 0.0729940240731,
814: 0.0693287858938,
815: 0.0693287858938,
816: 0.514338662174,
817: 0.347003766038,
818: 0.0313199947658,
819: 0.0313199947658,
820: 0.23235780058,
821: 0.705002209888,
822: 0.0952198798417,
823: 0.749664528222,
824: 0.0821215678634,
825: 0.0729940240731,
826: 0.0670742417521,
827: 0.528074388273,
828: 0.0578476039361,
829: 0.347003766038,
830: 0.0303014811743,
831: 0.23856305665,
832: 0.0261332522867,
833: 0.705002209888,
834: 0.0616960186091,
835: 0.485731727037,
836: 0.379578230281,
837: 0.0729940240731,
838: 0.0434595556538,
839: 0.342156357896,
840: 0.267380320412,
841: 0.347003766038,
842: 0.0196333029355,
843: 0.154572667042,
844: 0.120791820134,
845: 0.705002209888,
846: 0.0221843026408,
847: 0.174656645238,
848: 0.730165028048,
849: 0.0729940240731,
850: 0.0156269392579,
851: 0.12303063253,
852: 0.514338662174,
853: 0.347003766038,
854: 0.00705963113955,
855: 0.0555803583921,
856: 0.23235780058,
857: 0.705002209888};
859: /* Nodal basis function derivative evaluations */
860: static double BasisDerivatives[324] = {
861: -0.5,
862: -0.5,
863: -0.5,
864: 0.5,
865: 2.6964953901e-17,
866: 8.15881875835e-17,
867: 0.0,
868: 0.5,
869: 1.52600313755e-17,
870: 0.0,
871: 0.0,
872: 0.5,
873: -0.5,
874: -0.5,
875: -0.5,
876: 0.5,
877: 2.6938349868e-17,
878: 1.08228622783e-16,
879: 0.0,
880: 0.5,
881: -1.68114298577e-17,
882: 0.0,
883: 0.0,
884: 0.5,
885: -0.5,
886: -0.5,
887: -0.5,
888: 0.5,
889: 2.69035912409e-17,
890: 1.43034809879e-16,
891: 0.0,
892: 0.5,
893: -5.87133459397e-17,
894: 0.0,
895: 0.0,
896: 0.5,
897: -0.5,
898: -0.5,
899: -0.5,
900: 0.5,
901: 2.6964953901e-17,
902: 2.8079494593e-17,
903: 0.0,
904: 0.5,
905: 1.52600313755e-17,
906: 0.0,
907: 0.0,
908: 0.5,
909: -0.5,
910: -0.5,
911: -0.5,
912: 0.5,
913: 2.6938349868e-17,
914: 7.0536336094e-17,
915: 0.0,
916: 0.5,
917: -1.68114298577e-17,
918: 0.0,
919: 0.0,
920: 0.5,
921: -0.5,
922: -0.5,
923: -0.5,
924: 0.5,
925: 2.69035912409e-17,
926: 1.26006930238e-16,
927: 0.0,
928: 0.5,
929: -5.87133459397e-17,
930: 0.0,
931: 0.0,
932: 0.5,
933: -0.5,
934: -0.5,
935: -0.5,
936: 0.5,
937: 2.6964953901e-17,
938: -3.49866380524e-17,
939: 0.0,
940: 0.5,
941: 1.52600313755e-17,
942: 0.0,
943: 0.0,
944: 0.5,
945: -0.5,
946: -0.5,
947: -0.5,
948: 0.5,
949: 2.6938349868e-17,
950: 2.61116525673e-17,
951: 0.0,
952: 0.5,
953: -1.68114298577e-17,
954: 0.0,
955: 0.0,
956: 0.5,
957: -0.5,
958: -0.5,
959: -0.5,
960: 0.5,
961: 2.69035912409e-17,
962: 1.05937620823e-16,
963: 0.0,
964: 0.5,
965: -5.87133459397e-17,
966: 0.0,
967: 0.0,
968: 0.5,
969: -0.5,
970: -0.5,
971: -0.5,
972: 0.5,
973: 2.6964953901e-17,
974: 1.52807111565e-17,
975: 0.0,
976: 0.5,
977: 1.52600313755e-17,
978: 0.0,
979: 0.0,
980: 0.5,
981: -0.5,
982: -0.5,
983: -0.5,
984: 0.5,
985: 2.6938349868e-17,
986: 6.1520690456e-17,
987: 0.0,
988: 0.5,
989: -1.68114298577e-17,
990: 0.0,
991: 0.0,
992: 0.5,
993: -0.5,
994: -0.5,
995: -0.5,
996: 0.5,
997: 2.69035912409e-17,
998: 1.21934019246e-16,
999: 0.0,
1000: 0.5,
1001: -5.87133459397e-17,
1002: 0.0,
1003: 0.0,
1004: 0.5,
1005: -0.5,
1006: -0.5,
1007: -0.5,
1008: 0.5,
1009: 2.6964953901e-17,
1010: -1.48832491782e-17,
1011: 0.0,
1012: 0.5,
1013: 1.52600313755e-17,
1014: 0.0,
1015: 0.0,
1016: 0.5,
1017: -0.5,
1018: -0.5,
1019: -0.5,
1020: 0.5,
1021: 2.6938349868e-17,
1022: 4.0272766482e-17,
1023: 0.0,
1024: 0.5,
1025: -1.68114298577e-17,
1026: 0.0,
1027: 0.0,
1028: 0.5,
1029: -0.5,
1030: -0.5,
1031: -0.5,
1032: 0.5,
1033: 2.69035912409e-17,
1034: 1.1233505023e-16,
1035: 0.0,
1036: 0.5,
1037: -5.87133459397e-17,
1038: 0.0,
1039: 0.0,
1040: 0.5,
1041: -0.5,
1042: -0.5,
1043: -0.5,
1044: 0.5,
1045: 2.6964953901e-17,
1046: -5.04349365259e-17,
1047: 0.0,
1048: 0.5,
1049: 1.52600313755e-17,
1050: 0.0,
1051: 0.0,
1052: 0.5,
1053: -0.5,
1054: -0.5,
1055: -0.5,
1056: 0.5,
1057: 2.6938349868e-17,
1058: 1.52296507396e-17,
1059: 0.0,
1060: 0.5,
1061: -1.68114298577e-17,
1062: 0.0,
1063: 0.0,
1064: 0.5,
1065: -0.5,
1066: -0.5,
1067: -0.5,
1068: 0.5,
1069: 2.69035912409e-17,
1070: 1.01021564153e-16,
1071: 0.0,
1072: 0.5,
1073: -5.87133459397e-17,
1074: 0.0,
1075: 0.0,
1076: 0.5,
1077: -0.5,
1078: -0.5,
1079: -0.5,
1080: 0.5,
1081: 2.6964953901e-17,
1082: -5.10267652705e-17,
1083: 0.0,
1084: 0.5,
1085: 1.52600313755e-17,
1086: 0.0,
1087: 0.0,
1088: 0.5,
1089: -0.5,
1090: -0.5,
1091: -0.5,
1092: 0.5,
1093: 2.6938349868e-17,
1094: 1.4812758129e-17,
1095: 0.0,
1096: 0.5,
1097: -1.68114298577e-17,
1098: 0.0,
1099: 0.0,
1100: 0.5,
1101: -0.5,
1102: -0.5,
1103: -0.5,
1104: 0.5,
1105: 2.69035912409e-17,
1106: 1.00833228612e-16,
1107: 0.0,
1108: 0.5,
1109: -5.87133459397e-17,
1110: 0.0,
1111: 0.0,
1112: 0.5,
1113: -0.5,
1114: -0.5,
1115: -0.5,
1116: 0.5,
1117: 2.6964953901e-17,
1118: -5.78459929494e-17,
1119: 0.0,
1120: 0.5,
1121: 1.52600313755e-17,
1122: 0.0,
1123: 0.0,
1124: 0.5,
1125: -0.5,
1126: -0.5,
1127: -0.5,
1128: 0.5,
1129: 2.6938349868e-17,
1130: 1.00091968699e-17,
1131: 0.0,
1132: 0.5,
1133: -1.68114298577e-17,
1134: 0.0,
1135: 0.0,
1136: 0.5,
1137: -0.5,
1138: -0.5,
1139: -0.5,
1140: 0.5,
1141: 2.69035912409e-17,
1142: 9.86631702223e-17,
1143: 0.0,
1144: 0.5,
1145: -5.87133459397e-17,
1146: 0.0,
1147: 0.0,
1148: 0.5,
1149: -0.5,
1150: -0.5,
1151: -0.5,
1152: 0.5,
1153: 2.6964953901e-17,
1154: -6.58832349994e-17,
1155: 0.0,
1156: 0.5,
1157: 1.52600313755e-17,
1158: 0.0,
1159: 0.0,
1160: 0.5,
1161: -0.5,
1162: -0.5,
1163: -0.5,
1164: 0.5,
1165: 2.6938349868e-17,
1166: 4.34764891191e-18,
1167: 0.0,
1168: 0.5,
1169: -1.68114298577e-17,
1170: 0.0,
1171: 0.0,
1172: 0.5,
1173: -0.5,
1174: -0.5,
1175: -0.5,
1176: 0.5,
1177: 2.69035912409e-17,
1178: 9.61055074835e-17,
1179: 0.0,
1180: 0.5,
1181: -5.87133459397e-17,
1182: 0.0,
1183: 0.0,
1184: 0.5};
1186: #endif
1190: PetscErrorCode ElementGeometry(ALE::Obj<ALE::Two::Mesh> mesh, const ALE::Two::Mesh::point_type& e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
1191: {
1192: const double *coords = mesh->getCoordinates()->restrict(std::string("element"), e);
1193: int dim = mesh->getDimension();
1194: PetscReal det, invDet;
1197: if (debug) {
1198: MPI_Comm comm = mesh->comm();
1199: int rank = mesh->commRank();
1201: PetscSynchronizedPrintf(comm, "[%d]Element (%d, %d)\n", rank, e.prefix, e.index);
1202: PetscSynchronizedPrintf(comm, "[%d]Coordinates:\n[%d] ", rank, rank);
1203: for(int f = 0; f <= dim; f++) {
1204: PetscSynchronizedPrintf(comm, " (");
1205: for(int d = 0; d < dim; d++) {
1206: if (d > 0) PetscSynchronizedPrintf(comm, ", ");
1207: PetscSynchronizedPrintf(comm, "%g", coords[f*dim+d]);
1208: }
1209: PetscSynchronizedPrintf(comm, ")");
1210: }
1211: PetscSynchronizedPrintf(comm, "\n");
1212: }
1213: if (v0) {
1214: for(int d = 0; d < dim; d++) {
1215: v0[d] = coords[d];
1216: }
1217: }
1218: if (J) {
1219: for(int d = 0; d < dim; d++) {
1220: for(int f = 0; f < dim; f++) {
1221: J[d*dim+f] = 0.5*(coords[(f+1)*dim+d] - coords[0*dim+d]);
1222: }
1223: }
1224: if (debug) {
1225: MPI_Comm comm = mesh->comm();
1226: int rank = mesh->commRank();
1228: for(int d = 0; d < dim; d++) {
1229: if (d == 0) {
1230: PetscSynchronizedPrintf(comm, "[%d]J = /", rank);
1231: } else if (d == dim-1) {
1232: PetscSynchronizedPrintf(comm, "[%d] \\", rank);
1233: } else {
1234: PetscSynchronizedPrintf(comm, "[%d] |", rank);
1235: }
1236: for(int e = 0; e < dim; e++) {
1237: PetscSynchronizedPrintf(comm, " %g", J[d*dim+e]);
1238: }
1239: if (d == 0) {
1240: PetscSynchronizedPrintf(comm, " \\\n");
1241: } else if (d == dim-1) {
1242: PetscSynchronizedPrintf(comm, " /\n");
1243: } else {
1244: PetscSynchronizedPrintf(comm, " |\n");
1245: }
1246: }
1247: }
1248: if (dim == 2) {
1249: det = J[0]*J[3] - J[1]*J[2];
1250: } else if (dim == 3) {
1251: det = J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
1252: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
1253: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
1254: }
1255: invDet = 1.0/det;
1256: if (detJ) {
1257: if (det < 0) {SETERRQ(PETSC_ERR_ARG_WRONG, "Negative Matrix determinant");}
1258: *detJ = det;
1259: }
1260: if (invJ) {
1261: if (dim == 2) {
1262: invJ[0] = invDet*J[3];
1263: invJ[1] = -invDet*J[1];
1264: invJ[2] = -invDet*J[2];
1265: invJ[3] = invDet*J[0];
1266: } else if (dim == 3) {
1267: // FIX: This may be wrong
1268: invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
1269: invJ[0*3+1] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
1270: invJ[0*3+2] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
1271: invJ[1*3+0] = invDet*(J[0*3+1]*J[2*3+2] - J[0*3+2]*J[2*3+1]);
1272: invJ[1*3+1] = invDet*(J[0*3+2]*J[2*3+0] - J[0*3+0]*J[2*3+2]);
1273: invJ[1*3+2] = invDet*(J[0*3+0]*J[2*3+1] - J[0*3+1]*J[2*3+0]);
1274: invJ[2*3+0] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
1275: invJ[2*3+1] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
1276: invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
1277: }
1278: if (debug) {
1279: MPI_Comm comm = mesh->comm();
1280: int rank = mesh->commRank();
1282: for(int d = 0; d < dim; d++) {
1283: if (d == 0) {
1284: PetscSynchronizedPrintf(comm, "[%d]Jinv = /", rank);
1285: } else if (d == dim-1) {
1286: PetscSynchronizedPrintf(comm, "[%d] \\", rank);
1287: } else {
1288: PetscSynchronizedPrintf(comm, "[%d] |", rank);
1289: }
1290: for(int e = 0; e < dim; e++) {
1291: PetscSynchronizedPrintf(comm, " %g", invJ[d*dim+e]);
1292: }
1293: if (d == 0) {
1294: PetscSynchronizedPrintf(comm, " \\\n");
1295: } else if (d == dim-1) {
1296: PetscSynchronizedPrintf(comm, " /\n");
1297: } else {
1298: PetscSynchronizedPrintf(comm, " |\n");
1299: }
1300: }
1301: }
1302: }
1303: }
1304: return(0);
1305: }
1309: PetscErrorCode CheckElementGeometry(ALE::Obj<ALE::Two::Mesh> mesh)
1310: {
1311: ALE::Obj<ALE::Two::Mesh::sieve_type::traits::heightSequence> elements = mesh->getTopology()->heightStratum(0);
1312: PetscInt dim = mesh->getDimension();
1313: PetscReal *v0, *Jac;
1314: PetscReal detJ;
1318: PetscMalloc2(dim,PetscReal,&v0,dim*dim,PetscReal,&Jac);
1319: for(ALE::Two::Mesh::sieve_type::traits::heightSequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
1320: ElementGeometry(mesh, *e_iter, v0, Jac, PETSC_NULL, &detJ);
1321: }
1322: PetscSynchronizedFlush(mesh->comm());
1323: PetscFree2(v0,Jac);
1324: return(0);
1325: }
1329: PetscErrorCode ComputeRHS(DMMG dmmg, Vec b)
1330: {
1331: ALE::Obj<ALE::Two::Mesh> m;
1332: Mesh mesh = (Mesh) dmmg->dm;
1333: UserContext *user = (UserContext *) dmmg->user;
1334: MPI_Comm comm;
1335: PetscReal elementVec[NUM_BASIS_FUNCTIONS];
1336: PetscReal *v0, *Jac;
1337: PetscReal xi, eta, x_q, y_q, detJ, funcValue;
1338: PetscInt dim;
1339: PetscInt f, q;
1340: PetscErrorCode ierr;
1343: PetscObjectGetComm((PetscObject) mesh, &comm);
1344: MeshGetMesh(mesh, &m);
1345: dim = m->getDimension();
1346: PetscMalloc(dim * sizeof(PetscReal), &v0);
1347: PetscMalloc(dim*dim * sizeof(PetscReal), &Jac);
1348: ALE::Obj<ALE::Two::Mesh::field_type> field = m->getField("b");
1349: ALE::Obj<ALE::Two::Mesh::sieve_type::traits::heightSequence> elements = m->getTopology()->heightStratum(0);
1350: ALE::Two::Mesh::field_type::patch_type patch;
1351: for(ALE::Two::Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); e_itor++) {
1352: ElementGeometry(m, *e_itor, v0, Jac, PETSC_NULL, &detJ);
1353: /* Element integral */
1354: PetscMemzero(elementVec, NUM_BASIS_FUNCTIONS*sizeof(PetscScalar));
1355: for(q = 0; q < NUM_QUADRATURE_POINTS; q++) {
1356: xi = points[q*2+0] + 1.0;
1357: eta = points[q*2+1] + 1.0;
1358: x_q = Jac[0]*xi + Jac[1]*eta + v0[0];
1359: y_q = Jac[2]*xi + Jac[3]*eta + v0[1];
1360: funcValue = 2.0;
1361: //funcValue = 2.0*y_q*(2.0 - y_q) + 2.0*x_q*(2.0 - x_q);
1362: for(f = 0; f < NUM_BASIS_FUNCTIONS; f++) {
1363: elementVec[f] += Basis[q*NUM_BASIS_FUNCTIONS+f]*funcValue*weights[q]*detJ;
1364: }
1365: }
1366: if (debug) {PetscSynchronizedPrintf(comm, "elementVec = [%g %g %g]\n", elementVec[0], elementVec[1], elementVec[2]);}
1367: /* Assembly */
1368: field->updateAdd("element", *e_itor, elementVec);
1369: if (debug) {PetscSynchronizedFlush(comm);}
1370: }
1371: PetscFree(v0);
1372: PetscFree(Jac);
1374: Vec locB;
1375: VecCreateSeqWithArray(PETSC_COMM_SELF, field->getSize(patch), field->restrict(patch), &locB);
1376: VecScatterBegin(user->injection, locB, b, ADD_VALUES, SCATTER_FORWARD);
1377: VecScatterEnd(user->injection, locB, b, ADD_VALUES, SCATTER_FORWARD);
1378: VecDestroy(locB);
1380: /* force right hand side to be consistent for singular matrix */
1381: /* note this is really a hack, normally the model would provide you with a consistent right handside */
1382: if (user->bcType == NEUMANN) {
1383: MatNullSpace nullspace;
1385: KSPGetNullSpace(dmmg->ksp,&nullspace);
1386: MatNullSpaceRemove(nullspace,b,PETSC_NULL);
1387: }
1388: if (user->bcType == DIRICHLET) {
1389: /* Zero out BC rows */
1390: ALE::Two::Mesh::field_type::patch_type patch;
1391: ALE::Two::Mesh::field_type::patch_type bdPatch(0, 1);
1392: ALE::Obj<ALE::Two::Mesh::field_type> boundary = m->getBoundary();
1393: ALE::Obj<ALE::Two::Mesh::field_type::order_type::coneSequence> cone = boundary->getPatch(bdPatch);
1394: PetscScalar *boundaryValues;
1395: PetscInt *boundaryIndices;
1396: PetscInt numBoundaryIndices = 0;
1397: PetscInt k = 0;
1399: for(ALE::Two::Mesh::field_type::order_type::coneSequence::iterator p = cone->begin(); p != cone->end(); ++p) {
1400: numBoundaryIndices += field->getGlobalOrder()->getIndex(patch, *p).index;
1401: }
1402: PetscMalloc2(numBoundaryIndices,PetscInt,&boundaryIndices,numBoundaryIndices,PetscScalar,&boundaryValues);
1403: for(ALE::Two::Mesh::field_type::order_type::coneSequence::iterator p = cone->begin(); p != cone->end(); ++p) {
1404: const ALE::Two::Mesh::field_type::index_type& idx = field->getGlobalOrder()->getIndex(patch, *p);
1405: const double *data = boundary->restrict(bdPatch, *p);
1407: for(int i = 0; i < idx.index; i++) {
1408: boundaryIndices[k] = idx.prefix + i;
1409: boundaryValues[k] = data[i];
1410: k++;
1411: }
1412: }
1413: if (debug) {
1414: boundary->view("Boundary for rhs conditions");
1415: for(int i = 0; i < numBoundaryIndices; i++) {
1416: PetscSynchronizedPrintf(comm, "[%d]boundaryIndices[%d] = %d\n", m->commRank(), i, boundaryIndices[i]);
1417: }
1418: }
1419: PetscSynchronizedFlush(comm);
1420: VecSetValues(b, numBoundaryIndices, boundaryIndices, boundaryValues, INSERT_VALUES);
1421: PetscFree2(boundaryIndices, boundaryValues);
1422: }
1423: if (debug) {
1424: PetscPrintf(comm, "Rhs vector:");
1425: VecView(b, PETSC_VIEWER_STDOUT_WORLD);
1426: }
1427: return(0);
1428: }
1432: PetscErrorCode ComputeMatrix(DMMG dmmg, Mat J, Mat jac)
1433: {
1434: ALE::Obj<ALE::Two::Mesh> m;
1435: Mesh mesh = (Mesh) dmmg->dm;
1436: UserContext *user = (UserContext *) dmmg->user;
1437: MPI_Comm comm;
1438: PetscReal elementMat[NUM_BASIS_FUNCTIONS*NUM_BASIS_FUNCTIONS];
1439: PetscReal *v0, *Jac, *Jinv, *t_der, *b_der;
1440: PetscReal xi, eta, x_q, y_q, detJ;
1441: PetscInt dim;
1442: PetscInt f, g, q;
1443: PetscMPIInt rank;
1444: PetscErrorCode ierr;
1447: PetscObjectGetComm((PetscObject) mesh, &comm);
1448: MPI_Comm_rank(comm, &rank);
1449: MeshGetMesh(mesh, &m);
1450: dim = m->getDimension();
1451: PetscMalloc(dim * sizeof(PetscReal), &v0);
1452: PetscMalloc(dim * sizeof(PetscReal), &t_der);
1453: PetscMalloc(dim * sizeof(PetscReal), &b_der);
1454: PetscMalloc(dim*dim * sizeof(PetscReal), &Jac);
1455: PetscMalloc(dim*dim * sizeof(PetscReal), &Jinv);
1456: ALE::Obj<ALE::Two::Mesh::field_type> field = m->getField("u");
1457: ALE::Obj<ALE::Two::Mesh::sieve_type::traits::heightSequence> elements = m->getTopology()->heightStratum(0);
1458: for(ALE::Two::Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); e_itor++) {
1459: CHKMEMQ;
1460: ElementGeometry(m, *e_itor, v0, Jac, Jinv, &detJ);
1461: /* Element integral */
1462: PetscMemzero(elementMat, NUM_BASIS_FUNCTIONS*NUM_BASIS_FUNCTIONS*sizeof(PetscScalar));
1463: for(q = 0; q < NUM_QUADRATURE_POINTS; q++) {
1464: xi = points[q*2+0] + 1.0;
1465: eta = points[q*2+1] + 1.0;
1466: x_q = Jac[0]*xi + Jac[1]*eta + v0[0];
1467: y_q = Jac[2]*xi + Jac[3]*eta + v0[1];
1468: for(f = 0; f < NUM_BASIS_FUNCTIONS; f++) {
1469: t_der[0] = Jinv[0]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+0] + Jinv[2]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+1];
1470: t_der[1] = Jinv[1]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+0] + Jinv[3]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+1];
1471: for(g = 0; g < NUM_BASIS_FUNCTIONS; g++) {
1472: b_der[0] = Jinv[0]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+0] + Jinv[2]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+1];
1473: b_der[1] = Jinv[1]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+0] + Jinv[3]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+1];
1474: elementMat[f*NUM_BASIS_FUNCTIONS+g] += (t_der[0]*b_der[0] + t_der[1]*b_der[1])*weights[q]*detJ;
1475: }
1476: }
1477: }
1478: if (debug) {
1479: PetscSynchronizedPrintf(comm, "[%d]elementMat = [%g %g %g]\n [%g %g %g]\n [%g %g %g]\n",
1480: rank, elementMat[0], elementMat[1], elementMat[2], elementMat[3], elementMat[4],
1481: elementMat[5], elementMat[6], elementMat[7], elementMat[8]);
1482: }
1483: /* Assembly */
1484: updateOperator(jac, field, *e_itor, elementMat, ADD_VALUES);
1485: if (debug) {PetscSynchronizedFlush(comm);}
1486: }
1487: PetscFree(v0);
1488: PetscFree(t_der);
1489: PetscFree(b_der);
1490: PetscFree(Jac);
1491: PetscFree(Jinv);
1492: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
1493: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
1495: if (user->bcType == DIRICHLET) {
1496: /* Zero out BC rows */
1497: ALE::Two::Mesh::field_type::patch_type patch;
1498: ALE::Two::Mesh::field_type::patch_type bdPatch(0, 1);
1499: ALE::Obj<ALE::Two::Mesh::field_type> boundary = m->getBoundary();
1500: ALE::Obj<ALE::Two::Mesh::field_type::order_type::coneSequence> cone = boundary->getPatch(bdPatch);
1501: PetscInt *boundaryIndices;
1502: PetscInt numBoundaryIndices = 0;
1503: PetscInt k = 0;
1505: for(ALE::Two::Mesh::field_type::order_type::coneSequence::iterator p = cone->begin(); p != cone->end(); ++p) {
1506: numBoundaryIndices += field->getGlobalOrder()->getIndex(patch, *p).index;
1507: }
1508: PetscMalloc(numBoundaryIndices * sizeof(PetscInt), &boundaryIndices);
1509: for(ALE::Two::Mesh::field_type::order_type::coneSequence::iterator p = cone->begin(); p != cone->end(); ++p) {
1510: const ALE::Two::Mesh::field_type::index_type& idx = field->getGlobalOrder()->getIndex(patch, *p);
1512: for(int i = 0; i < idx.index; i++) {
1513: boundaryIndices[k++] = idx.prefix + i;
1514: }
1515: }
1516: if (debug) {
1517: for(int i = 0; i < numBoundaryIndices; i++) {
1518: PetscSynchronizedPrintf(comm, "[%d]boundaryIndices[%d] = %d\n", rank, i, boundaryIndices[i]);
1519: }
1520: }
1521: PetscSynchronizedFlush(comm);
1522: MatZeroRows(jac, numBoundaryIndices, boundaryIndices, 1.0);
1523: PetscFree(boundaryIndices);
1524: }
1525: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
1526: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
1527: return(0);
1528: }
1532: PetscErrorCode ComputeError(DMMG dmmg, Vec u, double *error)
1533: {
1534: ALE::Obj<ALE::Two::Mesh> m;
1535: Mesh mesh = (Mesh) dmmg->dm;
1536: UserContext *user = (UserContext *) dmmg->user;
1537: MPI_Comm comm;
1538: PetscReal *v0, *Jac;
1539: PetscReal detJ, totalError;
1540: PetscInt dim;
1541: PetscInt f, q;
1542: PetscErrorCode ierr;
1545: PetscObjectGetComm((PetscObject) mesh, &comm);
1546: MeshGetMesh(mesh, &m);
1547: dim = m->getDimension();
1548: PetscMalloc(dim * sizeof(PetscReal), &v0);
1549: PetscMalloc(dim*dim * sizeof(PetscReal), &Jac);
1550: ALE::Obj<ALE::Two::Mesh::field_type> field = m->getField("u");
1551: ALE::Obj<ALE::Two::Mesh::sieve_type::traits::heightSequence> elements = m->getTopology()->heightStratum(0);
1552: ALE::Two::Mesh::field_type::patch_type patch;
1554: Vec locU;
1555: VecCreateSeqWithArray(PETSC_COMM_SELF, field->getSize(patch), field->restrict(patch), &locU);
1556: VecScatterBegin(user->injection, u, locU, INSERT_VALUES, SCATTER_REVERSE);
1557: VecScatterEnd(user->injection, u, locU, INSERT_VALUES, SCATTER_REVERSE);
1558: VecDestroy(locU);
1560: totalError = 0.0;
1561: for(ALE::Two::Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); e_itor++) {
1562: const double *elementVec = field->restrict("element", *e_itor);
1563: double e;
1565: e = 0.0;
1566: ElementGeometry(m, *e_itor, v0, Jac, PETSC_NULL, &detJ);
1567: /* Element integral */
1568: for(q = 0; q < NUM_QUADRATURE_POINTS; q++) {
1569: PetscReal xi, eta, x_q, y_q, uExact, uComputed;
1571: xi = points[q*2+0] + 1.0;
1572: eta = points[q*2+1] + 1.0;
1573: x_q = Jac[0]*xi + Jac[1]*eta + v0[0];
1574: y_q = Jac[2]*xi + Jac[3]*eta + v0[1];
1575: uExact = x_q*(2.0 - x_q);
1576: //uExact = x_q*(2.0 - x_q)*y_q*(2.0 - y_q);
1577: uComputed = 0.0;
1578: for(f = 0; f < NUM_BASIS_FUNCTIONS; f++) {
1579: uComputed += elementVec[f]*Basis[q*NUM_BASIS_FUNCTIONS+f];
1580: }
1581: e += (uComputed - uExact)*(uComputed - uExact)*weights[q]*detJ;
1582: }
1583: totalError += e;
1584: if (debug) {PetscSynchronizedPrintf(comm, "elementError = %g\n", e);}
1585: if (debug) {PetscSynchronizedFlush(comm);}
1586: }
1587: PetscFree(v0);
1588: PetscFree(Jac);
1589: *error = sqrt(totalError);
1590: return(0);
1591: }