Actual source code: ex35.c
1: /*T
2: Concepts: KSP^solving a system of linear equations
3: Concepts: KSP^Laplacian, 2d
4: Processors: n
5: T*/
7: /*
10: Subduction Zone Benchmark
12: */
14: static char help[] = "Subduction Zone Benchmark\n\n";
16: #include petscmesh.h
17: #include petscksp.h
18: #include petscmg.h
19: #include petscdmmg.h
21: PetscErrorCode MeshView_Sieve_Newer(ALE::Obj<ALE::Two::Mesh>, PetscViewer);
22: PetscErrorCode CreateZoneBoundary(ALE::Obj<ALE::Two::Mesh>);
24: PetscErrorCode updateOperator(Mat, ALE::Obj<ALE::Two::Mesh::field_type>, const ALE::Two::Mesh::point_type&, PetscScalar [], InsertMode);
29: typedef enum {DIRICHLET, NEUMANN} BCType;
31: typedef struct {
32: PetscScalar nu;
33: BCType bcType;
34: VecScatter injection;
35: } UserContext;
37: PetscInt debug;
41: int main(int argc,char **argv)
42: {
43: MPI_Comm comm;
44: DMMG *dmmg;
45: UserContext user;
46: PetscViewer viewer;
47: const char *bcTypes[2] = {"dirichlet", "neumann"};
48: PetscReal refinementLimit, norm;
49: PetscInt dim, bc, l;
52: PetscInitialize(&argc,&argv,(char *)0,help);
53: comm = PETSC_COMM_WORLD;
55: PetscOptionsBegin(comm, "", "Options for the inhomogeneous Poisson equation", "DMMG");
56: debug = 0;
57: PetscOptionsInt("-debug", "The debugging flag", "ex33.c", 0, &debug, PETSC_NULL);
58: dim = 2;
59: PetscOptionsInt("-dim", "The mesh dimension", "ex33.c", 2, &dim, PETSC_NULL);
60: refinementLimit = 0.0;
61: PetscOptionsReal("-refinement_limit", "The area of the largest triangle in the mesh", "ex33.c", 1.0, &refinementLimit, PETSC_NULL);
62: user.nu = 0.1;
63: PetscOptionsScalar("-nu", "The width of the Gaussian source", "ex33.c", 0.1, &user.nu, PETSC_NULL);
64: bc = (PetscInt)DIRICHLET;
65: PetscOptionsEList("-bc_type","Type of boundary condition","ex33.c",bcTypes,2,bcTypes[0],&bc,PETSC_NULL);
66: user.bcType = (BCType) bc;
67: PetscOptionsEnd();
70: ALE::Obj<ALE::Two::Mesh> meshBoundary = ALE::Two::Mesh(comm, 1, debug);
71: ALE::Obj<ALE::Two::Mesh> mesh;
73: try {
74: ALE::LogStage stage = ALE::LogStageRegister("MeshCreation");
75: ALE::LogStagePush(stage);
76: PetscPrintf(comm, "Generating mesh\n");
77: CreateZoneBoundary(meshBoundary);
78: mesh = ALE::Two::Generator::generate(meshBoundary);
79: ALE::Obj<ALE::Two::Mesh::sieve_type> topology = mesh->getTopology();
80: PetscPrintf(comm, " Read %d elements\n", topology->heightStratum(0)->size());
81: PetscPrintf(comm, " Read %d vertices\n", topology->depthStratum(0)->size());
82: ALE::LogStagePop(stage);
84: stage = ALE::LogStageRegister("BndCreation");
85: ALE::LogStagePush(stage);
86: PetscPrintf(comm, "Creating boundary\n");
87: ALE::Obj<ALE::Two::Mesh::field_type> boundary = mesh->getBoundary();
88: ALE::Obj<ALE::Two::Mesh::sieve_type::traits::depthSequence> bdVertices = topology->depthStratum(0, 1);
89: ALE::Two::Mesh::field_type::patch_type patch;
91: boundary->setTopology(topology);
92: boundary->setPatch(topology->leaves(), patch);
93: //boundary->setFiberDimensionByDepth(patch, 0, 1);
94: for(ALE::Two::Mesh::sieve_type::traits::depthSequence::iterator v_iter = bdVertices->begin(); v_iter != bdVertices->end(); ++v_iter) {
95: boundary->setFiberDimension(patch, *v_iter, 1);
96: }
97: boundary->orderPatches();
98: for(ALE::Two::Mesh::sieve_type::traits::depthSequence::iterator v_iter = bdVertices->begin(); v_iter != bdVertices->end(); ++v_iter) {
99: //double *coords = mesh->getCoordinates()->restrict(patch, *v_iter);
100: double values[1] = {0.0};
102: boundary->update(patch, *v_iter, values);
103: }
104: boundary->view("Mesh Boundary");
105: ALE::LogStagePop(stage);
107: stage = ALE::LogStageRegister("MeshDistribution");
108: ALE::LogStagePush(stage);
109: PetscPrintf(comm, "Distributing mesh\n");
110: mesh = mesh->distribute();
111: ALE::LogStagePop(stage);
112: mesh->getBoundary()->view("Mesh Boundary");
114: /* Refine... ideally we will refine based on the distance from the singularity.
115: This "function call" based refinement has not yet been implemented. Matt's doing it... */
116: if (refinementLimit > 0.0) {
117: stage = ALE::LogStageRegister("MeshRefine");
118: ALE::LogStagePush(stage);
119: PetscPrintf(comm, "Refining mesh\n");
120: mesh = ALE::Two::Generator::refine(mesh, refinementLimit);
121: ALE::LogStagePop(stage);
122: }
123: topology = mesh->getTopology();
125: ALE::Obj<ALE::Two::Mesh::field_type> coords = mesh->getCoordinates();
126: ALE::Obj<ALE::Two::Mesh::field_type> u = mesh->getField("u");
127: u->setPatch(topology->leaves(), ALE::Two::Mesh::field_type::patch_type());
128: u->setFiberDimensionByDepth(patch, 0, 1);
129: u->orderPatches();
130: u->createGlobalOrder();
131: ALE::Obj<ALE::Two::Mesh::sieve_type::traits::heightSequence> elements = topology->heightStratum(0);
132: ALE::Obj<ALE::Two::Mesh::bundle_type> vertexBundle = mesh->getBundle(0);
133: std::string orderName("element");
135: for(ALE::Two::Mesh::sieve_type::traits::heightSequence::iterator e_iter = elements->begin(); e_iter != elements->end(); e_iter++) {
136: // setFiberDimensionByDepth() does not work here since we only want it to apply to the patch cone
137: // What we really need is the depthStratum relative to the patch
138: ALE::Obj<ALE::Two::Mesh::bundle_type::order_type::coneSequence> cone = vertexBundle->getPatch(orderName, *e_iter);
140: coords->setPatch(orderName, cone, *e_iter);
141: u->setPatch(orderName, cone, *e_iter);
142: for(ALE::Two::Mesh::bundle_type::order_type::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
143: u->setFiberDimension(orderName, *e_iter, *c_iter, 1);
144: }
145: }
146: u->orderPatches(orderName);
148: Mesh petscMesh;
149: MeshCreate(comm, &petscMesh);
150: MeshSetMesh(petscMesh, mesh);
151: DMMGCreate(comm,3,PETSC_NULL,&dmmg);
152: DMMGSetDM(dmmg, (DM) petscMesh);
153: MeshDestroy(petscMesh);
154: for (l = 0; l < DMMGGetLevels(dmmg); l++) {
155: DMMGSetUser(dmmg,l,&user);
156: }
158: DMMGSetKSP(dmmg,ComputeRHS,ComputeMatrix);
159: if (user.bcType == NEUMANN) {
160: DMMGSetNullSpace(dmmg,PETSC_TRUE,0,PETSC_NULL);
161: }
162: MeshGetGlobalScatter(mesh, "u", DMMGGetx(dmmg), &user.injection);
164: DMMGSolve(dmmg);
166: MatMult(DMMGGetJ(dmmg),DMMGGetx(dmmg),DMMGGetr(dmmg));
167: VecAXPY(DMMGGetr(dmmg),-1.0,DMMGGetRHS(dmmg));
168: VecNorm(DMMGGetr(dmmg),NORM_2,&norm);
169: PetscPrintf(comm,"Residual norm %g\n",norm);
170: VecAssemblyBegin(DMMGGetx(dmmg));
171: VecAssemblyEnd(DMMGGetx(dmmg));
173: stage = ALE::LogStageRegister("MeshOutput");
174: ALE::LogStagePush(stage);
175: PetscPrintf(comm, "Creating VTK mesh file\n");
176: PetscViewerCreate(comm, &viewer);
177: PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);
178: PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_VTK);
179: PetscViewerFileSetName(viewer, "poisson.vtk");
180: MeshView_Sieve_Newer(mesh, viewer);
181: //VecView(DMMGGetRHS(dmmg), viewer);
182: VecView(DMMGGetx(dmmg), viewer);
183: PetscViewerDestroy(viewer);
184: ALE::LogStagePop(stage);
186: DMMGDestroy(dmmg);
187: } catch (ALE::Exception e) {
188: std::cout << e << std::endl;
189: }
190: PetscFinalize();
191: return 0;
192: }
196: /*
197: Subduction zone benchmark boundary
199: Vertices Edges
201: 0---------1 -----7-----
202: |\ | | 12 14
203: | 2-------3 | \---8----|
204: | \ | 11 \ |
205: | \ | | 13 15
206: | \ | | \ |
207: 4-----5---6 --9-----10-
210: */
212: PetscErrorCode CreateZoneBoundary(ALE::Obj<ALE::Two::Mesh> mesh)
213: {
214: MPI_Comm comm = mesh->comm();
215: ALE::Obj<ALE::Two::Mesh::sieve_type> topology = mesh->getTopology();
216: PetscScalar coords[14] = { 0., 0.,
217: 660., 0.,
218: 50., 50.,
219: 660., 50.,
220: 0.,600.,
221: 600.,600.,
222: 660.,600.};
223: PetscInt conns[18] = { 0, 1,
224: 1, 3,
225: 3, 6,
226: 6, 5,
227: 5, 4,
228: 4, 0,
229: 0, 2,
230: 2, 3,
231: 2, 5};
234: PetscInt order = 0;
235: PetscMPIInt rank;
236: PetscErrorCode ierr;
239: MPI_Comm_rank(comm, &rank);
240: mesh->populate(9, conns, 7, coords, false);
242: /* Create boundary condition markers
243: 1 - surface
244: 2 - left inlet boundary
245: 3 - right corner flow boundary
246: etc
247: */
248: if (rank == 0) {
249: topology->setMarker(ALE::Two::Mesh::point_type(0, 0), 1);
250: topology->setMarker(ALE::Two::Mesh::point_type(0, 1), 1);
251: topology->setMarker(ALE::Two::Mesh::point_type(0, 3), 1);
252: topology->setMarker(ALE::Two::Mesh::point_type(0, 4), 1);
253: topology->setMarker(ALE::Two::Mesh::point_type(0, 5), 1);
254: topology->setMarker(ALE::Two::Mesh::point_type(0, 6), 1);
256: topology->setMarker(ALE::Two::Mesh::point_type(0, 7), 1);
257: topology->setMarker(ALE::Two::Mesh::point_type(0, 8), 1);
258: topology->setMarker(ALE::Two::Mesh::point_type(0, 9), 1);
259: topology->setMarker(ALE::Two::Mesh::point_type(0, 10), 1);
260: topology->setMarker(ALE::Two::Mesh::point_type(0, 11), 1);
261: topology->setMarker(ALE::Two::Mesh::point_type(0, 12), 1);
262: }
263: return(0);
264: }
267: #ifndef MESH_3D
269: #define NUM_QUADRATURE_POINTS 9
271: /* Quadrature points */
272: static double points[18] = {
273: -0.794564690381,
274: -0.822824080975,
275: -0.866891864322,
276: -0.181066271119,
277: -0.952137735426,
278: 0.575318923522,
279: -0.0885879595127,
280: -0.822824080975,
281: -0.409466864441,
282: -0.181066271119,
283: -0.787659461761,
284: 0.575318923522,
285: 0.617388771355,
286: -0.822824080975,
287: 0.0479581354402,
288: -0.181066271119,
289: -0.623181188096,
290: 0.575318923522};
292: /* Quadrature weights */
293: static double weights[9] = {
294: 0.223257681932,
295: 0.2547123404,
296: 0.0775855332238,
297: 0.357212291091,
298: 0.407539744639,
299: 0.124136853158,
300: 0.223257681932,
301: 0.2547123404,
302: 0.0775855332238};
304: #define NUM_BASIS_FUNCTIONS 3
306: /* Nodal basis function evaluations */
307: static double Basis[27] = {
308: 0.808694385678,
309: 0.10271765481,
310: 0.0885879595127,
311: 0.52397906772,
312: 0.0665540678392,
313: 0.409466864441,
314: 0.188409405952,
315: 0.0239311322871,
316: 0.787659461761,
317: 0.455706020244,
318: 0.455706020244,
319: 0.0885879595127,
320: 0.29526656778,
321: 0.29526656778,
322: 0.409466864441,
323: 0.10617026912,
324: 0.10617026912,
325: 0.787659461761,
326: 0.10271765481,
327: 0.808694385678,
328: 0.0885879595127,
329: 0.0665540678392,
330: 0.52397906772,
331: 0.409466864441,
332: 0.0239311322871,
333: 0.188409405952,
334: 0.787659461761};
336: /* Nodal basis function derivative evaluations */
337: static double BasisDerivatives[54] = {
338: -0.5,
339: -0.5,
340: 0.5,
341: 4.74937635818e-17,
342: 0.0,
343: 0.5,
344: -0.5,
345: -0.5,
346: 0.5,
347: 4.74937635818e-17,
348: 0.0,
349: 0.5,
350: -0.5,
351: -0.5,
352: 0.5,
353: 4.74937635818e-17,
354: 0.0,
355: 0.5,
356: -0.5,
357: -0.5,
358: 0.5,
359: 4.74937635818e-17,
360: 0.0,
361: 0.5,
362: -0.5,
363: -0.5,
364: 0.5,
365: 4.74937635818e-17,
366: 0.0,
367: 0.5,
368: -0.5,
369: -0.5,
370: 0.5,
371: 4.74937635818e-17,
372: 0.0,
373: 0.5,
374: -0.5,
375: -0.5,
376: 0.5,
377: 4.74937635818e-17,
378: 0.0,
379: 0.5,
380: -0.5,
381: -0.5,
382: 0.5,
383: 4.74937635818e-17,
384: 0.0,
385: 0.5,
386: -0.5,
387: -0.5,
388: 0.5,
389: 4.74937635818e-17,
390: 0.0,
391: 0.5};
393: #else
395: #define NUM_QUADRATURE_POINTS 27
397: /* Quadrature points */
398: static double points[81] = {
399: -0.809560240317,
400: -0.835756864273,
401: -0.854011951854,
402: -0.865851516496,
403: -0.884304792128,
404: -0.305992467923,
405: -0.939397037651,
406: -0.947733495427,
407: 0.410004419777,
408: -0.876607962782,
409: -0.240843539439,
410: -0.854011951854,
411: -0.913080888692,
412: -0.465239359176,
413: -0.305992467923,
414: -0.960733394129,
415: -0.758416359732,
416: 0.410004419777,
417: -0.955631394718,
418: 0.460330056095,
419: -0.854011951854,
420: -0.968746121484,
421: 0.0286773243482,
422: -0.305992467923,
423: -0.985880737721,
424: -0.53528439884,
425: 0.410004419777,
426: -0.155115591937,
427: -0.835756864273,
428: -0.854011951854,
429: -0.404851369974,
430: -0.884304792128,
431: -0.305992467923,
432: -0.731135462175,
433: -0.947733495427,
434: 0.410004419777,
435: -0.452572254354,
436: -0.240843539439,
437: -0.854011951854,
438: -0.61438408645,
439: -0.465239359176,
440: -0.305992467923,
441: -0.825794030022,
442: -0.758416359732,
443: 0.410004419777,
444: -0.803159052121,
445: 0.460330056095,
446: -0.854011951854,
447: -0.861342428212,
448: 0.0286773243482,
449: -0.305992467923,
450: -0.937360010468,
451: -0.53528439884,
452: 0.410004419777,
453: 0.499329056443,
454: -0.835756864273,
455: -0.854011951854,
456: 0.0561487765469,
457: -0.884304792128,
458: -0.305992467923,
459: -0.522873886699,
460: -0.947733495427,
461: 0.410004419777,
462: -0.0285365459258,
463: -0.240843539439,
464: -0.854011951854,
465: -0.315687284208,
466: -0.465239359176,
467: -0.305992467923,
468: -0.690854665916,
469: -0.758416359732,
470: 0.410004419777,
471: -0.650686709523,
472: 0.460330056095,
473: -0.854011951854,
474: -0.753938734941,
475: 0.0286773243482,
476: -0.305992467923,
477: -0.888839283216,
478: -0.53528439884,
479: 0.410004419777};
481: /* Quadrature weights */
482: static double weights[27] = {
483: 0.0701637994372,
484: 0.0653012061324,
485: 0.0133734490519,
486: 0.0800491405774,
487: 0.0745014590358,
488: 0.0152576273199,
489: 0.0243830167241,
490: 0.022693189565,
491: 0.0046474825267,
492: 0.1122620791,
493: 0.104481929812,
494: 0.021397518483,
495: 0.128078624924,
496: 0.119202334457,
497: 0.0244122037118,
498: 0.0390128267586,
499: 0.0363091033041,
500: 0.00743597204272,
501: 0.0701637994372,
502: 0.0653012061324,
503: 0.0133734490519,
504: 0.0800491405774,
505: 0.0745014590358,
506: 0.0152576273199,
507: 0.0243830167241,
508: 0.022693189565,
509: 0.0046474825267};
511: #define NUM_BASIS_FUNCTIONS 4
513: /* Nodal basis function evaluations */
514: static double Basis[108] = {
515: 0.749664528222,
516: 0.0952198798417,
517: 0.0821215678634,
518: 0.0729940240731,
519: 0.528074388273,
520: 0.0670742417521,
521: 0.0578476039361,
522: 0.347003766038,
523: 0.23856305665,
524: 0.0303014811743,
525: 0.0261332522867,
526: 0.705002209888,
527: 0.485731727037,
528: 0.0616960186091,
529: 0.379578230281,
530: 0.0729940240731,
531: 0.342156357896,
532: 0.0434595556538,
533: 0.267380320412,
534: 0.347003766038,
535: 0.154572667042,
536: 0.0196333029355,
537: 0.120791820134,
538: 0.705002209888,
539: 0.174656645238,
540: 0.0221843026408,
541: 0.730165028048,
542: 0.0729940240731,
543: 0.12303063253,
544: 0.0156269392579,
545: 0.514338662174,
546: 0.347003766038,
547: 0.0555803583921,
548: 0.00705963113955,
549: 0.23235780058,
550: 0.705002209888,
551: 0.422442204032,
552: 0.422442204032,
553: 0.0821215678634,
554: 0.0729940240731,
555: 0.297574315013,
556: 0.297574315013,
557: 0.0578476039361,
558: 0.347003766038,
559: 0.134432268912,
560: 0.134432268912,
561: 0.0261332522867,
562: 0.705002209888,
563: 0.273713872823,
564: 0.273713872823,
565: 0.379578230281,
566: 0.0729940240731,
567: 0.192807956775,
568: 0.192807956775,
569: 0.267380320412,
570: 0.347003766038,
571: 0.0871029849888,
572: 0.0871029849888,
573: 0.120791820134,
574: 0.705002209888,
575: 0.0984204739396,
576: 0.0984204739396,
577: 0.730165028048,
578: 0.0729940240731,
579: 0.0693287858938,
580: 0.0693287858938,
581: 0.514338662174,
582: 0.347003766038,
583: 0.0313199947658,
584: 0.0313199947658,
585: 0.23235780058,
586: 0.705002209888,
587: 0.0952198798417,
588: 0.749664528222,
589: 0.0821215678634,
590: 0.0729940240731,
591: 0.0670742417521,
592: 0.528074388273,
593: 0.0578476039361,
594: 0.347003766038,
595: 0.0303014811743,
596: 0.23856305665,
597: 0.0261332522867,
598: 0.705002209888,
599: 0.0616960186091,
600: 0.485731727037,
601: 0.379578230281,
602: 0.0729940240731,
603: 0.0434595556538,
604: 0.342156357896,
605: 0.267380320412,
606: 0.347003766038,
607: 0.0196333029355,
608: 0.154572667042,
609: 0.120791820134,
610: 0.705002209888,
611: 0.0221843026408,
612: 0.174656645238,
613: 0.730165028048,
614: 0.0729940240731,
615: 0.0156269392579,
616: 0.12303063253,
617: 0.514338662174,
618: 0.347003766038,
619: 0.00705963113955,
620: 0.0555803583921,
621: 0.23235780058,
622: 0.705002209888};
624: /* Nodal basis function derivative evaluations */
625: static double BasisDerivatives[324] = {
626: -0.5,
627: -0.5,
628: -0.5,
629: 0.5,
630: 2.6964953901e-17,
631: 8.15881875835e-17,
632: 0.0,
633: 0.5,
634: 1.52600313755e-17,
635: 0.0,
636: 0.0,
637: 0.5,
638: -0.5,
639: -0.5,
640: -0.5,
641: 0.5,
642: 2.6938349868e-17,
643: 1.08228622783e-16,
644: 0.0,
645: 0.5,
646: -1.68114298577e-17,
647: 0.0,
648: 0.0,
649: 0.5,
650: -0.5,
651: -0.5,
652: -0.5,
653: 0.5,
654: 2.69035912409e-17,
655: 1.43034809879e-16,
656: 0.0,
657: 0.5,
658: -5.87133459397e-17,
659: 0.0,
660: 0.0,
661: 0.5,
662: -0.5,
663: -0.5,
664: -0.5,
665: 0.5,
666: 2.6964953901e-17,
667: 2.8079494593e-17,
668: 0.0,
669: 0.5,
670: 1.52600313755e-17,
671: 0.0,
672: 0.0,
673: 0.5,
674: -0.5,
675: -0.5,
676: -0.5,
677: 0.5,
678: 2.6938349868e-17,
679: 7.0536336094e-17,
680: 0.0,
681: 0.5,
682: -1.68114298577e-17,
683: 0.0,
684: 0.0,
685: 0.5,
686: -0.5,
687: -0.5,
688: -0.5,
689: 0.5,
690: 2.69035912409e-17,
691: 1.26006930238e-16,
692: 0.0,
693: 0.5,
694: -5.87133459397e-17,
695: 0.0,
696: 0.0,
697: 0.5,
698: -0.5,
699: -0.5,
700: -0.5,
701: 0.5,
702: 2.6964953901e-17,
703: -3.49866380524e-17,
704: 0.0,
705: 0.5,
706: 1.52600313755e-17,
707: 0.0,
708: 0.0,
709: 0.5,
710: -0.5,
711: -0.5,
712: -0.5,
713: 0.5,
714: 2.6938349868e-17,
715: 2.61116525673e-17,
716: 0.0,
717: 0.5,
718: -1.68114298577e-17,
719: 0.0,
720: 0.0,
721: 0.5,
722: -0.5,
723: -0.5,
724: -0.5,
725: 0.5,
726: 2.69035912409e-17,
727: 1.05937620823e-16,
728: 0.0,
729: 0.5,
730: -5.87133459397e-17,
731: 0.0,
732: 0.0,
733: 0.5,
734: -0.5,
735: -0.5,
736: -0.5,
737: 0.5,
738: 2.6964953901e-17,
739: 1.52807111565e-17,
740: 0.0,
741: 0.5,
742: 1.52600313755e-17,
743: 0.0,
744: 0.0,
745: 0.5,
746: -0.5,
747: -0.5,
748: -0.5,
749: 0.5,
750: 2.6938349868e-17,
751: 6.1520690456e-17,
752: 0.0,
753: 0.5,
754: -1.68114298577e-17,
755: 0.0,
756: 0.0,
757: 0.5,
758: -0.5,
759: -0.5,
760: -0.5,
761: 0.5,
762: 2.69035912409e-17,
763: 1.21934019246e-16,
764: 0.0,
765: 0.5,
766: -5.87133459397e-17,
767: 0.0,
768: 0.0,
769: 0.5,
770: -0.5,
771: -0.5,
772: -0.5,
773: 0.5,
774: 2.6964953901e-17,
775: -1.48832491782e-17,
776: 0.0,
777: 0.5,
778: 1.52600313755e-17,
779: 0.0,
780: 0.0,
781: 0.5,
782: -0.5,
783: -0.5,
784: -0.5,
785: 0.5,
786: 2.6938349868e-17,
787: 4.0272766482e-17,
788: 0.0,
789: 0.5,
790: -1.68114298577e-17,
791: 0.0,
792: 0.0,
793: 0.5,
794: -0.5,
795: -0.5,
796: -0.5,
797: 0.5,
798: 2.69035912409e-17,
799: 1.1233505023e-16,
800: 0.0,
801: 0.5,
802: -5.87133459397e-17,
803: 0.0,
804: 0.0,
805: 0.5,
806: -0.5,
807: -0.5,
808: -0.5,
809: 0.5,
810: 2.6964953901e-17,
811: -5.04349365259e-17,
812: 0.0,
813: 0.5,
814: 1.52600313755e-17,
815: 0.0,
816: 0.0,
817: 0.5,
818: -0.5,
819: -0.5,
820: -0.5,
821: 0.5,
822: 2.6938349868e-17,
823: 1.52296507396e-17,
824: 0.0,
825: 0.5,
826: -1.68114298577e-17,
827: 0.0,
828: 0.0,
829: 0.5,
830: -0.5,
831: -0.5,
832: -0.5,
833: 0.5,
834: 2.69035912409e-17,
835: 1.01021564153e-16,
836: 0.0,
837: 0.5,
838: -5.87133459397e-17,
839: 0.0,
840: 0.0,
841: 0.5,
842: -0.5,
843: -0.5,
844: -0.5,
845: 0.5,
846: 2.6964953901e-17,
847: -5.10267652705e-17,
848: 0.0,
849: 0.5,
850: 1.52600313755e-17,
851: 0.0,
852: 0.0,
853: 0.5,
854: -0.5,
855: -0.5,
856: -0.5,
857: 0.5,
858: 2.6938349868e-17,
859: 1.4812758129e-17,
860: 0.0,
861: 0.5,
862: -1.68114298577e-17,
863: 0.0,
864: 0.0,
865: 0.5,
866: -0.5,
867: -0.5,
868: -0.5,
869: 0.5,
870: 2.69035912409e-17,
871: 1.00833228612e-16,
872: 0.0,
873: 0.5,
874: -5.87133459397e-17,
875: 0.0,
876: 0.0,
877: 0.5,
878: -0.5,
879: -0.5,
880: -0.5,
881: 0.5,
882: 2.6964953901e-17,
883: -5.78459929494e-17,
884: 0.0,
885: 0.5,
886: 1.52600313755e-17,
887: 0.0,
888: 0.0,
889: 0.5,
890: -0.5,
891: -0.5,
892: -0.5,
893: 0.5,
894: 2.6938349868e-17,
895: 1.00091968699e-17,
896: 0.0,
897: 0.5,
898: -1.68114298577e-17,
899: 0.0,
900: 0.0,
901: 0.5,
902: -0.5,
903: -0.5,
904: -0.5,
905: 0.5,
906: 2.69035912409e-17,
907: 9.86631702223e-17,
908: 0.0,
909: 0.5,
910: -5.87133459397e-17,
911: 0.0,
912: 0.0,
913: 0.5,
914: -0.5,
915: -0.5,
916: -0.5,
917: 0.5,
918: 2.6964953901e-17,
919: -6.58832349994e-17,
920: 0.0,
921: 0.5,
922: 1.52600313755e-17,
923: 0.0,
924: 0.0,
925: 0.5,
926: -0.5,
927: -0.5,
928: -0.5,
929: 0.5,
930: 2.6938349868e-17,
931: 4.34764891191e-18,
932: 0.0,
933: 0.5,
934: -1.68114298577e-17,
935: 0.0,
936: 0.0,
937: 0.5,
938: -0.5,
939: -0.5,
940: -0.5,
941: 0.5,
942: 2.69035912409e-17,
943: 9.61055074835e-17,
944: 0.0,
945: 0.5,
946: -5.87133459397e-17,
947: 0.0,
948: 0.0,
949: 0.5};
951: #endif
955: PetscErrorCode ElementGeometry(ALE::Obj<ALE::Two::Mesh> mesh, const ALE::Two::Mesh::point_type& e, PetscReal v0[], PetscReal J[], PetscReal invJ[], PetscReal *detJ)
956: {
957: const double *coords = mesh->getCoordinates()->restrict(std::string("element"), e);
958: int dim = mesh->getDimension();
959: PetscReal det, invDet;
963: if (v0) {
964: for(int d = 0; d < dim; d++) {
965: v0[d] = coords[d];
966: }
967: }
968: if (J) {
969: for(int d = 0; d < dim; d++) {
970: for(int f = 0; f < dim; f++) {
971: J[d*dim+f] = 0.5*(coords[(f+1)*dim+d] - coords[0*dim+d]);
972: }
973: }
974: if (debug) {
975: MPI_Comm comm = mesh->comm();
976: PetscMPIInt rank;
978: MPI_Comm_rank(comm, &rank);
979: for(int d = 0; d < dim; d++) {
980: if (d == 0) {
981: PetscSynchronizedPrintf(comm, "[%d]J = /", rank);
982: } else if (d == dim-1) {
983: PetscSynchronizedPrintf(comm, "[%d] \\", rank);
984: } else {
985: PetscSynchronizedPrintf(comm, "[%d] |", rank);
986: }
987: for(int e = 0; e < dim; e++) {
988: PetscSynchronizedPrintf(comm, " %g", J[d*dim+e]);
989: }
990: if (d == 0) {
991: PetscSynchronizedPrintf(comm, " \\\n");
992: } else if (d == dim-1) {
993: PetscSynchronizedPrintf(comm, " /\n");
994: } else {
995: PetscSynchronizedPrintf(comm, " |\n");
996: }
997: }
998: }
999: if (dim == 2) {
1000: det = J[0]*J[3] - J[1]*J[2];
1001: } else if (dim == 3) {
1002: det = J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
1003: J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
1004: J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
1005: }
1006: invDet = 1.0/det;
1007: if (detJ) {
1008: if (det < 0) {SETERRQ(PETSC_ERR_ARG_WRONG, "Negative Matrix determinant");}
1009: *detJ = det;
1010: }
1011: if (invJ) {
1012: if (dim == 2) {
1013: invJ[0] = invDet*J[3];
1014: invJ[1] = -invDet*J[1];
1015: invJ[2] = -invDet*J[2];
1016: invJ[3] = invDet*J[0];
1017: } else if (dim == 3) {
1018: // FIX: This may be wrong
1019: invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
1020: invJ[0*3+1] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
1021: invJ[0*3+2] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
1022: invJ[1*3+0] = invDet*(J[0*3+1]*J[2*3+2] - J[0*3+2]*J[2*3+1]);
1023: invJ[1*3+1] = invDet*(J[0*3+2]*J[2*3+0] - J[0*3+0]*J[2*3+2]);
1024: invJ[1*3+2] = invDet*(J[0*3+0]*J[2*3+1] - J[0*3+1]*J[2*3+0]);
1025: invJ[2*3+0] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
1026: invJ[2*3+1] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
1027: invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
1028: }
1029: if (debug) {
1030: MPI_Comm comm = mesh->comm();
1031: PetscMPIInt rank;
1033: MPI_Comm_rank(comm, &rank);
1034: for(int d = 0; d < dim; d++) {
1035: if (d == 0) {
1036: PetscSynchronizedPrintf(comm, "[%d]Jinv = /", rank);
1037: } else if (d == dim-1) {
1038: PetscSynchronizedPrintf(comm, "[%d] \\", rank);
1039: } else {
1040: PetscSynchronizedPrintf(comm, "[%d] |", rank);
1041: }
1042: for(int e = 0; e < dim; e++) {
1043: PetscSynchronizedPrintf(comm, " %g", invJ[d*dim+e]);
1044: }
1045: if (d == 0) {
1046: PetscSynchronizedPrintf(comm, " \\\n");
1047: } else if (d == dim-1) {
1048: PetscSynchronizedPrintf(comm, " /\n");
1049: } else {
1050: PetscSynchronizedPrintf(comm, " |\n");
1051: }
1052: }
1053: }
1054: }
1055: }
1056: return(0);
1057: }
1061: PetscErrorCode ComputeRho(PetscReal x, PetscReal y, PetscScalar *rho)
1062: {
1064: if ((x > 1.0/3.0) && (x < 2.0/3.0) && (y > 1.0/3.0) && (y < 2.0/3.0)) {
1065: //*rho = 100.0;
1066: *rho = 1.0;
1067: } else {
1068: *rho = 1.0;
1069: }
1070: return(0);
1071: }
1075: PetscErrorCode ComputeRHS(DMMG dmmg, Vec b)
1076: {
1077: ALE::Obj<ALE::Two::Mesh> m;
1078: Mesh mesh = (Mesh) dmmg->dm;
1079: UserContext *user = (UserContext *) dmmg->user;
1080: MPI_Comm comm;
1081: PetscReal elementVec[NUM_BASIS_FUNCTIONS];
1082: PetscReal *v0, *Jac;
1083: PetscReal xi, eta, x_q, y_q, detJ, funcValue;
1084: PetscInt dim;
1085: PetscInt f, q;
1086: PetscErrorCode ierr;
1089: PetscObjectGetComm((PetscObject) mesh, &comm);
1090: MeshGetMesh(mesh, &m);
1091: dim = m->getDimension();
1092: PetscMalloc(dim * sizeof(PetscReal), &v0);
1093: PetscMalloc(dim*dim * sizeof(PetscReal), &Jac);
1094: ALE::Obj<ALE::Two::Mesh::field_type> field = m->getField("u");
1095: ALE::Obj<ALE::Two::Mesh::sieve_type::traits::heightSequence> elements = m->getTopology()->heightStratum(0);
1096: ALE::Two::Mesh::field_type::patch_type patch;
1097: for(ALE::Two::Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); e_itor++) {
1098: ElementGeometry(m, *e_itor, v0, Jac, PETSC_NULL, &detJ);
1099: /* Element integral */
1100: PetscMemzero(elementVec, NUM_BASIS_FUNCTIONS*sizeof(PetscScalar));
1101: for(q = 0; q < NUM_QUADRATURE_POINTS; q++) {
1102: xi = points[q*2+0] + 1.0;
1103: eta = points[q*2+1] + 1.0;
1104: x_q = Jac[0]*xi + Jac[1]*eta + v0[0];
1105: y_q = Jac[2]*xi + Jac[3]*eta + v0[1];
1106: funcValue = PetscExpScalar(-(x_q*x_q)/user->nu)*PetscExpScalar(-(y_q*y_q)/user->nu);
1107: for(f = 0; f < NUM_BASIS_FUNCTIONS; f++) {
1108: elementVec[f] += Basis[q*NUM_BASIS_FUNCTIONS+f]*funcValue*weights[q]*detJ;
1109: }
1110: }
1111: if (debug) {PetscSynchronizedPrintf(comm, "elementVec = [%g %g %g]\n", elementVec[0], elementVec[1], elementVec[2]);}
1112: /* Assembly */
1113: field->updateAdd("element", *e_itor, elementVec);
1114: if (debug) {PetscSynchronizedFlush(comm);}
1115: }
1116: PetscFree(v0);
1117: PetscFree(Jac);
1119: Vec locB;
1120: VecCreateSeqWithArray(PETSC_COMM_SELF, field->getSize(patch), field->restrict(patch), &locB);
1121: VecScatterBegin(user->injection, locB, b, ADD_VALUES, SCATTER_FORWARD);
1122: VecScatterEnd(user->injection, locB, b, ADD_VALUES, SCATTER_FORWARD);
1123: VecDestroy(locB);
1125: /* force right hand side to be consistent for singular matrix */
1126: /* note this is really a hack, normally the model would provide you with a consistent right handside */
1127: if (user->bcType == NEUMANN) {
1128: MatNullSpace nullspace;
1130: KSPGetNullSpace(dmmg->ksp,&nullspace);
1131: MatNullSpaceRemove(nullspace,b,PETSC_NULL);
1132: }
1133: return(0);
1134: }
1138: PetscErrorCode ComputeMatrix(DMMG dmmg, Mat J, Mat jac)
1139: {
1140: ALE::Obj<ALE::Two::Mesh> m;
1141: Mesh mesh = (Mesh) dmmg->dm;
1142: UserContext *user = (UserContext *) dmmg->user;
1143: MPI_Comm comm;
1144: PetscReal elementMat[NUM_BASIS_FUNCTIONS*NUM_BASIS_FUNCTIONS];
1145: PetscReal *v0, *Jac, *Jinv, *t_der, *b_der;
1146: PetscReal xi, eta, x_q, y_q, detJ, rho;
1147: PetscInt dim;
1148: PetscInt f, g, q;
1149: PetscMPIInt rank;
1150: PetscErrorCode ierr;
1153: PetscObjectGetComm((PetscObject) mesh, &comm);
1154: MPI_Comm_rank(comm, &rank);
1155: MeshGetMesh(mesh, &m);
1156: dim = m->getDimension();
1157: PetscMalloc(dim * sizeof(PetscReal), &v0);
1158: PetscMalloc(dim * sizeof(PetscReal), &t_der);
1159: PetscMalloc(dim * sizeof(PetscReal), &b_der);
1160: PetscMalloc(dim*dim * sizeof(PetscReal), &Jac);
1161: PetscMalloc(dim*dim * sizeof(PetscReal), &Jinv);
1162: ALE::Obj<ALE::Two::Mesh::field_type> field = m->getField("u");
1163: ALE::Obj<ALE::Two::Mesh::sieve_type::traits::heightSequence> elements = m->getTopology()->heightStratum(0);
1164: for(ALE::Two::Mesh::sieve_type::traits::heightSequence::iterator e_itor = elements->begin(); e_itor != elements->end(); e_itor++) {
1165: CHKMEMQ;
1166: ElementGeometry(m, *e_itor, v0, Jac, Jinv, &detJ);
1167: /* Element integral */
1168: PetscMemzero(elementMat, NUM_BASIS_FUNCTIONS*NUM_BASIS_FUNCTIONS*sizeof(PetscScalar));
1169: for(q = 0; q < NUM_QUADRATURE_POINTS; q++) {
1170: xi = points[q*2+0] + 1.0;
1171: eta = points[q*2+1] + 1.0;
1172: x_q = Jac[0]*xi + Jac[1]*eta + v0[0];
1173: y_q = Jac[2]*xi + Jac[3]*eta + v0[1];
1174: ComputeRho(x_q, y_q, &rho);
1175: for(f = 0; f < NUM_BASIS_FUNCTIONS; f++) {
1176: t_der[0] = Jinv[0]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+0] + Jinv[2]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+1];
1177: t_der[1] = Jinv[1]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+0] + Jinv[3]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+f)*2+1];
1178: for(g = 0; g < NUM_BASIS_FUNCTIONS; g++) {
1179: b_der[0] = Jinv[0]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+0] + Jinv[2]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+1];
1180: b_der[1] = Jinv[1]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+0] + Jinv[3]*BasisDerivatives[(q*NUM_BASIS_FUNCTIONS+g)*2+1];
1181: elementMat[f*NUM_BASIS_FUNCTIONS+g] += rho*(t_der[0]*b_der[0] + t_der[1]*b_der[1])*weights[q]*detJ;
1182: }
1183: }
1184: }
1185: if (debug) {
1186: PetscSynchronizedPrintf(comm, "[%d]elementMat = [%g %g %g]\n [%g %g %g]\n [%g %g %g]\n",
1187: rank, elementMat[0], elementMat[1], elementMat[2], elementMat[3], elementMat[4],
1188: elementMat[5], elementMat[6], elementMat[7], elementMat[8]);
1189: }
1190: /* Assembly */
1191: updateOperator(jac, field, *e_itor, elementMat, ADD_VALUES);
1192: if (debug) {PetscSynchronizedFlush(comm);}
1193: }
1194: PetscFree(v0);
1195: PetscFree(t_der);
1196: PetscFree(b_der);
1197: PetscFree(Jac);
1198: PetscFree(Jinv);
1199: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
1200: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
1202: if (user->bcType == DIRICHLET) {
1203: /* Zero out BC rows */
1204: ALE::Obj<ALE::Two::Mesh::field_type> boundary = m->getBoundary();
1205: ALE::Obj<ALE::Two::Mesh::sieve_type::traits::depthSequence> vertices = m->getTopology()->depthStratum(0);
1206: ALE::Two::Mesh::field_type::patch_type patch;
1207: PetscInt *boundaryIndices;
1208: PetscInt numBoundaryIndices = 0;
1209: PetscInt k = 0;
1211: for(ALE::Two::Mesh::sieve_type::traits::depthSequence::iterator p = vertices->begin(); p != vertices->end(); ++p) {
1212: if (boundary->getIndex(patch, *p).index > 0) {
1213: const ALE::Two::Mesh::field_type::index_type& idx = field->getGlobalOrder()->getIndex(patch, *p);
1215: if (idx.index > 0) {
1216: numBoundaryIndices += idx.index;
1217: }
1218: }
1219: }
1220: PetscMalloc(numBoundaryIndices * sizeof(PetscInt), &boundaryIndices);
1221: for(ALE::Two::Mesh::sieve_type::traits::depthSequence::iterator p = vertices->begin(); p != vertices->end(); ++p) {
1222: if (boundary->getIndex(patch, *p).index > 0) {
1223: const ALE::Two::Mesh::field_type::index_type& idx = field->getGlobalOrder()->getIndex(patch, *p);
1225: for(int i = 0; i < idx.index; i++) {
1226: boundaryIndices[k++] = idx.prefix + i;
1227: }
1228: }
1229: }
1230: if (debug) {
1231: for(int i = 0; i < numBoundaryIndices; i++) {
1232: PetscSynchronizedPrintf(comm, "[%d]boundaryIndices[%d] = %d\n", rank, i, boundaryIndices[i]);
1233: }
1234: }
1235: PetscSynchronizedFlush(comm);
1236: MatZeroRows(jac, numBoundaryIndices, boundaryIndices, 1.0);
1237: PetscFree(boundaryIndices);
1238: }
1239: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
1240: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
1241: return(0);
1242: }