Actual source code: mathematica.c
1: #ifdef PETSC_RCS_HEADER
2: static char vcid[] = "$Id: mathematica.c,v 1.9 2000/01/26 15:46:22 baggag Exp $";
3: #endif
5: /*
6: Written by Matt Knepley, knepley@cs.purdue.edu 7/23/97
7: Major overhall for interactivity 11/14/97
8: Reorganized 11/8/98
9: */
10: #include src/sys/src/viewer/viewerimpl.h
11: #include src/gvec/gvecimpl.h
12: #include "src/mesh/impls/triangular/2d/2dimpl.h"
13: #include src/sles/pc/pcimpl.h
14: #include src/sles/pc/impls/ml/ml.h
15: #include src/mat/impls/aij/seq/aij.h
16: #include mathematica.h
17: #include "petscfix.h"
19: #if defined (PETSC_HAVE__SNPRINTF)
20: #define snprintf _snprintf
21: #endif
23: PetscViewer PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE = PETSC_NULL;
24: #ifdef PETSC_HAVE_MATHEMATICA
25: static void *mathematicaEnv = PETSC_NULL;
26: #endif
28: #undef __FUNCT__
30: /*@C
31: PetscViewerMathematicaInitializePackage - This function initializes everything in the Petsc interface to Mathematica. It is
32: called from PetscDLLibraryRegister() when using dynamic libraries, and on the call to PetscInitialize()
33: when using static libraries.
35: Input Parameter:
36: path - The dynamic library path, or PETSC_NULL
38: Level: developer
40: .keywords: Petsc, initialize, package, PLAPACK
41: .seealso: PetscInitializePackage(), PetscInitialize()
42: @*/
43: int PetscViewerMathematicaInitializePackage(char *path) {
44: static PetscTruth initialized = PETSC_FALSE;
47: if (initialized == PETSC_TRUE) return(0);
48: initialized = PETSC_TRUE;
49: #ifdef PETSC_HAVE_MATHEMATICA
50: mathematicaEnv = (void *) MLInitialize(0);
51: #endif
52: return(0);
53: }
55: #undef __FUNCT__
57: /*@C
58: PetscViewerMathematicaDestroyPackage - This function destroys everything in the Petsc interface to Mathematica. It is
59: called from PetscFinalize().
61: Level: developer
63: .keywords: Petsc, destroy, package, mathematica
64: .seealso: PetscFinalize()
65: @*/
66: int PetscViewerMathematicaFinalizePackage(void) {
68: #ifdef PETSC_HAVE_MATHEMATICA
69: if (mathematicaEnv != PETSC_NULL) MLDeinitialize((MLEnvironment) mathematicaEnv);
70: #endif
71: return(0);
72: }
74: #undef __FUNCT__
76: int PetscViewerInitializeMathematicaWorld_Private()
77: {
81: if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) return(0);
82: PetscViewerMathematicaOpen(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_NULL, PETSC_NULL, &PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);
83:
84: return(0);
85: }
87: #undef __FUNCT__
89: static int PetscViewerDestroy_Mathematica(PetscViewer viewer)
90: {
91: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
92: int ierr;
95: #ifdef PETSC_HAVE_MATHEMATICA
96: MLClose(vmath->link);
97: #endif
98: if (vmath->linkname != PETSC_NULL) {
99: PetscFree(vmath->linkname);
100: }
101: if (vmath->linkhost != PETSC_NULL) {
102: PetscFree(vmath->linkhost);
103: }
104: PetscFree(vmath);
105: return(0);
106: }
108: #undef __FUNCT__
110: int PetscViewerDestroyMathematica_Private(void)
111: {
115: if (PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE) {
116: PetscViewerDestroy(PETSC_VIEWER_MATHEMATICA_WORLD_PRIVATE);
117: }
118: return(0);
119: }
121: #undef __FUNCT__
123: int PetscViewerMathematicaSetupConnection_Private(PetscViewer v) {
124: #ifdef PETSC_HAVE_MATHEMATICA
125: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) v->data;
126: #ifdef MATHEMATICA_3_0
127: int argc = 6;
128: char *argv[6];
129: #else
130: int argc = 5;
131: char *argv[5];
132: #endif
133: char hostname[256];
134: long lerr;
135: int ierr;
136: #endif
139: #ifdef PETSC_HAVE_MATHEMATICA
140: /* Link name */
141: argv[0] = "-linkname";
142: if (vmath->linkname == PETSC_NULL) {
143: argv[1] = "math -mathlink";
144: } else {
145: argv[1] = vmath->linkname;
146: }
148: /* Link host */
149: argv[2] = "-linkhost";
150: if (vmath->linkhost == PETSC_NULL) {
151: PetscGetHostName(hostname, 255);
152: argv[3] = hostname;
153: } else {
154: argv[3] = vmath->linkhost;
155: }
157: /* Link mode */
158: #ifdef MATHEMATICA_3_0
159: argv[4] = "-linkmode";
160: switch(vmath->linkmode) {
161: case MATHEMATICA_LINK_CREATE:
162: argv[5] = "Create";
163: break;
164: case MATHEMATICA_LINK_CONNECT:
165: argv[5] = "Connect";
166: break;
167: case MATHEMATICA_LINK_LAUNCH:
168: argv[5] = "Launch";
169: break;
170: }
171: #else
172: switch(vmath->linkmode) {
173: case MATHEMATICA_LINK_CREATE:
174: argv[4] = "-linkcreate";
175: break;
176: case MATHEMATICA_LINK_CONNECT:
177: argv[4] = "-linkconnect";
178: break;
179: case MATHEMATICA_LINK_LAUNCH:
180: argv[4] = "-linklaunch";
181: break;
182: }
183: #endif
185: vmath->link = MLOpenInEnv(mathematicaEnv, argc, argv, &lerr);
186: #endif
188: return(0);
189: }
191: EXTERN_C_BEGIN
192: #undef __FUNCT__
194: int PetscViewerCreate_Mathematica(PetscViewer v) {
195: PetscViewer_Mathematica *vmath;
196: int ierr;
200: PetscNew(PetscViewer_Mathematica, &vmath);
201: v->data = (void *) vmath;
202: v->ops->destroy = PetscViewerDestroy_Mathematica;
203: v->ops->flush = 0;
204: PetscStrallocpy(PETSC_VIEWER_MATHEMATICA, &v->type_name);
206: vmath->linkname = PETSC_NULL;
207: vmath->linkhost = PETSC_NULL;
208: vmath->linkmode = MATHEMATICA_LINK_CONNECT;
209: vmath->graphicsType = GRAPHICS_MOTIF;
210: vmath->plotType = MATHEMATICA_TRIANGULATION_PLOT;
211: vmath->objName = PETSC_NULL;
213: PetscViewerMathematicaSetFromOptions(v);
214: PetscViewerMathematicaSetupConnection_Private(v);
215: return(0);
216: }
217: EXTERN_C_END
219: #undef __FUNCT__
221: int PetscViewerMathematicaParseLinkMode_Private(char *modename, LinkMode *mode) {
222: PetscTruth isCreate, isConnect, isLaunch;
223: int ierr;
226: PetscStrcasecmp(modename, "Create", &isCreate);
227: PetscStrcasecmp(modename, "Connect", &isConnect);
228: PetscStrcasecmp(modename, "Launch", &isLaunch);
229: if (isCreate == PETSC_TRUE) {
230: *mode = MATHEMATICA_LINK_CREATE;
231: } else if (isConnect == PETSC_TRUE) {
232: *mode = MATHEMATICA_LINK_CONNECT;
233: } else if (isLaunch == PETSC_TRUE) {
234: *mode = MATHEMATICA_LINK_LAUNCH;
235: } else {
236: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid Mathematica link mode: %s", modename);
237: }
238: return(0);
239: }
241: #undef __FUNCT__
243: int PetscViewerMathematicaSetFromOptions(PetscViewer v) {
244: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) v->data;
245: char linkname[256];
246: char modename[256];
247: char hostname[256];
248: char type[256];
249: int numPorts;
250: int *ports;
251: int numHosts;
252: char **hosts;
253: int size, rank;
254: int h;
255: PetscTruth opt;
256: int ierr;
259: MPI_Comm_size(v->comm, &size);
260: MPI_Comm_rank(v->comm, &rank);
262: /* Get link name */
263: PetscOptionsGetString("viewer_", "-math_linkname", linkname, 255, &opt);
264: if (opt == PETSC_TRUE) {
265: PetscViewerMathematicaSetLinkName(v, linkname);
266: }
267: /* Get link port */
268: numPorts = size;
269: PetscMalloc(size * sizeof(int), &ports);
270: PetscOptionsGetIntArray("viewer_", "-math_linkport", ports, &numPorts, &opt);
271: if (opt == PETSC_TRUE) {
272: if (numPorts > rank) {
273: snprintf(linkname, 255, "%6d", ports[rank]);
274: } else {
275: snprintf(linkname, 255, "%6d", ports[0]);
276: }
277: PetscViewerMathematicaSetLinkName(v, linkname);
278: }
279: PetscFree(ports);
280: /* Get link host */
281: numHosts = size;
282: PetscMalloc(size * sizeof(char *), &hosts);
283: PetscOptionsGetStringArray("viewer_", "-math_linkhost", hosts, &numHosts, &opt);
284: if (opt == PETSC_TRUE) {
285: if (numHosts > rank) {
286: PetscStrncpy(hostname, hosts[rank], 255);
287: } else {
288: PetscStrncpy(hostname, hosts[0], 255);
289: }
290: PetscViewerMathematicaSetLinkHost(v, hostname);
291: }
292: for(h = 0; h < numHosts; h++) {
293: PetscFree(hosts[h]);
294: }
295: PetscFree(hosts);
296: /* Get link mode */
297: PetscOptionsGetString("viewer_", "-math_linkmode", modename, 255, &opt);
298: if (opt == PETSC_TRUE) {
299: LinkMode mode;
301: PetscViewerMathematicaParseLinkMode_Private(modename, &mode);
302: PetscViewerMathematicaSetLinkMode(v, mode);
303: }
304: /* Get graphics type */
305: PetscOptionsGetString("viewer_", "-math_graphics", type, 255, &opt);
306: if (opt == PETSC_TRUE) {
307: PetscTruth isMotif, isPS, isPSFile;
309: PetscStrcasecmp(type, "Motif", &isMotif);
310: PetscStrcasecmp(type, "PS", &isPS);
311: PetscStrcasecmp(type, "PSFile", &isPSFile);
312: if (isMotif == PETSC_TRUE) {
313: vmath->graphicsType = GRAPHICS_MOTIF;
314: } else if (isPS == PETSC_TRUE) {
315: vmath->graphicsType = GRAPHICS_PS_STDOUT;
316: } else if (isPSFile == PETSC_TRUE) {
317: vmath->graphicsType = GRAPHICS_PS_FILE;
318: }
319: }
320: /* Get plot type */
321: PetscOptionsGetString("viewer_", "-math_type", type, 255, &opt);
322: if (opt == PETSC_TRUE) {
323: PetscTruth isTri, isVecTri, isVec, isSurface;
325: PetscStrcasecmp(type, "Triangulation", &isTri);
326: PetscStrcasecmp(type, "VectorTriangulation", &isVecTri);
327: PetscStrcasecmp(type, "Vector", &isVec);
328: PetscStrcasecmp(type, "Surface", &isSurface);
329: if (isTri == PETSC_TRUE) {
330: vmath->plotType = MATHEMATICA_TRIANGULATION_PLOT;
331: } else if (isVecTri == PETSC_TRUE) {
332: vmath->plotType = MATHEMATICA_VECTOR_TRIANGULATION_PLOT;
333: } else if (isVec == PETSC_TRUE) {
334: vmath->plotType = MATHEMATICA_VECTOR_PLOT;
335: } else if (isSurface == PETSC_TRUE) {
336: vmath->plotType = MATHEMATICA_SURFACE_PLOT;
337: }
338: }
339: return(0);
340: }
342: #undef __FUNCT__
344: int PetscViewerMathematicaSetLinkName(PetscViewer v, const char *name) {
345: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) v->data;
346: int ierr;
350: PetscStrallocpy(name, &vmath->linkname);
351: return(0);
352: }
354: #undef __FUNCT__
356: int PetscViewerMathematicaSetLinkPort(PetscViewer v, int port) {
357: char name[16];
358: int ierr;
361: snprintf(name, 16, "%6d", port);
362: PetscViewerMathematicaSetLinkName(v, name);
363: return(0);
364: }
366: #undef __FUNCT__
368: int PetscViewerMathematicaSetLinkHost(PetscViewer v, const char *host) {
369: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) v->data;
370: int ierr;
374: PetscStrallocpy(host, &vmath->linkhost);
375: return(0);
376: }
378: #undef __FUNCT__
380: int PetscViewerMathematicaSetLinkMode(PetscViewer v, LinkMode mode) {
381: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) v->data;
384: vmath->linkmode = mode;
385: return(0);
386: }
388: /*----------------------------------------- Public Functions --------------------------------------------------------*/
389: #undef __FUNCT__
391: /*@C
392: PetscViewerMathematicaOpen - Communicates with Mathemtica using MathLink.
394: Collective on comm
396: Input Parameters:
397: + comm - The MPI communicator
398: . port - [optional] The port to connect on, or PETSC_DECIDE
399: . machine - [optional] The machine to run Mathematica on, or PETSC_NULL
400: - mode - [optional] The connection mode, or PETSC_NULL
402: Output Parameter:
403: . viewer - The Mathematica viewer
405: Level: intermediate
407: Notes:
408: Most users should employ the following commands to access the
409: Mathematica viewers
410: $
411: $ PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
412: $ MatView(Mat matrix, PetscViewer viewer)
413: $
414: $ or
415: $
416: $ PetscViewerMathematicaOpen(MPI_Comm comm, int port, char *machine, char *mode, PetscViewer &viewer)
417: $ VecView(Vec vector, PetscViewer viewer)
419: Options Database Keys:
420: $ -viewer_math_linkhost <machine> - The host machine for the kernel
421: $ -viewer_math_linkname <name> - The full link name for the connection
422: $ -viewer_math_linkport <port> - The port for the connection
423: $ -viewer_math_mode <mode> - The mode, e.g. Launch, Connect
424: $ -viewer_math_type <type> - The plot type, e.g. Triangulation, Vector
425: $ -viewer_math_graphics <output> - The output type, e.g. Motif, PS, PSFile
427: .keywords: PetscViewer, Mathematica, open
429: .seealso: MatView(), VecView()
430: @*/
431: int PetscViewerMathematicaOpen(MPI_Comm comm, int port, const char machine[], const char mode[], PetscViewer *v) {
432: int ierr;
435: PetscViewerCreate(comm, v);
436: #if 0
437: LinkMode linkmode;
438: PetscViewerMathematicaSetLinkPort(*v, port);
439: PetscViewerMathematicaSetLinkHost(*v, machine);
440: PetscViewerMathematicaParseLinkMode_Private(mode, &linkmode);
441: PetscViewerMathematicaSetLinkMode(*v, linkmode);
442: #endif
443: PetscViewerSetType(*v, PETSC_VIEWER_MATHEMATICA);
444: return(0);
445: }
447: #ifdef PETSC_HAVE_MATHEMATICA
448: #undef __FUNCT__
450: /*@C
451: PetscViewerMathematicaGetLink - Returns the link to Mathematica
453: Input Parameters:
454: . viewer - The Mathematica viewer
455: . link - The link to Mathematica
457: Level: intermediate
459: .keywords PetscViewer, Mathematica, link
460: .seealso PetscViewerMathematicaOpen()
461: @*/
462: int PetscViewerMathematicaGetLink(PetscViewer viewer, MLINK *link)
463: {
464: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
468: *link = vmath->link;
469: return(0);
470: }
471: #endif
473: #undef __FUNCT__
475: /*@C
476: PetscViewerMathematicaSkipPackets - Discard packets sent by Mathematica until a certain packet type is received
478: Input Parameters:
479: . viewer - The Mathematica viewer
480: . type - The packet type to search for, e.g RETURNPKT
482: Level: advanced
484: .keywords PetscViewer, Mathematica, packets
485: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaGetVector()
486: @*/
487: int PetscViewerMathematicaSkipPackets(PetscViewer viewer, int type)
488: {
489: #ifdef PETSC_HAVE_MATHEMATICA
490: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
491: MLINK link = vmath->link; /* The link to Mathematica */
492: int pkt; /* The packet type */
495: while((pkt = MLNextPacket(link)) && (pkt != type))
496: MLNewPacket(link);
497: if (!pkt) {
498: MLClearError(link);
499: SETERRQ(PETSC_ERR_LIB, (char *) MLErrorMessage(link));
500: }
501: return(0);
502: #endif
504: return(0);
505: }
507: #undef __FUNCT__
509: /*@C
510: PetscViewerMathematicaLoadGraphics - Change the type of graphics output for Mathematica
512: Input Parameters:
513: . viewer - The Mathematica viewer
514: . type - The type of graphics, e.g. GRAPHICS_MOTIF, GRAPHICS_PS_FILE, GRAPHICS_PS_STDOUT
516: Level: intermediate
518: .keywords PetscViewer, Mathematica, packets
519: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaSkipPackets()
520: @*/
521: int PetscViewerMathematicaLoadGraphics(PetscViewer viewer, GraphicsType type)
522: {
523: #ifdef PETSC_HAVE_MATHEMATICA
524: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
525: MLINK link = vmath->link; /* The link to Mathematica */
526: char programName[256];
527: int ierr;
530: /* Load graphics package */
531: MLPutFunction(link, "EvaluatePacket", 1);
532: switch(type)
533: {
534: case GRAPHICS_MOTIF:
535: MLPutFunction(link, "Get", 1);
536: MLPutString(link, "Motif.m");
537: break;
538: case GRAPHICS_PS_FILE:
539: MLPutFunction(link, "CompoundExpression", 4);
540: MLPutFunction(link, "Set", 2);
541: MLPutSymbol(link, "PetscGraphicsCounter");
542: MLPutInteger(link, 0);
543: MLPutFunction(link, "SetDelayed", 2);
544: MLPutSymbol(link, "$Display");
545: MLPutSymbol(link, "$FileDisplay");
546: MLPutFunction(link, "Set", 2);
547: MLPutSymbol(link, "$FileDisplay");
548: MLPutFunction(link, "OpenWrite", 1);
549: MLPutFunction(link, "StringJoin", 3);
550: if (!PetscGetProgramName(programName, 255))
551: MLPutString(link, programName);
552: else
553: MLPutString(link, "GVec");
554: MLPutFunction(link, "ToString", 1);
555: MLPutSymbol(link, "PetscGraphicsCounter");
556: MLPutString(link, ".ps");
557: MLPutFunction(link, "Set", 2);
558: MLPutSymbol(link, "$DisplayFunction");
559: MLPutFunction(link, "Function", 1);
560: MLPutFunction(link, "CompoundExpression", 2);
561: MLPutFunction(link, "Display", 3);
562: MLPutSymbol(link, "$Display");
563: MLPutFunction(link, "Slot", 1);
564: MLPutInteger(link, 1);
565: MLPutString(link, "EPS");
566: MLPutFunction(link, "Increment", 1);
567: MLPutSymbol(link, "PetscGraphicsCounter");
568: break;
569: case GRAPHICS_PS_STDOUT:
570: MLPutFunction(link, "Get", 1);
571: MLPutString(link, "PSDirect.m");
572: break;
573: default:
574: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid graphics type: %d", type);
575: }
576: MLEndPacket(link);
578: /* Skip packets until ReturnPacket */
579: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
580: /* Skip ReturnPacket */
581: MLNewPacket(link);
583: /* Load PlotField.m for vector plots */
584: MLPutFunction(link, "EvaluatePacket", 1);
585: MLPutFunction(link, "Get", 1);
586: MLPutString(link, "Graphics/PlotField.m");
587: MLEndPacket(link);
589: /* Skip packets until ReturnPacket */
590: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
591: /* Skip ReturnPacket */
592: MLNewPacket(link);
594: return(0);
595: #else
597: switch(type)
598: {
599: case GRAPHICS_MOTIF:
600: case GRAPHICS_PS_FILE:
601: case GRAPHICS_PS_STDOUT:
602: break;
603: default:
604: SETERRQ1(PETSC_ERR_ARG_WRONG, "Invalid graphics type: %d", type);
605: }
606: return(0);
607: #endif
608: }
610: #undef __FUNCT__
612: /*@C
613: PetscViewerMathematicaGetName - Retrieve the default name for objects communicated to Mathematica
615: Input Parameter:
616: . viewer - The Mathematica viewer
618: Output Parameter:
619: . name - The name for new objects created in Mathematica
621: Level: intermediate
623: .keywords PetscViewer, Mathematica, name
624: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaClearName()
625: @*/
626: int PetscViewerMathematicaGetName(PetscViewer viewer, const char **name)
627: {
628: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
633: *name = vmath->objName;
634: return(0);
635: }
637: #undef __FUNCT__
639: /*@C
640: PetscViewerMathematicaSetName - Override the default name for objects communicated to Mathematica
642: Input Parameters:
643: . viewer - The Mathematica viewer
644: . name - The name for new objects created in Mathematica
646: Level: intermediate
648: .keywords PetscViewer, Mathematica, name
649: .seealso PetscViewerMathematicaSetName(), PetscViewerMathematicaClearName()
650: @*/
651: int PetscViewerMathematicaSetName(PetscViewer viewer, const char name[])
652: {
653: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
658: vmath->objName = name;
659: return(0);
660: }
662: #undef __FUNCT__
664: /*@C
665: PetscViewerMathematicaClearName - Use the default name for objects communicated to Mathematica
667: Input Parameter:
668: . viewer - The Mathematica viewer
670: Level: intermediate
672: .keywords PetscViewer, Mathematica, name
673: .seealso PetscViewerMathematicaGetName(), PetscViewerMathematicaSetName()
674: @*/
675: int PetscViewerMathematicaClearName(PetscViewer viewer)
676: {
677: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
681: vmath->objName = PETSC_NULL;
682: return(0);
683: }
685: #undef __FUNCT__
687: /*@C
688: PetscViewerMathematicaGetVector - Retrieve a vector from Mathematica
690: Input Parameter:
691: . viewer - The Mathematica viewer
693: Output Parameter:
694: . v - The vector
696: Level: intermediate
698: .keywords PetscViewer, Mathematica, vector
699: .seealso VecView(), PetscViewerMathematicaPutVector()
700: @*/
701: int PetscViewerMathematicaGetVector(PetscViewer viewer, Vec v) {
702: #ifdef PETSC_HAVE_MATHEMATICA
703: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
704: MLINK link; /* The link to Mathematica */
705: char *name;
706: PetscScalar *mArray;
707: PetscScalar *array;
708: long mSize;
709: int size;
710: int ierr;
716: /* Determine the object name */
717: if (vmath->objName == PETSC_NULL) {
718: name = "vec";
719: } else {
720: name = (char *) vmath->objName;
721: }
723: link = vmath->link;
724: VecGetLocalSize(v, &size);
725: VecGetArray(v, &array);
726: MLPutFunction(link, "EvaluatePacket", 1);
727: MLPutSymbol(link, name);
728: MLEndPacket(link);
729: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
730: MLGetRealList(link, &mArray, &mSize);
731: if (size != mSize) SETERRQ(PETSC_ERR_ARG_WRONG, "Incompatible vector sizes");
732: PetscMemcpy(array, mArray, mSize * sizeof(double));
733: MLDisownRealList(link, mArray, mSize);
734: VecRestoreArray(v, &array);
736: return(0);
737: #endif
739: return(0);
740: }
742: #undef __FUNCT__
744: /*@C
745: PetscViewerMathematicaPutVector - Send a vector to Mathematica
747: Input Parameters:
748: + viewer - The Mathematica viewer
749: - v - The vector
751: Level: intermediate
753: .keywords PetscViewer, Mathematica, vector
754: .seealso VecView(), PetscViewerMathematicaGetVector()
755: @*/
756: int PetscViewerMathematicaPutVector(PetscViewer viewer, Vec v) {
757: #ifdef PETSC_HAVE_MATHEMATICA
758: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
759: MLINK link = vmath->link; /* The link to Mathematica */
760: char *name;
761: PetscScalar *array;
762: int size;
763: int ierr;
766: /* Determine the object name */
767: if (vmath->objName == PETSC_NULL) {
768: name = "vec";
769: } else {
770: name = (char *) vmath->objName;
771: }
772: VecGetLocalSize(v, &size);
773: VecGetArray(v, &array);
775: /* Send the Vector object */
776: MLPutFunction(link, "EvaluatePacket", 1);
777: MLPutFunction(link, "Set", 2);
778: MLPutSymbol(link, name);
779: MLPutRealList(link, array, size);
780: MLEndPacket(link);
781: /* Skip packets until ReturnPacket */
782: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
783: /* Skip ReturnPacket */
784: MLNewPacket(link);
786: VecRestoreArray(v, &array);
787: return(0);
788: #endif
790: return(0);
791: }
793: int PetscViewerMathematicaPutMatrix(PetscViewer viewer, int m, int n, PetscReal *a) {
794: #ifdef PETSC_HAVE_MATHEMATICA
795: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
796: MLINK link = vmath->link; /* The link to Mathematica */
797: char *name;
798: int ierr;
801: /* Determine the object name */
802: if (vmath->objName == PETSC_NULL) {
803: name = "mat";
804: } else {
805: name = (char *) vmath->objName;
806: }
808: /* Send the dense matrix object */
809: MLPutFunction(link, "EvaluatePacket", 1);
810: MLPutFunction(link, "Set", 2);
811: MLPutSymbol(link, name);
812: MLPutFunction(link, "Transpose", 1);
813: MLPutFunction(link, "Partition", 2);
814: MLPutRealList(link, a, m*n);
815: MLPutInteger(link, m);
816: MLEndPacket(link);
817: /* Skip packets until ReturnPacket */
818: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
819: /* Skip ReturnPacket */
820: MLNewPacket(link);
822: return(0);
823: #endif
825: return(0);
826: }
828: int PetscViewerMathematicaPutCSRMatrix(PetscViewer viewer, int m, int n, int *i, int *j, PetscReal *a) {
829: #ifdef PETSC_HAVE_MATHEMATICA
830: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
831: MLINK link = vmath->link; /* The link to Mathematica */
832: const char *symbol;
833: char *name;
834: PetscTruth match;
835: int ierr;
838: /* Determine the object name */
839: if (vmath->objName == PETSC_NULL) {
840: name = "mat";
841: } else {
842: name = (char *) vmath->objName;
843: }
845: /* Make sure Mathematica recognizes sparse matrices */
846: MLPutFunction(link, "EvaluatePacket", 1);
847: MLPutFunction(link, "Needs", 1);
848: MLPutString(link, "LinearAlgebra`CSRMatrix`");
849: MLEndPacket(link);
850: /* Skip packets until ReturnPacket */
851: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
852: /* Skip ReturnPacket */
853: MLNewPacket(link);
855: /* Send the CSRMatrix object */
856: MLPutFunction(link, "EvaluatePacket", 1);
857: MLPutFunction(link, "Set", 2);
858: MLPutSymbol(link, name);
859: MLPutFunction(link, "CSRMatrix", 5);
860: MLPutInteger(link, m);
861: MLPutInteger(link, n);
862: MLPutFunction(link, "Plus", 2);
863: MLPutIntegerList(link, i, m+1);
864: MLPutInteger(link, 1);
865: MLPutFunction(link, "Plus", 2);
866: MLPutIntegerList(link, j, i[m]);
867: MLPutInteger(link, 1);
868: MLPutRealList(link, a, i[m]);
869: MLEndPacket(link);
870: /* Skip packets until ReturnPacket */
871: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
872: /* Skip ReturnPacket */
873: MLNewPacket(link);
875: /* Check that matrix is valid */
876: MLPutFunction(link, "EvaluatePacket", 1);
877: MLPutFunction(link, "ValidQ", 1);
878: MLPutSymbol(link, name);
879: MLEndPacket(link);
880: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
881: MLGetSymbol(link, &symbol);
882: PetscStrcmp("True", (char *) symbol, &match);
883: if (match == PETSC_FALSE) {
884: MLDisownSymbol(link, symbol);
885: SETERRQ(PETSC_ERR_PLIB, "Invalid CSR matrix in Mathematica");
886: }
887: MLDisownSymbol(link, symbol);
888: /* Skip ReturnPacket */
889: MLNewPacket(link);
891: return(0);
892: #endif
894: return(0);
895: }
897: /*------------------------------------------- ML Functions ----------------------------------------------------------*/
898: #undef __FUNCT__
900: int PetscViewerMathematicaPartitionMesh(PetscViewer viewer, int numElements, int numVertices, double *vertices, int **mesh,
901: int *numParts, int *colPartition)
902: {
903: #ifdef PETSC_HAVE_MATHEMATICA
904: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
905: MLINK link; /* The link to Mathematica */
906: const char *symbol;
907: int numOptions;
908: int partSize;
909: int *part;
910: long numCols;
911: int col;
912: PetscTruth match, opt;
913: int ierr;
917: link = vmath->link;
919: /* Make sure that the reduce.m package is loaded */
920: MLPutFunction(link, "EvaluatePacket", 1);
921: MLPutFunction(link, "Get", 1);
922: MLPutFunction(link, "StringJoin", 2);
923: MLPutFunction(link, "Environment", 1);
924: MLPutString(link, "PETSC_DIR");
925: MLPutString(link, "/src/pc/impls/ml/reduce.m");
926: MLEndPacket(link);
927: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
928: MLGetSymbol(link, &symbol);
929: PetscStrcmp("$Failed", (char *) symbol, &match);
930: if (match == PETSC_TRUE) {
931: MLDisownSymbol(link, symbol);
932: SETERRQ(PETSC_ERR_FILE_OPEN, "Unable to load package reduce.m");
933: }
934: MLDisownSymbol(link, symbol);
936: /* Partition the mesh */
937: numOptions = 0;
938: partSize = 0;
939: PetscOptionsGetInt(PETSC_NULL, "-pc_ml_partition_size", &partSize, &opt);
940: if ((opt == PETSC_TRUE) && (partSize > 0))
941: numOptions++;
942: MLPutFunction(link, "EvaluatePacket", 1);
943: MLPutFunction(link, "PartitionMesh", 1 + numOptions);
944: MLPutFunction(link, "MeshData", 5);
945: MLPutInteger(link, numElements);
946: MLPutInteger(link, numVertices);
947: MLPutInteger(link, numVertices);
948: MLPutFunction(link, "Partition", 2);
949: MLPutRealList(link, vertices, numVertices*2);
950: MLPutInteger(link, 2);
951: MLPutFunction(link, "Partition", 2);
952: MLPutIntegerList(link, mesh[MESH_ELEM], numElements*3);
953: MLPutInteger(link, 3);
954: if (partSize > 0)
955: {
956: MLPutFunction(link, "Rule", 2);
957: MLPutSymbol(link, "PartitionSize");
958: MLPutInteger(link, partSize);
959: }
960: MLEndPacket(link);
961: /* Skip packets until ReturnPacket */
962: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
964: /* Get the vertex partiton */
965: MLGetIntegerList(link, &part, &numCols);
966: if (numCols != numVertices) SETERRQ(PETSC_ERR_PLIB, "Invalid partition");
967: for(col = 0, *numParts = 0; col < numCols; col++) {
968: colPartition[col] = part[col]-1;
969: *numParts = PetscMax(*numParts, part[col]);
970: }
971: MLDisownIntegerList(link, part, numCols);
973: return(0);
974: #endif
976: return(0);
977: }
979: #undef __FUNCT__
981: int PetscViewerMathematicaReduce(PetscViewer viewer, PC pc, int thresh)
982: {
983: #ifdef PETSC_HAVE_MATHEMATICA
984: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
985: MLINK link; /* The link to Mathematica */
986: PC_Multilevel *ml;
987: int *range;
988: long numRange;
989: int *null;
990: long numNull;
991: const char *symbol;
992: int numOptions;
993: int partSize;
994: int row, col;
995: PetscTruth match, opt;
996: int ierr;
1001: link = vmath->link;
1002: ml = (PC_Multilevel *) pc->data;
1004: /* Make sure that the reduce.m package is loaded */
1005: MLPutFunction(link, "EvaluatePacket", 1);
1006: MLPutFunction(link, "Get", 1);
1007: MLPutFunction(link, "StringJoin", 2);
1008: MLPutFunction(link, "Environment", 1);
1009: MLPutString(link, "PETSC_DIR");
1010: MLPutString(link, "/src/pc/impls/ml/reduce.m");
1011: MLEndPacket(link);
1012: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1013: MLGetSymbol(link, &symbol);
1014: PetscStrcmp("$Failed", (char *) symbol, &match);
1015: if (match == PETSC_TRUE) {
1016: MLDisownSymbol(link, symbol);
1017: SETERRQ(PETSC_ERR_FILE_OPEN, "Unable to load package reduce.m");
1018: }
1019: MLDisownSymbol(link, symbol);
1021: /* mesh = MeshData[] */
1022: MLPutFunction(link, "EvaluatePacket", 1);
1023: MLPutFunction(link, "Set", 2);
1024: MLPutSymbol(link, "mesh");
1025: MLPutFunction(link, "MeshData", 5);
1026: MLPutInteger(link, ml->numElements[0]);
1027: MLPutInteger(link, ml->numVertices[0]);
1028: MLPutInteger(link, ml->numVertices[0]);
1029: MLPutFunction(link, "Partition", 2);
1030: MLPutRealList(link, ml->vertices, ml->numVertices[0]*2);
1031: MLPutInteger(link, 2);
1032: MLPutFunction(link, "Partition", 2);
1033: MLPutIntegerList(link, ml->meshes[0][MESH_ELEM], ml->numElements[0]*3);
1034: MLPutInteger(link, 3);
1035: MLEndPacket(link);
1036: /* Skip packets until ReturnPacket */
1037: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1038: /* Skip ReturnPacket */
1039: MLNewPacket(link);
1040: /* Check that mesh is valid */
1041: MLPutFunction(link, "EvaluatePacket", 1);
1042: MLPutFunction(link, "ValidQ", 1);
1043: MLPutSymbol(link, "mesh");
1044: MLEndPacket(link);
1045: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1046: MLGetSymbol(link, &symbol);
1047: PetscStrcmp("True", (char *) symbol, &match);
1048: if (match == PETSC_FALSE) {
1049: MLDisownSymbol(link, symbol);
1050: SETERRQ(PETSC_ERR_PLIB, "Invalid mesh in Mathematica");
1051: }
1052: MLDisownSymbol(link, symbol);
1054: /* mat = gradient matrix */
1055: MatView(ml->locB, viewer);
1057: /* mattML = ReduceMatrix[mat,ml->minNodes] */
1058: numOptions = 0;
1059: partSize = 0;
1060: PetscOptionsGetInt(PETSC_NULL, "-pc_ml_partition_size", &partSize, &opt);
1061: if ((opt == PETSC_TRUE) && (partSize > 0))
1062: numOptions++;
1063: MLPutFunction(link, "EvaluatePacket", 1);
1064: MLPutFunction(link, "Set", 2);
1065: MLPutSymbol(link, "mattML");
1066: MLPutFunction(link, "ReduceMatrix", 3 + numOptions);
1067: MLPutSymbol(link, "mesh");
1068: MLPutSymbol(link, "mat");
1069: MLPutInteger(link, thresh);
1070: if (partSize > 0) {
1071: MLPutFunction(link, "Rule", 2);
1072: MLPutSymbol(link, "PartitionSize");
1073: MLPutInteger(link, partSize);
1074: }
1075: MLEndPacket(link);
1076: /* Skip packets until ReturnPacket */
1077: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1078: /* Skip ReturnPacket */
1079: MLNewPacket(link);
1080: /* Check that mattML is valid */
1081: MLPutFunction(link, "EvaluatePacket", 1);
1082: MLPutFunction(link, "ValidQ", 1);
1083: MLPutSymbol(link, "mattML");
1084: MLEndPacket(link);
1085: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1086: MLGetSymbol(link, &symbol);
1087: PetscStrcmp("True", (char *) symbol, &match);
1088: if (match == PETSC_FALSE) {
1089: MLDisownSymbol(link, symbol);
1090: SETERRQ(PETSC_ERR_PLIB, "Invalid MLData in Mathematica");
1091: }
1092: MLDisownSymbol(link, symbol);
1094: /* Copy information to the preconditioner */
1095: MLPutFunction(link, "EvaluatePacket", 1);
1096: MLPutFunction(link, "Part", 2);
1097: MLPutSymbol(link, "mattML");
1098: MLPutInteger(link, 3);
1099: MLEndPacket(link);
1100: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1101: MLGetInteger(link, &ml->numLevels);
1103: /* Create lists of the range and nullspace columns */
1104: MLPutFunction(link, "EvaluatePacket", 1);
1105: MLPutFunction(link, "GetRange", 1);
1106: MLPutSymbol(link, "mattML");
1107: MLEndPacket(link);
1108: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1109: MLGetIntegerList(link, &range, &numRange);
1110: if (numRange > ml->sOrder->numLocVars) SETERRQ(PETSC_ERR_PLIB, "Invalid size for range space");
1111: ml->rank = numRange;
1112: ml->globalRank = ml->rank;
1113: if (ml->rank > 0) {
1114: PetscMalloc(numRange * sizeof(int), &ml->range);
1115: for(row = 0; row < numRange; row++)
1116: ml->range[row] = range[row]-1;
1117: }
1118: MLDisownIntegerList(link, range, numRange);
1120: MLPutFunction(link, "EvaluatePacket", 1);
1121: MLPutFunction(link, "GetNullColumns", 1);
1122: MLPutSymbol(link, "mattML");
1123: MLEndPacket(link);
1124: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1125: MLGetIntegerList(link, &null, &numNull);
1126: if (numRange + numNull != ml->sOrder->numLocVars) SETERRQ(PETSC_ERR_PLIB, "Invalid size for range and null spaces");
1127: ml->numLocNullCols = numNull;
1128: if (numNull > 0)
1129: {
1130: PetscMalloc(numNull * sizeof(int), &ml->nullCols);
1131: for(col = 0; col < numNull; col++)
1132: ml->nullCols[col] = null[col] - 1;
1133: }
1134: MLDisownIntegerList(link, null, numNull);
1136: return(0);
1137: #endif
1139: return(0);
1140: }
1142: #undef __FUNCT__
1144: int PetscViewerMathematicaMultiLevelConvert(PetscViewer viewer, PC pc)
1145: {
1146: #ifdef PETSC_HAVE_MATHEMATICA
1147: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
1148: MLINK link; /* The link to Mathematica */
1149: PC_Multilevel *ml;
1150: Mat_SeqAIJ *grad;
1151: int *numPartitions;
1152: int *numPartitionCols, *cols;
1153: int *numPartitionRows, *rows;
1154: double *U, *singVals, *V;
1155: long *Udims, *Vdims;
1156: char **Uheads, **Vheads;
1157: int *nnz;
1158: int *offsets;
1159: double *vals;
1160: long numLevels, numParts, numCols, numRows, Udepth, numSingVals, Vdepth, len;
1161: int numBdRows, numBdCols;
1162: int mesh, level, part, col, row;
1163: int ierr;
1168: link = vmath->link;
1169: ml = (PC_Multilevel *) pc->data;
1171: /* ml->numLevels = ml[[3]] */
1172: MLPutFunction(link, "EvaluatePacket", 1);
1173: MLPutFunction(link, "Part", 2);
1174: MLPutSymbol(link, "mattML");
1175: MLPutInteger(link, 3);
1176: MLEndPacket(link);
1177: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1178: MLGetInteger(link, &ml->numLevels);
1180: /* ml->numMeshes = Length[ml[[4]]] */
1181: MLPutFunction(link, "EvaluatePacket", 1);
1182: MLPutFunction(link, "Length", 1);
1183: MLPutFunction(link, "Part", 2);
1184: MLPutSymbol(link, "mattML");
1185: MLPutInteger(link, 4);
1186: MLEndPacket(link);
1187: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1188: MLGetInteger(link, &ml->numMeshes);
1190: /* ml->numElements[level] = ml[[4,level,1]] */
1191: PetscMalloc(ml->numMeshes * sizeof(int), &ml->numElements);
1193: /* ml->numVertices[level] = ml[[4,level,2]] */
1194: PetscMalloc(ml->numMeshes * sizeof(int), &ml->numVertices);
1196: /* ml->meshes = ml[[4]] */
1197: PetscMalloc(ml->numMeshes * sizeof(int **), &ml->meshes);
1198: for(mesh = 0; mesh < ml->numMeshes; mesh++) {
1199: PetscMalloc(NUM_MESH_DIV * sizeof(int *), &ml->meshes[mesh]);
1200: /* Here we should get meshes */
1201: PetscMalloc(1 * sizeof(int), &ml->meshes[mesh][MESH_OFFSETS]);
1202: PetscMalloc(1 * sizeof(int), &ml->meshes[mesh][MESH_ADJ]);
1203: PetscMalloc(1 * sizeof(int), &ml->meshes[mesh][MESH_ELEM]);
1204: }
1206: /* ml->numPartitions = Map[Length,Drop[ml[[5]],-1]] */
1207: MLPutFunction(link, "EvaluatePacket", 1);
1208: MLPutFunction(link, "Map", 2);
1209: MLPutSymbol(link, "Length");
1210: MLPutFunction(link, "Drop", 2);
1211: MLPutFunction(link, "Part", 2);
1212: MLPutSymbol(link, "mattML");
1213: MLPutInteger(link, 5);
1214: MLPutInteger(link, -1);
1215: MLEndPacket(link);
1216: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1217: MLGetIntegerList(link, &numPartitions, &numLevels);
1218: if (numLevels != ml->numLevels) SETERRQ(PETSC_ERR_PLIB, "Invalid node partition in MLData object");
1219: if (numLevels > 0) {
1220: PetscMalloc(numLevels * sizeof(int), &ml->numPartitions);
1221: PetscMemcpy(ml->numPartitions, numPartitions, numLevels * sizeof(int));
1222: }
1223: MLDisownIntegerList(link, numPartitions, numLevels);
1225: if (ml->numLevels > 0) {
1226: /* ml->numPartitionCols = Map[Length,ml[[5,level]]] */
1227: PetscMalloc(ml->numLevels * sizeof(int *), &ml->numPartitionCols);
1228: for(level = 0; level < ml->numLevels; level++) {
1229: MLPutFunction(link, "EvaluatePacket", 1);
1230: MLPutFunction(link, "Map", 2);
1231: MLPutSymbol(link, "Length");
1232: MLPutFunction(link, "Part", 3);
1233: MLPutSymbol(link, "mattML");
1234: MLPutInteger(link, 5);
1235: MLPutInteger(link, level+1);
1236: MLEndPacket(link);
1237: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1238: MLGetIntegerList(link, &numPartitionCols, &numParts);
1239: if (numParts != ml->numPartitions[level]) SETERRQ(PETSC_ERR_PLIB, "Invalid node partition in MLData object");
1240: if (numParts > 0) {
1241: PetscMalloc(numParts * sizeof(int), &ml->numPartitionCols[level]);
1242: PetscMemcpy(ml->numPartitionCols[level], numPartitionCols, numParts * sizeof(int));
1243: }
1244: MLDisownIntegerList(link, numPartitionCols, numParts);
1245: }
1247: /* ml->colPartition[level][part] = ml[[5,level,part]] */
1248: PetscMalloc(ml->numLevels * sizeof(int **), &ml->colPartition);
1249: for(level = 0; level < ml->numLevels; level++) {
1250: if (ml->numPartitions[level] == 0) continue;
1251: PetscMalloc(ml->numPartitions[level] * sizeof(int *), &ml->colPartition[level]);
1252: for(part = 0; part < ml->numPartitions[level]; part++) {
1253: MLPutFunction(link, "EvaluatePacket", 1);
1254: MLPutFunction(link, "Part", 4);
1255: MLPutSymbol(link, "mattML");
1256: MLPutInteger(link, 5);
1257: MLPutInteger(link, level+1);
1258: MLPutInteger(link, part+1);
1259: MLEndPacket(link);
1260: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1261: MLGetIntegerList(link, &cols, &numCols);
1262: if (numCols != ml->numPartitionCols[level][part]) SETERRQ(PETSC_ERR_PLIB, "Invalid node partition in MLData object");
1263: if (numCols > 0) {
1264: PetscMalloc(numCols * sizeof(int), &ml->colPartition[level][part]);
1265: /* Convert to zero-based indices */
1266: for(col = 0; col < numCols; col++) ml->colPartition[level][part][col] = cols[col] - 1;
1267: }
1268: MLDisownIntegerList(link, cols, numCols);
1269: }
1270: }
1272: /* ml->numPartitionRows = Map[Length,FlattenAt[ml[[6,level]],1]] */
1273: PetscMalloc(ml->numLevels * sizeof(int *), &ml->numPartitionRows);
1274: for(level = 0; level < ml->numLevels; level++) {
1275: MLPutFunction(link, "EvaluatePacket", 1);
1276: MLPutFunction(link, "Map", 2);
1277: MLPutSymbol(link, "Length");
1278: MLPutFunction(link, "FlattenAt", 2);
1279: MLPutFunction(link, "Part", 3);
1280: MLPutSymbol(link, "mattML");
1281: MLPutInteger(link, 6);
1282: MLPutInteger(link, level+1);
1283: MLPutInteger(link, 1);
1284: MLEndPacket(link);
1285: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1286: MLGetIntegerList(link, &numPartitionRows, &numParts);
1287: if (numParts != ml->numPartitions[level] + NUM_PART_ROW_DIV - 1) {
1288: SETERRQ(PETSC_ERR_PLIB, "Invalid edge partition in MLData object");
1289: }
1290: PetscMalloc(numParts * sizeof(int), &ml->numPartitionRows[level]);
1291: PetscMemcpy(ml->numPartitionRows[level], numPartitionRows, numParts * sizeof(int));
1292: MLDisownIntegerList(link, numPartitionRows, numParts);
1293: }
1295: /* ml->rowPartition[level][PART_ROW_INT][part] = ml[[6,level,1,part]]
1296: ml->rowPartition[level][PART_ROW_BD] = ml[[6,level,2]]
1297: ml->rowPartition[level][PART_ROW_RES] = ml[[6,level,3]] */
1298: PetscMalloc(ml->numLevels * sizeof(int ***), &ml->rowPartition);
1299: for(level = 0; level < ml->numLevels; level++) {
1300: PetscMalloc(NUM_PART_ROW_DIV * sizeof(int **), &ml->rowPartition[level]);
1301: /* Interior rows */
1302: if (ml->numPartitions[level] > 0) {
1303: PetscMalloc(ml->numPartitions[level] * sizeof(int *), &ml->rowPartition[level][PART_ROW_INT]);
1304: for(part = 0; part < ml->numPartitions[level]; part++) {
1305: MLPutFunction(link, "EvaluatePacket", 1);
1306: MLPutFunction(link, "Part", 5);
1307: MLPutSymbol(link, "mattML");
1308: MLPutInteger(link, 6);
1309: MLPutInteger(link, level+1);
1310: MLPutInteger(link, 1);
1311: MLPutInteger(link, part+1);
1312: MLEndPacket(link);
1313: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1314: MLGetIntegerList(link, &rows, &numRows);
1315: if (numRows != ml->numPartitionRows[level][part]) {
1316: SETERRQ(PETSC_ERR_PLIB, "Invalid edge partition in MLData object");
1317: }
1318: if (numRows > 0) {
1319: PetscMalloc(numRows * sizeof(int), &ml->rowPartition[level][PART_ROW_INT][part]);
1320: /* Convert to zero-based indices */
1321: for(row = 0; row < numRows; row++) {
1322: ml->rowPartition[level][PART_ROW_INT][part][row] = rows[row] - 1;
1323: }
1324: }
1325: MLDisownIntegerList(link, rows, numRows);
1326: }
1327: }
1328: /* Boundary rows */
1329: PetscMalloc(1 * sizeof(int *), &ml->rowPartition[level][PART_ROW_BD]);
1330: MLPutFunction(link, "EvaluatePacket", 1);
1331: MLPutFunction(link, "Part", 4);
1332: MLPutSymbol(link, "mattML");
1333: MLPutInteger(link, 6);
1334: MLPutInteger(link, level+1);
1335: MLPutInteger(link, 2);
1336: MLEndPacket(link);
1337: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1338: MLGetIntegerList(link, &rows, &numRows);
1339: if (numRows != ml->numPartitionRows[level][ml->numPartitions[level]]) {
1340: SETERRQ(PETSC_ERR_PLIB, "Invalid edge partition in MLData object");
1341: }
1342: if (numRows > 0) {
1343: PetscMalloc(numRows * sizeof(int), &ml->rowPartition[level][PART_ROW_BD][0]);
1344: /* Convert to zero-based indices */
1345: for(row = 0; row < numRows; row++) {
1346: ml->rowPartition[level][PART_ROW_BD][0][row] = rows[row] - 1;
1347: }
1348: }
1349: MLDisownIntegerList(link, rows, numRows);
1350: /* Residual rows*/
1351: PetscMalloc(1 * sizeof(int *), &ml->rowPartition[level][PART_ROW_RES]);
1352: MLPutFunction(link, "EvaluatePacket", 1);
1353: MLPutFunction(link, "Part", 4);
1354: MLPutSymbol(link, "mattML");
1355: MLPutInteger(link, 6);
1356: MLPutInteger(link, level+1);
1357: MLPutInteger(link, 3);
1358: MLEndPacket(link);
1359: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1360: MLGetIntegerList(link, &rows, &numRows);
1361: if (numRows != ml->numPartitionRows[level][ml->numPartitions[level]+1]) {
1362: SETERRQ(PETSC_ERR_PLIB, "Invalid edge partition in MLData object");
1363: }
1364: if (numRows > 0) {
1365: PetscMalloc(numRows * sizeof(int), &ml->rowPartition[level][PART_ROW_RES][0]);
1366: /* Convert to zero-based indices */
1367: for(row = 0; row < numRows; row++) {
1368: ml->rowPartition[level][PART_ROW_RES][0][row] = rows[row] - 1;
1369: }
1370: }
1371: MLDisownIntegerList(link, rows, numRows);
1372: }
1373: } else {
1374: PetscMalloc(1 * sizeof(int), &ml->numPartitions);
1375: PetscMalloc(1 * sizeof(int *), &ml->numPartitionCols);
1376: PetscMalloc(1 * sizeof(int **), &ml->colPartition);
1377: PetscMalloc(1 * sizeof(int *), &ml->numPartitionRows);
1378: PetscMalloc(1 * sizeof(int ***), &ml->rowPartition);
1379: }
1381: /* ml->numRows = ml[[1]] */
1382: MLPutFunction(link, "EvaluatePacket", 1);
1383: MLPutFunction(link, "Part", 2);
1384: MLPutSymbol(link, "mattML");
1385: MLPutInteger(link, 1);
1386: MLEndPacket(link);
1387: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1388: MLGetInteger(link, &ml->numRows);
1390: /* ml->numCols = ml[[2]] */
1391: MLPutFunction(link, "EvaluatePacket", 1);
1392: MLPutFunction(link, "Part", 2);
1393: MLPutSymbol(link, "mattML");
1394: MLPutInteger(link, 2);
1395: MLEndPacket(link);
1396: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1397: MLGetInteger(link, &ml->numCols);
1399: /* ml->zeroTol = ml[[9]] */
1400: MLPutFunction(link, "EvaluatePacket", 1);
1401: MLPutFunction(link, "N", 1);
1402: MLPutFunction(link, "Part", 2);
1403: MLPutSymbol(link, "mattML");
1404: MLPutInteger(link, 9);
1405: MLEndPacket(link);
1406: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1407: MLGetReal(link, &ml->zeroTol);
1409: if (ml->numLevels > 0) {
1410: /* ml->factors[level][part][FACT_U] = ml[[7,level,part,1]]
1411: ml->factors[level][part][FACT_DINV] = Divide[1,Select[ml[[7,level,part,2]],(#>ml[[9]])&]]
1412: ml->factors[level][part][FACT_V] = ml[[7,level,part,3]] */
1413: PetscMalloc(ml->numLevels * sizeof(double ***), &ml->factors);
1414: for(level = 0; level < ml->numLevels; level++) {
1415: if (ml->numPartitions[level] == 0) continue;
1416: PetscMalloc(ml->numPartitions[level] * sizeof(double **), &ml->factors[level]);
1417: for(part = 0; part < ml->numPartitions[level]; part++) {
1418: PetscMalloc(NUM_FACT_DIV * sizeof(double *), &ml->factors[level][part]);
1419: /* U, or U^T in LAPACK terms */
1420: MLPutFunction(link, "EvaluatePacket", 1);
1421: MLPutFunction(link, "Part", 5);
1422: MLPutSymbol(link, "mattML");
1423: MLPutInteger(link, 7);
1424: MLPutInteger(link, level+1);
1425: MLPutInteger(link, part+1);
1426: MLPutInteger(link, 1);
1427: MLEndPacket(link);
1428: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1429: MLGetRealArray(link, &U, &Udims, &Uheads, &Udepth);
1430: if (Udepth > 1) {
1431: if (Udepth != 2) SETERRQ(PETSC_ERR_PLIB, "Invalid U matrix");
1432: if ((Udims[0] != ml->numPartitionRows[level][part]) || (Udims[0] != Udims[1])) {
1433: SETERRQ(PETSC_ERR_PLIB, "Incompatible dimensions for U matrix");
1434: }
1435: PetscMalloc(Udims[0]*Udims[0] * sizeof(double), &ml->factors[level][part][FACT_U]);
1436: /* Notice that LAPACK will think that this is U^T, or U in LAPACK terms */
1437: PetscMemcpy(ml->factors[level][part][FACT_U], U, Udims[0]*Udims[0] * sizeof(double));
1438: } else if (ml->numPartitionRows[level][part] != 0) {
1439: SETERRQ(PETSC_ERR_PLIB, "Missing U matrix");
1440: }
1441: MLDisownRealArray(link, U, Udims, Uheads, Udepth);
1442: /* D^{-1} */
1443: MLPutFunction(link, "EvaluatePacket", 1);
1444: MLPutFunction(link, "Divide", 2);
1445: MLPutReal(link, 1.0);
1446: MLPutFunction(link, "Select", 2);
1447: MLPutFunction(link, "Part", 5);
1448: MLPutSymbol(link, "mattML");
1449: MLPutInteger(link, 7);
1450: MLPutInteger(link, level+1);
1451: MLPutInteger(link, part+1);
1452: MLPutInteger(link, 2);
1453: MLPutFunction(link, "Function", 2);
1454: MLPutSymbol(link, "x");
1455: MLPutFunction(link, "Greater", 2);
1456: MLPutSymbol(link, "x");
1457: MLPutFunction(link, "Part", 2);
1458: MLPutSymbol(link, "mattML");
1459: MLPutInteger(link, 9);
1460: MLEndPacket(link);
1461: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1462: MLGetRealList(link, &singVals, &numSingVals);
1463: if (numSingVals > ml->numPartitionCols[level][part]) {
1464: SETERRQ(PETSC_ERR_PLIB, "Invalid factor list in MLData object");
1465: }
1466: if (ml->numPartitionCols[level][part] > 0) {
1467: PetscMalloc(ml->numPartitionCols[level][part] * sizeof(double), &ml->factors[level][part][FACT_DINV]);
1468: PetscMemzero(ml->factors[level][part][FACT_DINV], ml->numPartitionCols[level][part] * sizeof(double));
1469: PetscMemcpy(ml->factors[level][part][FACT_DINV], singVals, numSingVals * sizeof(double));
1470: }
1471: MLDisownRealList(link, singVals, numSingVals);
1472: /* V^T, but V in LAPACK terms */
1473: MLPutFunction(link, "EvaluatePacket", 1);
1474: MLPutFunction(link, "Transpose", 1);
1475: MLPutFunction(link, "Part", 5);
1476: MLPutSymbol(link, "mattML");
1477: MLPutInteger(link, 7);
1478: MLPutInteger(link, level+1);
1479: MLPutInteger(link, part+1);
1480: MLPutInteger(link, 3);
1481: MLEndPacket(link);
1482: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1483: MLGetRealArray(link, &V, &Vdims, &Vheads, &Vdepth);
1484: if (Vdepth > 1) {
1485: if (Vdepth != 2) SETERRQ(PETSC_ERR_PLIB, "Invalid V matrix");
1486: if ((Vdims[0] != ml->numPartitionCols[level][part]) || (Vdims[0] != Vdims[1])) {
1487: SETERRQ(PETSC_ERR_PLIB, "Incompatible dimensions for U matrix");
1488: }
1489: PetscMalloc(Vdims[0]*Vdims[0] * sizeof(double), &ml->factors[level][part][FACT_V]);
1490: /* Notice that LAPACK will think that this is V, or V^T in LAPACK terms */
1491: PetscMemcpy(ml->factors[level][part][FACT_V], V, Vdims[0]*Vdims[0] * sizeof(double));
1492: } else if (ml->numPartitionCols[level][part] != 0) {
1493: SETERRQ(PETSC_ERR_PLIB, "Missing V matrix");
1494: }
1495: MLDisownRealArray(link, V, Vdims, Vheads, Vdepth);
1496: }
1497: }
1499: /* ml->grads = ml[[8]] */
1500: PetscMalloc(ml->numLevels * sizeof(Mat), &ml->grads);
1501: PetscMalloc(ml->numLevels * sizeof(Vec), &ml->bdReduceVecs);
1502: PetscMalloc(ml->numLevels * sizeof(Vec), &ml->colReduceVecs);
1503: PetscMalloc(ml->numLevels * sizeof(Vec), &ml->colReduceVecs2);
1504: for(level = 0; level < ml->numLevels; level++) {
1505: if (ml->numPartitions[level] == 0) continue;
1506: /* numRows = ml[[8,level,1]] */
1507: MLPutFunction(link, "EvaluatePacket", 1);
1508: MLPutFunction(link, "Part", 4);
1509: MLPutSymbol(link, "mattML");
1510: MLPutInteger(link, 8);
1511: MLPutInteger(link, level+1);
1512: MLPutInteger(link, 1);
1513: MLEndPacket(link);
1514: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1515: MLGetInteger(link, &numBdRows);
1516: VecCreateSeq(PETSC_COMM_SELF, numBdRows, &ml->bdReduceVecs[level]);
1517: /* numCols = ml[[8,level,2]] */
1518: MLPutFunction(link, "EvaluatePacket", 1);
1519: MLPutFunction(link, "Part", 4);
1520: MLPutSymbol(link, "mattML");
1521: MLPutInteger(link, 8);
1522: MLPutInteger(link, level+1);
1523: MLPutInteger(link, 2);
1524: MLEndPacket(link);
1525: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1526: MLGetInteger(link, &numBdCols);
1527: VecCreateSeq(PETSC_COMM_SELF, numBdCols, &ml->colReduceVecs[level]);
1528: VecDuplicate(ml->colReduceVecs[level], &ml->colReduceVecs2[level]);
1529: /* nnz = Map[Apply[Subtract,Sort[#,Greater]]&, Partition[ml[[8,level,3]],2,1]] */
1530: MLPutFunction(link, "EvaluatePacket", 1);
1531: MLPutFunction(link, "Map", 2);
1532: MLPutFunction(link, "Function", 2);
1533: MLPutSymbol(link, "x");
1534: MLPutFunction(link, "Apply", 2);
1535: MLPutSymbol(link, "Subtract");
1536: MLPutFunction(link, "Sort", 2);
1537: MLPutSymbol(link, "x");
1538: MLPutSymbol(link, "Greater");
1539: MLPutFunction(link, "Partition", 3);
1540: MLPutFunction(link, "Part", 4);
1541: MLPutSymbol(link, "mattML");
1542: MLPutInteger(link, 8);
1543: MLPutInteger(link, level+1);
1544: MLPutInteger(link, 3);
1545: MLPutInteger(link, 2);
1546: MLPutInteger(link, 1);
1547: MLEndPacket(link);
1548: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1549: MLGetIntegerList(link, &nnz, &len);
1550: if (len != numBdRows) SETERRQ(PETSC_ERR_PLIB, "Invalid boundary gradient matrix");
1551: MatCreateSeqAIJ(PETSC_COMM_SELF, numBdRows, numBdCols, PETSC_DEFAULT, nnz, &ml->grads[level]);
1552: grad = (Mat_SeqAIJ *) ml->grads[level]->data;
1553: PetscMemcpy(grad->ilen, nnz, numBdRows * sizeof(int));
1554: MLDisownIntegerList(link, nnz, len);
1555: /* ml->grads[level]->i = ml[[8,level,3]] */
1556: MLPutFunction(link, "EvaluatePacket", 1);
1557: MLPutFunction(link, "Part", 4);
1558: MLPutSymbol(link, "mattML");
1559: MLPutInteger(link, 8);
1560: MLPutInteger(link, level+1);
1561: MLPutInteger(link, 3);
1562: MLEndPacket(link);
1563: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1564: MLGetIntegerList(link, &offsets, &len);
1565: if (len != numBdRows+1) SETERRQ(PETSC_ERR_PLIB, "Invalid boundary gradient matrix");
1566: for(row = 0; row <= numBdRows; row++)
1567: grad->i[row] = offsets[row]-1;
1568: MLDisownIntegerList(link, offsets, len);
1569: /* ml->grads[level]->j = ml[[8,level,4]] */
1570: MLPutFunction(link, "EvaluatePacket", 1);
1571: MLPutFunction(link, "Part", 4);
1572: MLPutSymbol(link, "mattML");
1573: MLPutInteger(link, 8);
1574: MLPutInteger(link, level+1);
1575: MLPutInteger(link, 4);
1576: MLEndPacket(link);
1577: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1578: MLGetIntegerList(link, &cols, &len);
1579: if (len != grad->i[numBdRows]) SETERRQ(PETSC_ERR_PLIB, "Invalid boundary gradient matrix");
1580: for(col = 0; col < len; col++)
1581: grad->j[col] = cols[col]-1;
1582: MLDisownIntegerList(link, cols, len);
1583: /* ml->grads[level]->i = ml[[8,level,5]] */
1584: MLPutFunction(link, "EvaluatePacket", 1);
1585: MLPutFunction(link, "Part", 4);
1586: MLPutSymbol(link, "mattML");
1587: MLPutInteger(link, 8);
1588: MLPutInteger(link, level+1);
1589: MLPutInteger(link, 5);
1590: MLEndPacket(link);
1591: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1592: MLGetRealList(link, &vals, &len);
1593: if (len != grad->i[numBdRows]) SETERRQ(PETSC_ERR_PLIB, "Invalid boundary gradient matrix");
1594: PetscMemcpy(grad->a, vals, len * sizeof(double));
1595: MLDisownRealList(link, vals, len);
1596: /* Fix up all the members */
1597: grad->nz += len;
1598: MatAssemblyBegin(ml->grads[level], MAT_FINAL_ASSEMBLY);
1599: MatAssemblyEnd(ml->grads[level], MAT_FINAL_ASSEMBLY);
1600: }
1601: } else {
1602: PetscMalloc(1 * sizeof(double ***), &ml->factors);
1603: PetscMalloc(1 * sizeof(Mat), &ml->grads);
1604: PetscMalloc(1 * sizeof(Vec), &ml->bdReduceVecs);
1605: PetscMalloc(1 * sizeof(Vec), &ml->colReduceVecs);
1606: PetscMalloc(1 * sizeof(Vec), &ml->colReduceVecs2);
1607: }
1609: ml->interiorWorkLen = 1;
1610: for(level = 0; level < ml->numLevels; level++) {
1611: for(part = 0; part < ml->numPartitions[level]; part++)
1612: ml->interiorWorkLen = PetscMax(ml->interiorWorkLen, ml->numPartitionRows[level][part]);
1613: }
1614: PetscMalloc(ml->interiorWorkLen * sizeof(double), &ml->interiorWork);
1615: PetscMalloc(ml->interiorWorkLen * sizeof(double), &ml->interiorWork2);
1617: return(0);
1618: #endif
1620: return(0);
1621: }
1623: #if 0
1624: /*------------------------------ Functions for Triangular 2d Meshes -------------------------------------------------*/
1625: #undef __FUNCT__
1627: int PetscViewerMathematicaCreateSamplePoints_Triangular_2D(PetscViewer viewer, GVec v)
1628: {
1629: #ifdef PETSC_HAVE_MATHEMATICA
1630: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
1631: MLINK link = vmath->link; /* The link to Mathematica */
1632: Grid grid;
1633: int numNodes;
1634: int *classes;
1635: int *offsets;
1636: double *array;
1637: int *localStart;
1638: int localOffset, comp;
1639: int node, nclass;
1640: int ierr;
1643: ierr = GVecGetGrid(v, &grid);
1644: numNodes = grid->mesh->numNodes;
1645: comp = grid->viewComp;
1646: offsets = grid->order->offsets;
1647: localStart = grid->order->localStart[grid->viewField];
1648: classes = grid->cm->classes;
1650: /* Attach a value to each point in the mesh -- Extra values are put in for fields not
1651: defined on some nodes, but these values are never used */
1652: VecGetArray(v, &array);
1653: MLPutFunction(link, "ReplaceAll", 2);
1654: MLPutFunction(link, "Thread", 1);
1655: MLPutFunction(link, "f", 2);
1656: MLPutSymbol(link, "nodes");
1657: MLPutFunction(link, "List", numNodes);
1658: for(node = 0; node < numNodes; node++)
1659: {
1660: nclass = classes[node];
1661: localOffset = localStart[nclass] + comp;
1662: MLPutReal(link, array[offsets[node] + localOffset]);
1663: }
1664: MLPutFunction(link, "Rule", 2);
1665: MLPutSymbol(link, "f");
1666: MLPutSymbol(link, "Append");
1667: VecRestoreArray(v, &array);
1669: return(0);
1670: #endif
1672: return(0);
1673: }
1675: #undef __FUNCT__
1677: int PetscViewerMathematicaCreateVectorSamplePoints_Triangular_2D(PetscViewer viewer, GVec v)
1678: {
1679: #ifdef PETSC_HAVE_MATHEMATICA
1680: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
1681: MLINK link = vmath->link; /* The link to Mathematica */
1682: Grid grid;
1683: int numNodes;
1684: int *classes;
1685: int *offsets;
1686: int *fieldClasses;
1687: double *array;
1688: int *localStart;
1689: int localOffset;
1690: int node, nclass;
1691: int ierr;
1694: ierr = GVecGetGrid(v, &grid);
1695: numNodes = grid->mesh->numNodes;
1696: fieldClasses = grid->cm->fieldClasses[grid->viewField];
1697: offsets = grid->order->offsets;
1698: localStart = grid->order->localStart[grid->viewField];
1699: classes = grid->cm->classes;
1701: /* Make a list {{{x_0,y_0},{f^0_x,f^0_y}},...} */
1702: VecGetArray(v, &array);
1703: MLPutFunction(link, "ReplaceAll", 2);
1704: MLPutFunction(link, "Thread", 1);
1705: MLPutFunction(link, "f", 2);
1706: MLPutSymbol(link, "nodes");
1707: MLPutFunction(link, "List", numNodes);
1708: for(node = 0; node < numNodes; node++)
1709: {
1710: nclass = classes[node];
1711: if (fieldClasses[nclass])
1712: {
1713: localOffset = localStart[nclass];
1714: MLPutFunction(link, "List", 2);
1715: MLPutReal(link, array[offsets[node] + localOffset]);
1716: MLPutReal(link, array[offsets[node] + localOffset + 1]);
1717: }
1718: else
1719: {
1720: /* Vectors of length zero are ignored */
1721: MLPutFunction(link, "List", 2);
1722: MLPutReal(link, 0.0);
1723: MLPutReal(link, 0.0);
1724: }
1725: }
1726: MLPutFunction(link, "Rule", 2);
1727: MLPutSymbol(link, "f");
1728: MLPutSymbol(link, "List");
1729: VecRestoreArray(v, &array);
1731: return(0);
1732: #endif
1734: return(0);
1735: }
1737: #undef __FUNCT__
1739: int PetscViewerMathematicaCreateInterpolatedSamplePoints_Triangular_2D(PetscViewer viewer, GVec v, int vComp)
1740: {
1741: #ifdef PETSC_HAVE_MATHEMATICA
1742: #if 0
1743: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
1744: MLINK link = vmath->link; /* The link to Mathematica */
1745: InterpolationContext iCtx;
1746: Grid grid;
1747: Mesh mesh;
1748: double *x, *y, *values;
1749: long m, n;
1750: double startX, endX, incX;
1751: double startY, endY, incY;
1752: int comp;
1753: int proc, row, col;
1754: PetscTruth opt;
1755: int ierr;
1756: #endif
1759: #if 0
1760: ierr = GVecGetGrid(v, &grid);
1761: ierr = GridGetMesh(grid, &mesh);
1762: comp = grid->comp[grid->viewField];
1764: /* This sucks, but I will fix it later (It is for GridReduceInterpolationElementVec_Triangular_2D) */
1765: grid->classesOld = grid->classes;
1767: /* Setup InterpolationContext */
1768: iCtx.vec = v;
1769: iCtx.mesh = mesh;
1770: iCtx.order = grid->order;
1771: iCtx.ghostVec = grid->ghostVec;
1772: iCtx.field = grid->viewField;
1773: iCtx.numProcs = mesh->part->numProcs;
1774: iCtx.rank = mesh->part->rank;
1775: PetscMalloc(iCtx.numProcs * sizeof(int), &iCtx.activeProcs);
1776: PetscMalloc(iCtx.numProcs * sizeof(int), &iCtx.calcProcs);
1777: PetscMalloc(iCtx.numProcs*3 * sizeof(PetscScalar), &iCtx.coords);
1778: PetscMalloc(iCtx.numProcs * sizeof(PetscScalar *), &iCtx.values);
1779: for(proc = 0; proc < iCtx.numProcs; proc++) {
1780: PetscMalloc(comp * sizeof(PetscScalar), &iCtx.values[proc]);
1781: }
1783: /* Setup domain */
1784: startX = 0.0;
1785: startY = 0.0;
1786: endX = 1.0;
1787: endY = 1.0;
1788: ierr = PetscOptionsGetDouble("viewer_", "-math_start_x", &startX, &opt);
1789: ierr = PetscOptionsGetDouble("viewer_", "-math_start_y", &startY, &opt);
1790: ierr = PetscOptionsGetDouble("viewer_", "-math_end_x", &endX, &opt);
1791: ierr = PetscOptionsGetDouble("viewer_", "-math_end_y", &endY, &opt);
1792: ierr = PetscOptionsGetInt("viewer_", "-math_div_x", (int *) &n, &opt);
1793: ierr = PetscOptionsGetInt("viewer_", "-math_div_y", (int *) &m, &opt);
1794: ierr = PetscMalloc((n+1) * sizeof(double), &x);
1795: ierr = PetscMalloc((n+1) * sizeof(double), &y);
1796: ierr = PetscMalloc((n+1)*comp * sizeof(double), &values);
1797: incX = (endX - startX)/n;
1798: incY = (endY - startY)/m;
1800: x[0] = startX;
1801: for(col = 1; col <= n; col++)
1802: x[col] = x[col-1] + incX;
1804: MLPutFunction(link, "List", m+1);
1805: for(row = 0; row <= m; row++)
1806: {
1807: PetscMemzero(values, (n+1)*comp * sizeof(double));
1808: for(col = 0; col <= n; col++)
1809: y[col] = startY + incY*row;
1810: PointFunctionInterpolateField(n+1, comp, x, y, x, values, &iCtx);
1811: if (vComp >= 0)
1812: {
1813: for(col = 0; col <= n; col++)
1814: values[col] = values[col*comp+vComp];
1815: MLPutRealList(link, values, n+1);
1816: }
1817: else
1818: {
1819: MLPutFunction(link, "Transpose", 1);
1820: MLPutFunction(link, "List", 2);
1821: MLPutFunction(link, "Transpose", 1);
1822: MLPutFunction(link, "List", 2);
1823: MLPutRealList(link, x, n+1);
1824: MLPutRealList(link, y, n+1);
1825: MLPutFunction(link, "Partition", 2);
1826: MLPutRealList(link, values, (n+1)*comp);
1827: MLPutInteger(link, comp);
1828: }
1829: }
1831: /* This sucks, but I will fix it later (It is for GridReduceInterpolationElementVec_Triangular_2D) */
1832: grid->classesOld = PETSC_NULL;
1834: /* Cleanup */
1835: PetscFree(x);
1836: PetscFree(y);
1837: PetscFree(values);
1838: PetscFree(iCtx.activeProcs);
1839: PetscFree(iCtx.calcProcs);
1840: PetscFree(iCtx.coords);
1841: for(proc = 0; proc < iCtx.numProcs; proc++) {
1842: PetscFree(iCtx.values[proc]);
1843: }
1844: PetscFree(iCtx.values);
1845: #endif
1847: return(0);
1848: #endif
1850: return(0);
1851: }
1853: #undef __FUNCT__
1855: int PetscViewerMathematica_Mesh_Triangular_2D(PetscViewer viewer, Mesh mesh)
1856: {
1857: #ifdef PETSC_HAVE_MATHEMATICA
1858: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
1859: MLINK link = vmath->link;
1860: Mesh_Triangular *tri = (Mesh_Triangular *) mesh->data;
1861: int numCorners = mesh->numCorners;
1862: int numFaces = mesh->numFaces;
1863: int *faces = tri->faces;
1864: int numNodes = mesh->numNodes;
1865: double *nodes = tri->nodes;
1866: int node, face, corner;
1867: int ierr;
1870: /* Load package to view graphics in a window */
1871: PetscViewerMathematicaLoadGraphics(viewer, vmath->graphicsType);
1873: /* Send the node coordinates */
1874: MLPutFunction(link, "EvaluatePacket", 1);
1875: MLPutFunction(link, "Set", 2);
1876: MLPutSymbol(link, "nodes");
1877: MLPutFunction(link, "List", numNodes);
1878: for(node = 0; node < numNodes; node++) {
1879: MLPutFunction(link, "List", 2);
1880: MLPutReal(link, nodes[node*2]);
1881: MLPutReal(link, nodes[node*2+1]);
1882: }
1883: MLEndPacket(link);
1884: /* Skip packets until ReturnPacket */
1885: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1886: /* Skip ReturnPacket */
1887: MLNewPacket(link);
1889: /* Send the faces */
1890: MLPutFunction(link, "EvaluatePacket", 1);
1891: MLPutFunction(link, "Set", 2);
1892: MLPutSymbol(link, "faces");
1893: MLPutFunction(link, "List", numFaces);
1894: for(face = 0; face < numFaces; face++) {
1895: MLPutFunction(link, "List", numCorners);
1896: for(corner = 0; corner < numCorners; corner++) {
1897: MLPutReal(link, faces[face*numCorners+corner]);
1898: }
1899: }
1900: MLEndPacket(link);
1901: /* Skip packets until ReturnPacket */
1902: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1903: /* Skip ReturnPacket */
1904: MLNewPacket(link);
1906: return(0);
1907: #endif
1909: return(0);
1910: }
1912: #undef __FUNCT__
1914: int PetscViewerMathematicaCheckMesh_Triangular_2D(PetscViewer viewer, Mesh mesh)
1915: {
1916: #ifdef PETSC_HAVE_MATHEMATICA
1917: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
1918: MLINK link = vmath->link;
1919: int numCorners = mesh->numCorners;
1920: int numFaces = mesh->numFaces;
1921: int numNodes = mesh->numNodes;
1922: const char *symbol;
1923: long args;
1924: int dim, type;
1925: PetscTruth match;
1926: int ierr;
1929: /* Check that nodes are defined */
1930: MLPutFunction(link, "EvaluatePacket", 1);
1931: MLPutFunction(link, "ValueQ", 1);
1932: MLPutSymbol(link, "nodes");
1933: MLEndPacket(link);
1934: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1935: MLGetSymbol(link, &symbol);
1936: PetscStrcmp("True", (char *) symbol, &match);
1937: if (match == PETSC_FALSE) {
1938: MLDisownSymbol(link, symbol);
1939: goto redefineMesh;
1940: }
1941: MLDisownSymbol(link, symbol);
1943: /* Check the dimensions */
1944: MLPutFunction(link, "EvaluatePacket", 1);
1945: MLPutFunction(link, "MatrixQ", 2);
1946: MLPutSymbol(link, "nodes");
1947: MLPutSymbol(link, "NumberQ");
1948: MLEndPacket(link);
1949: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1950: MLGetSymbol(link, &symbol);
1951: PetscStrcmp("True", (char *) symbol, &match);
1952: if (match == PETSC_FALSE) {
1953: MLDisownSymbol(link, symbol);
1954: goto redefineMesh;
1955: }
1956: MLDisownSymbol(link, symbol);
1957: MLPutFunction(link, "EvaluatePacket", 1);
1958: MLPutFunction(link, "Dimensions", 1);
1959: MLPutSymbol(link, "nodes");
1960: MLEndPacket(link);
1961: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1962: args = 0;
1963: type = MLGetNext(link);
1964: MLGetArgCount(link, &args);
1965: if (args != 2) {
1966: MLNewPacket(link);
1967: goto redefineMesh;
1968: }
1969: MLGetSymbol(link, &symbol);
1970: MLDisownSymbol(link, symbol);
1971: MLGetInteger(link, &dim);
1972: if (dim != numNodes) {
1973: MLNewPacket(link);
1974: goto redefineMesh;
1975: }
1976: MLGetInteger(link, &dim);
1977: if (dim != mesh->dim) {
1978: MLNewPacket(link);
1979: goto redefineMesh;
1980: }
1982: /* Check that faces are defined */
1983: MLPutFunction(link, "EvaluatePacket", 1);
1984: MLPutFunction(link, "ValueQ", 1);
1985: MLPutSymbol(link, "faces");
1986: MLEndPacket(link);
1987: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
1988: MLGetSymbol(link, &symbol);
1989: PetscStrcmp("True", (char *) symbol, &match);
1990: if (match == PETSC_FALSE) {
1991: MLDisownSymbol(link, symbol);
1992: goto redefineMesh;
1993: }
1994: MLDisownSymbol(link, symbol);
1996: /* Check the dimensions */
1997: MLPutFunction(link, "EvaluatePacket", 1);
1998: MLPutFunction(link, "MatrixQ", 2);
1999: MLPutSymbol(link, "faces");
2000: MLPutSymbol(link, "NumberQ");
2001: MLEndPacket(link);
2002: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
2003: MLGetSymbol(link, &symbol);
2004: PetscStrcmp("True", (char *) symbol, &match);
2005: if (match == PETSC_FALSE) {
2006: MLDisownSymbol(link, symbol);
2007: goto redefineMesh;
2008: }
2009: MLDisownSymbol(link, symbol);
2010: MLPutFunction(link, "EvaluatePacket", 1);
2011: MLPutFunction(link, "Dimensions", 1);
2012: MLPutSymbol(link, "faces");
2013: MLEndPacket(link);
2014: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
2015: args = 0;
2016: type = MLGetNext(link);
2017: MLGetArgCount(link, &args);
2018: if (args != 2) {
2019: MLNewPacket(link);
2020: goto redefineMesh;
2021: }
2022: MLGetSymbol(link, &symbol);
2023: MLDisownSymbol(link, symbol);
2024: MLGetInteger(link, &dim);
2025: if (dim != numFaces) {
2026: MLNewPacket(link);
2027: goto redefineMesh;
2028: }
2029: MLGetInteger(link, &dim);
2030: if (dim != numCorners) {
2031: MLNewPacket(link);
2032: goto redefineMesh;
2033: }
2035: return(0);
2037: redefineMesh:
2038: MeshView(mesh, viewer);
2039: return(0);
2040: #endif
2042: return(0);
2043: }
2045: #undef __FUNCT__
2047: int PetscViewerMathematicaTriangulationPlot_Triangular_2D(PetscViewer viewer, GVec v)
2048: {
2049: #ifdef PETSC_HAVE_MATHEMATICA
2050: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
2051: MLINK link = vmath->link; /* The link to Mathematica */
2052: Mesh mesh;
2053: Grid grid;
2054: int numCorners;
2055: int ierr;
2058: ierr = GVecGetGrid(v, &grid);
2059: mesh = grid->mesh;
2060: numCorners = mesh->numCorners;
2062: MLPutFunction(link, "Show", 2);
2063: MLPutFunction(link, "Graphics3D", 1);
2064: MLPutFunction(link, "Flatten", 1);
2065: MLPutFunction(link, "Map", 2);
2066: MLPutFunction(link, "Function", 1);
2067: if ((numCorners == 6) && ((grid->cm->numClasses == 1) || (grid->cm->fieldClasses[grid->viewField][1])))
2068: {
2069: MLPutFunction(link, "List", 4);
2070: /* Triangle 0--5--4 */
2071: MLPutFunction(link, "Polygon", 1);
2072: MLPutFunction(link, "Part", 2);
2073: MLPutSymbol(link, "points");
2074: MLPutFunction(link, "Plus", 2);
2075: MLPutFunction(link, "Part", 2);
2076: MLPutFunction(link, "Slot", 1);
2077: MLPutInteger(link, 1);
2078: MLPutFunction(link, "List", 3);
2079: MLPutInteger(link, 1);
2080: MLPutInteger(link, 6);
2081: MLPutInteger(link, 5);
2082: MLPutInteger(link, 1);
2083: /* Triangle 1--3--5 */
2084: MLPutFunction(link, "Polygon", 1);
2085: MLPutFunction(link, "Part", 2);
2086: MLPutSymbol(link, "points");
2087: MLPutFunction(link, "Plus", 2);
2088: MLPutFunction(link, "Part", 2);
2089: MLPutFunction(link, "Slot", 1);
2090: MLPutInteger(link, 1);
2091: MLPutFunction(link, "List", 3);
2092: MLPutInteger(link, 2);
2093: MLPutInteger(link, 4);
2094: MLPutInteger(link, 6);
2095: MLPutInteger(link, 1);
2096: /* Triangle 2--4--3 */
2097: MLPutFunction(link, "Polygon", 1);
2098: MLPutFunction(link, "Part", 2);
2099: MLPutSymbol(link, "points");
2100: MLPutFunction(link, "Plus", 2);
2101: MLPutFunction(link, "Part", 2);
2102: MLPutFunction(link, "Slot", 1);
2103: MLPutInteger(link, 1);
2104: MLPutFunction(link, "List", 3);
2105: MLPutInteger(link, 3);
2106: MLPutInteger(link, 5);
2107: MLPutInteger(link, 4);
2108: MLPutInteger(link, 1);
2109: /* Triangle 3--4--5 */
2110: MLPutFunction(link, "Polygon", 1);
2111: MLPutFunction(link, "Part", 2);
2112: MLPutSymbol(link, "points");
2113: MLPutFunction(link, "Plus", 2);
2114: MLPutFunction(link, "Part", 2);
2115: MLPutFunction(link, "Slot", 1);
2116: MLPutInteger(link, 1);
2117: MLPutFunction(link, "List", 3);
2118: MLPutInteger(link, 4);
2119: MLPutInteger(link, 5);
2120: MLPutInteger(link, 6);
2121: MLPutInteger(link, 1);
2122: }
2123: else if ((numCorners == 3) || (!grid->cm->fieldClasses[grid->viewField][1]))
2124: {
2125: /* Triangle 0--1--2 */
2126: MLPutFunction(link, "Polygon", 1);
2127: MLPutFunction(link, "Part", 2);
2128: MLPutSymbol(link, "points");
2129: MLPutFunction(link, "Plus", 2);
2130: MLPutFunction(link, "Part", 2);
2131: MLPutFunction(link, "Slot", 1);
2132: MLPutInteger(link, 1);
2133: MLPutFunction(link, "List", 3);
2134: MLPutInteger(link, 1);
2135: MLPutInteger(link, 2);
2136: MLPutInteger(link, 3);
2137: MLPutInteger(link, 1);
2138: } else {
2139: SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid number of local nodes");
2140: }
2141: MLPutSymbol(link, "faces");
2142: MLPutFunction(link, "Rule", 2);
2143: MLPutSymbol(link, "AspectRatio");
2144: MLPutReal(link, 1.0);
2145: return(0);
2146: #endif
2148: return(0);
2149: }
2151: #undef __FUNCT__
2153: int PetscViewerMathematicaVectorPlot_Triangular_2D(PetscViewer viewer, GVec v)
2154: {
2155: #ifdef PETSC_HAVE_MATHEMATICA
2156: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
2157: MLINK link = vmath->link; /* The link to Mathematica */
2158: Grid grid;
2159: Mesh mesh;
2160: PetscReal scale;
2161: PetscTruth opt;
2162: int ierr;
2165: GVecGetGrid(v, &grid);
2166: GridGetMesh(grid, &mesh);
2168: MLPutFunction(link, "ListPlotVectorField", 3);
2169: MLPutSymbol(link, "points");
2170: MLPutFunction(link, "Rule", 2);
2171: MLPutSymbol(link, "AspectRatio");
2172: MLPutReal(link, mesh->sizeY/mesh->sizeX);
2173: MLPutFunction(link, "Rule", 2);
2174: MLPutSymbol(link, "ScaleFactor");
2175: PetscOptionsGetReal("viewer_", "-math_scale", &scale, &opt);
2176: if (opt == PETSC_TRUE) {
2177: MLPutReal(link, scale);
2178: } else {
2179: MLPutSymbol(link, "None");
2180: }
2181: return(0);
2182: #endif
2184: return(0);
2185: }
2187: #undef __FUNCT__
2189: int PetscViewerMathematicaSurfacePlot_Triangular_2D(PetscViewer viewer, GVec v)
2190: {
2191: #ifdef PETSC_HAVE_MATHEMATICA
2192: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
2193: MLINK link = vmath->link; /* The link to Mathematica */
2194: PetscReal startX, endX;
2195: PetscReal startY, endY;
2196: PetscTruth opt;
2197: int ierr;
2200: /* Setup domain */
2201: startX = 0.0;
2202: startY = 0.0;
2203: endX = 1.0;
2204: endY = 1.0;
2205: ierr = PetscOptionsGetReal("viewer_", "-math_start_x", &startX, &opt);
2206: ierr = PetscOptionsGetReal("viewer_", "-math_start_y", &startY, &opt);
2207: ierr = PetscOptionsGetReal("viewer_", "-math_end_x", &endX, &opt);
2208: ierr = PetscOptionsGetReal("viewer_", "-math_end_y", &endY, &opt);
2210: MLPutFunction(link, "Show", 1);
2211: MLPutFunction(link, "SurfaceGraphics", 6);
2212: MLPutSymbol(link, "points");
2213: MLPutFunction(link, "Rule", 2);
2214: MLPutSymbol(link, "ClipFill");
2215: MLPutSymbol(link, "None");
2216: MLPutFunction(link, "Rule", 2);
2217: MLPutSymbol(link, "Axes");
2218: MLPutSymbol(link, "True");
2219: MLPutFunction(link, "Rule", 2);
2220: MLPutSymbol(link, "PlotLabel");
2221: MLPutString(link, vmath->objName);
2222: MLPutFunction(link, "Rule", 2);
2223: MLPutSymbol(link, "MeshRange");
2224: MLPutFunction(link, "List", 2);
2225: MLPutFunction(link, "List", 2);
2226: MLPutReal(link, startX);
2227: MLPutReal(link, endX);
2228: MLPutFunction(link, "List", 2);
2229: MLPutReal(link, startY);
2230: MLPutReal(link, endY);
2231: MLPutFunction(link, "Rule", 2);
2232: MLPutSymbol(link, "AspectRatio");
2233: MLPutReal(link, (endY - startY)/(endX - startX));
2234: return(0);
2235: #endif
2237: return(0);
2238: }
2240: #undef __FUNCT__
2242: int PetscViewerMathematica_GVec_Triangular_2D(PetscViewer viewer, GVec v)
2243: {
2244: PetscViewer_Mathematica *vmath = (PetscViewer_Mathematica *) viewer->data;
2245: #ifdef PETSC_HAVE_MATHEMATICA
2246: MLINK link = vmath->link; /* The link to Mathematica */
2247: Mesh mesh;
2248: Grid grid;
2249: Mat P;
2250: GVec w;
2251: int numCorners;
2252: int ierr;
2255: ierr = GVecGetGrid(v, &grid);
2256: mesh = grid->mesh;
2257: numCorners = mesh->numCorners;
2259: /* Check that a field component has been specified */
2260: if ((grid->viewField < 0) || (grid->viewField >= grid->numFields)) return(0);
2262: if (grid->isConstrained) {
2263: GVecCreate(grid, &w);
2264: GridGetConstraintMatrix(grid, &P);
2265: MatMult(P, v, w);
2266: } else {
2267: w = v;
2268: }
2270: /* Check that the mesh is defined correctly */
2271: PetscViewerMathematicaCheckMesh_Triangular_2D(viewer, mesh);
2273: /* Send the first component of the first field */
2274: MLPutFunction(link, "EvaluatePacket", 1);
2275: switch(vmath->plotType)
2276: {
2277: case MATHEMATICA_TRIANGULATION_PLOT:
2278: MLPutFunction(link, "Module", 2);
2279: /* Make temporary points with each value of the field component in the vector */
2280: MLPutFunction(link, "List", 2);
2281: MLPutSymbol(link, "f");
2282: MLPutFunction(link, "Set", 2);
2283: MLPutSymbol(link, "points");
2284: PetscViewerMathematicaCreateSamplePoints_Triangular_2D(viewer, w);
2285: /* Display the picture */
2286: PetscViewerMathematicaTriangulationPlot_Triangular_2D(viewer, w);
2287: break;
2288: case MATHEMATICA_VECTOR_TRIANGULATION_PLOT:
2289: if (grid->fields[grid->viewField].numComp != 2) {
2290: SETERRQ(PETSC_ERR_ARG_WRONG, "Field must be a 2D vector field for this plot type");
2291: }
2292: MLPutFunction(link, "Module", 2);
2293: /* Make temporary list with points and 2D vector field values */
2294: MLPutFunction(link, "List", 2);
2295: MLPutSymbol(link, "f");
2296: MLPutFunction(link, "Set", 2);
2297: MLPutSymbol(link, "points");
2298: PetscViewerMathematicaCreateVectorSamplePoints_Triangular_2D(viewer, w);
2299: /* Display the picture */
2300: PetscViewerMathematicaVectorPlot_Triangular_2D(viewer, w);
2301: break;
2302: case MATHEMATICA_VECTOR_PLOT:
2303: if (grid->fields[grid->viewField].numComp != 2) {
2304: SETERRQ(PETSC_ERR_ARG_WRONG, "Field must be a 2D vector field for this plot type");
2305: }
2306: MLPutFunction(link, "Module", 2);
2307: /* Make temporary list with points and 2D vector field values */
2308: MLPutFunction(link, "List", 2);
2309: MLPutSymbol(link, "f");
2310: MLPutFunction(link, "Set", 2);
2311: MLPutSymbol(link, "points");
2312: PetscViewerMathematicaCreateInterpolatedSamplePoints_Triangular_2D(viewer, w, -1);
2313: /* Display the picture */
2314: PetscViewerMathematicaVectorPlot_Triangular_2D(viewer, w);
2315: break;
2316: case MATHEMATICA_SURFACE_PLOT:
2317: if (grid->fields[grid->viewField].numComp != 2) {
2318: SETERRQ(PETSC_ERR_ARG_WRONG, "Field must be a 2D vector field for this plot type");
2319: }
2320: MLPutFunction(link, "Module", 2);
2321: /* Make temporary list with interpolated field values on a square mesh */
2322: MLPutFunction(link, "List", 1);
2323: MLPutFunction(link, "Set", 2);
2324: MLPutSymbol(link, "points");
2325: PetscViewerMathematicaCreateInterpolatedSamplePoints_Triangular_2D(viewer, w, grid->viewComp);
2326: /* Display the picture */
2327: PetscViewerMathematicaSurfacePlot_Triangular_2D(viewer, w);
2328: break;
2329: default:
2330: SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid plot type");
2331: }
2332: MLEndPacket(link);
2333: /* Skip packets until ReturnPacket */
2334: PetscViewerMathematicaSkipPackets(viewer, RETURNPKT);
2335: /* Skip ReturnPacket */
2336: MLNewPacket(link);
2338: if (grid->isConstrained) {
2339: VecDestroy(w);
2340: }
2341: #else
2343: switch(vmath->plotType)
2344: {
2345: case MATHEMATICA_TRIANGULATION_PLOT:
2346: break;
2347: case MATHEMATICA_VECTOR_TRIANGULATION_PLOT:
2348: break;
2349: case MATHEMATICA_VECTOR_PLOT:
2350: break;
2351: case MATHEMATICA_SURFACE_PLOT:
2352: break;
2353: default:
2354: SETERRQ(PETSC_ERR_ARG_WRONG, "Invalid plot type");
2355: }
2356: #endif
2357: return(0);
2358: }
2359: #endif