A new directory for ParFUM example programs, with its own simple2D example which...
authorIsaac Dooley <idooley2@illinois.edu>
Fri, 9 Dec 2005 18:11:54 +0000 (18:11 +0000)
committerIsaac Dooley <idooley2@illinois.edu>
Fri, 9 Dec 2005 18:11:54 +0000 (18:11 +0000)
examples/ParFUM/simple2D/Makefile [new file with mode: 0644]
examples/ParFUM/simple2D/README [new file with mode: 0644]
examples/ParFUM/simple2D/cst_nl.C [new file with mode: 0644]
examples/ParFUM/simple2D/pgm.C [new file with mode: 0644]
examples/ParFUM/simple2D/pgm.h [new file with mode: 0644]
examples/ParFUM/simple2D/vector2d.h [new file with mode: 0644]
examples/ParFUM/simple2D/xxx.1.ele [new file with mode: 0644]
examples/ParFUM/simple2D/xxx.1.node [new file with mode: 0644]
examples/ParFUM/simple2D/xxx.poly [new file with mode: 0644]

diff --git a/examples/ParFUM/simple2D/Makefile b/examples/ParFUM/simple2D/Makefile
new file mode 100644 (file)
index 0000000..67fab12
--- /dev/null
@@ -0,0 +1,22 @@
+CHARMC=../../../bin/charmc $(OPTS) 
+
+all: pgm
+
+pgm: pgm.o cst_nl.o
+       $(CHARMC) -o pgm pgm.o cst_nl.o -language ParFUM -module netfem 
+
+pgm.o: pgm.C
+       $(CHARMC) -c pgm.C
+
+cst_nl.o: cst_nl.C
+       $(CHARMC) -c cst_nl.C
+
+
+test: pgm
+       ./charmrun ./pgm +vp4 +p2
+
+bgtest: pgm
+       ./charmrun ./pgm +vp4 +p2 +x2 +y2 +z1
+
+clean:
+       rm -f pgm fpgm *.o conv-host charmrun
diff --git a/examples/ParFUM/simple2D/README b/examples/ParFUM/simple2D/README
new file mode 100644 (file)
index 0000000..4070a13
--- /dev/null
@@ -0,0 +1,50 @@
+This is a simple 2D structural dynamics program
+parallelized using the Charm++ FEM Framework.
+
+
+INPUT
+
+Several very important values like the material 
+properties, boundary conditions, and timestep 
+are hardcoded, which is very silly.
+
+The input mesh for this program is read from 
+files named "xxx.1.ele" and "xxx.1.node".
+The mesh format is compatible with Shewchuk's
+"triangle" 2D meshing program, part of the 
+Archimedes toolkit:
+       http://www-2.cs.cmu.edu/~quake/triangle.html
+
+For example, to generate a new mesh, use:
+       triangle -pqVAa xxx.poly
+
+
+PROCESSING
+
+Time loop:
+  Compute internal force vector for linear-strain triangles
+     in cst_nl.C
+  Communicate to sum internal force vector with other partitions.
+  Advance nodes based on net force and velocity.
+  Occasionally migrate or output results.
+
+
+OUTPUT
+
+This program exports its solution data via NetFEM.
+You can run the program so NetFEM will connect to it
+like:
+       ./charmrun ./pgm ++server ++server-port 1234 +p4
+You'd then connect the NetFEM client to yourhostname:1234.
+
+
+CREDITS
+
+The physics code in simple2D is derived from f90 codes
+by Philippe Geubelle and Michael Scot Breitenfeld.
+
+The translation to C and rather silly "simplification"
+were performed by Orion Sky Lawlor, olawlor@acm.org,
+at the University of Illinois at Urbana-Champaign
+Parallel Programming Lab.  Any errors introduced are mine.
+
diff --git a/examples/ParFUM/simple2D/cst_nl.C b/examples/ParFUM/simple2D/cst_nl.C
new file mode 100644 (file)
index 0000000..f98acd4
--- /dev/null
@@ -0,0 +1,85 @@
+#include "pgm.h"
+
+//Compute forces on constant-strain triangles:
+void CST_NL(const vector2d *coor,const connRec *lm,vector2d *R_net,
+           const vector2d *d,const double *c,
+           int numnp,int numel,
+           double *S11o,double *S22o,double *S12o)
+{
+  int n1,n2,n3,i;
+  double S11,S22,S12,u1,u2,u3,v1,v2,v3,x21,y21,x31,y31,x32,y32;
+  double aa,B1,B2,B3,B4,B5,B6,dudx,dvdy,dudy,dvdx;
+  double E11,E22,E12;
+
+  for (i=0;i<numel;i++) {
+    n1=lm[i][0];
+    n2=lm[i][1];
+    n3=lm[i][2];
+          u1 = d[n1].x;
+          u2 = d[n2].x;
+          u3 = d[n3].x;
+          v1 = d[n1].y;
+          v2 = d[n2].y;
+          v3 = d[n3].y;
+
+          x21 = coor[n2].x-coor[n1].x;
+          y21 = coor[n2].y-coor[n1].y;
+          x31 = coor[n3].x-coor[n1].x;
+          y31 = coor[n3].y-coor[n1].y;
+          x32 = coor[n3].x-coor[n2].x;
+          y32 = coor[n3].y-coor[n2].y;
+
+          aa = x21*y31-x31*y21;
+          B1 = -y32/aa;
+          B2 = x32/aa;
+          B3 = y31/aa;
+          B4 = -x31/aa;
+          B5 = -y21/aa;
+          B6 = x21/aa;
+
+          dudx = B1*u1 + B3*u2 + B5*u3;
+          dvdy = B2*v1 + B4*v2 + B6*v3;
+          dudy = B2*u1 + B4*u2 + B6*u3;
+          dvdx = B1*v1 + B3*v2 + B5*v3;
+          E11 = dudx + 0.5*(dudx*dudx + dvdx*dvdx);
+          E22 = dvdy + 0.5*(dvdy*dvdy + dudy*dudy);
+          E12 = dudy + dvdx + dudx*dudy + dvdy*dvdx;
+
+          // Calculate CST stresses
+          S11 = E11*c[0] + E22*c[1];
+          S22 = E11*c[1] + E22*c[2];
+          S12 = E12*c[3];
+          S11o[i]=S11;
+          S22o[i]=S22;
+          S12o[i]=S12;
+         
+          // Calculate R-internal force vector
+          R_net[n1] -= aa*0.5*vector2d(
+               S11*B1*(1.0+dudx) +                 
+               S22*B2*dudy +                        
+               S12*(B2*(1.0+dudx) + B1*dudy)
+           ,
+               S11*B1*dvdx +                        
+               S22*B2*(1.0+dvdy) +                 
+               S12*(B1*(1.0+dvdy)+B2*dvdx)
+           );
+          R_net[n2] -= aa*0.5*vector2d(   
+               S11*B3*(1.0+dudx) +                 
+               S22*B4*dudy +                        
+               S12*(B4*(1.0+dudx) + B3*dudy)
+           ,
+              S11*B3*dvdx +                        
+               S22*B4*(1.0+dvdy) +                 
+               S12*(B3*(1.0+dvdy)+B4*dvdx)
+           );
+          R_net[n3] -= aa*0.5*vector2d(   
+               S11*B5*(1.0+dudx) +                 
+               S22*B6*dudy +                        
+               S12*(B6*(1.0+dudx) + B5*dudy)
+           ,
+               S11*B5*dvdx +                        
+               S22*B6*(1.0+dvdy) +                 
+               S12*(B5*(1.0+dvdy)+B6*dvdx)
+           ); 
+  }
+}
diff --git a/examples/ParFUM/simple2D/pgm.C b/examples/ParFUM/simple2D/pgm.C
new file mode 100644 (file)
index 0000000..e9d7b87
--- /dev/null
@@ -0,0 +1,329 @@
+/*
+ Charm++ Finite-Element Framework Program:
+   Performs simple 2D structural simulation on Triangle-style inputs.
+   
+ init: 
+    read node input file
+    pass nodes into ParFUM framework
+    read element input file
+    pass elements into ParFUM framework
+ driver:
+    extract mesh chunk from framework
+    calculate masses of nodes
+    timeloop
+      compute forces within my chunk
+      communicate
+      apply boundary conditions and integrate
+      pass data to NetFEM
+ Among the hideous things about this program are:
+   -Hardcoded material properties, timestep, and boundary conditions
+   -Limited to 2D
+
+ Converted from f90 Structural Materials program by 
+       Orion Sky Lawlor, 2001
+
+ Updated to new ParFUM interface by
+    Isaac Dooley, 2005
+
+ */
+
+#include "pgm.h"
+
+
+extern "C" void
+init(void)
+{
+  CkPrintf("init started\n");
+  double startTime=CmiWallTimer();
+  const char *eleName="xxx.1.ele";
+  const char *nodeName="xxx.1.node";
+  int nPts=0; //Number of nodes
+  vector2d *pts=0; //Node coordinates
+
+  CkPrintf("Reading node coordinates from %s\n",nodeName);
+  //Open and read the node coordinate file
+  {
+    char line[1024];
+    FILE *f=fopen(nodeName,"r");
+    if (f==NULL) die("Can't open node file!");
+    fgets(line,1024,f);
+    if (1!=sscanf(line,"%d",&nPts)) die("Can't read number of points!");
+    pts=new vector2d[nPts];
+    for (int i=0;i<nPts;i++) {
+      int ptNo;
+      if (NULL==fgets(line,1024,f)) die("Can't read node input line!");
+      if (3!=sscanf(line,"%d%lf%lf",&ptNo,&pts[i].x,&pts[i].y)) 
+       die("Can't parse node input line!");
+    }
+    fclose(f);
+  }
+  int nEle=0;
+  connRec *ele=NULL;
+  CkPrintf("Reading elements from %s\n",eleName);
+  //Open and read the element connectivity file
+  {
+    char line[1024];
+    FILE *f=fopen(eleName,"r");
+    if (f==NULL) die("Can't open element file!");
+    fgets(line,1024,f);
+    if (1!=sscanf(line,"%d",&nEle)) die("Can't read number of elements!");
+    ele=new connRec[nEle];
+    for (int i=0;i<nEle;i++) {
+      int elNo;
+      if (NULL==fgets(line,1024,f)) die("Can't read element input line!");
+      if (4!=sscanf(line,"%d%d%d%d",&elNo,&ele[i][0],&ele[i][1],&ele[i][2])) 
+       die("Can't parse element input line!");  
+      ele[i][0]--; //Fortran to C indexing
+      ele[i][1]--; //Fortran to C indexing
+      ele[i][2]--; //Fortran to C indexing
+    }
+    fclose(f);
+  }
+
+
+  int mesh=FEM_Mesh_default_write(); // Tell framework we are writing to the mesh
+
+  CkPrintf("Passing node coords to framework\n");
+
+  /*   Old versions used FEM_Set_node() and FEM_Set_node_data()
+   *   New versions use the more flexible FEM_Set_Data()
+   */
+
+  FEM_Mesh_data(mesh,        // Add nodes to the current mesh
+                FEM_NODE,        // We are registering nodes
+                FEM_DATA+0,      // Register the point locations which are normally 
+                                 // the first data elements for an FEM_NODE
+                (double *)pts,   // The array of point locations
+                0,               // 0 based indexing
+                nPts,            // The number of points
+                FEM_DOUBLE,      // Coordinates are doubles
+                2);              // Points have dimension 2 (x,y)
+
+  CkPrintf("Passing elements to framework\n");
+
+  /*   Old versions used FEM_Set_elem() and FEM_Set_elem_conn() 
+   *   New versions use the more flexible FEM_Set_Data()
+   */
+
+  FEM_Mesh_data(mesh,      // Add nodes to the current mesh
+                FEM_ELEM+0,      // We are registering elements with type 0
+                                 // The next type of element could be registered with FEM_ELEM+1
+                FEM_CONN,        // Register the connectivity table for this
+                                 // data elements for this type of FEM entity
+                (int *)ele,      // The array of point locations
+                0,               // 0 based indexing
+                nEle,            // The number of elements
+                FEM_INDEX_0,     // We use zero based node numbering
+                3);              // Elements have degree 3, since triangles are defined 
+                                 // by three nodes
+   
+  delete[] ele;
+  delete[] pts;
+
+  CkPrintf("Finished with init (Reading took %.f s)\n",CmiWallTimer()-startTime);
+
+}
+
+
+// A driver() function 
+// driver() is required in all FEM programs
+extern "C" void
+driver(void)
+{
+  int nnodes,nelems,ignored;
+  int i, myId=FEM_My_partition();
+  myGlobals g;
+  FEM_Register(&g,(FEM_PupFn)pup_myGlobals);
+  
+  
+  int mesh=FEM_Mesh_default_read(); // Tell framework we are reading data from the mesh
+  
+  // Get node data
+  nnodes=FEM_Mesh_get_length(mesh,FEM_NODE); // Get number of nodes
+  g.coord=new vector2d[nnodes];
+  // Get node positions
+  FEM_Mesh_data(mesh, FEM_NODE, FEM_DATA+0, (double*)g.coord, 0, nnodes, FEM_DOUBLE, 2);  
+
+
+  // Get element data
+  nelems=FEM_Mesh_get_length(mesh,FEM_ELEM+0); // Get number of elements
+  g.nnodes=nnodes; g.nelems=nelems;
+  g.conn=new connRec[nelems];
+  g.S11=new double[nelems];
+  g.S22=new double[nelems];
+  g.S12=new double[nelems];
+  // Get connectivity for elements
+  FEM_Mesh_data(mesh, FEM_ELEM+0, FEM_CONN, (int *)g.conn, 0, nelems, FEM_INDEX_0, 3);  
+
+
+  //Initialize associated data
+  g.R_net=new vector2d[nnodes]; //Net force
+  g.d=new vector2d[nnodes];//Node displacement
+  g.v=new vector2d[nnodes];//Node velocity
+  g.a=new vector2d[nnodes];//Node acceleration
+  for (i=0;i<nnodes;i++)
+    g.R_net[i]=g.d[i]=g.v[i]=g.a[i]=vector2d(0.0);
+
+//Apply a small initial perturbation to positions
+  for (i=0;i<nnodes;i++) {
+         const double max=1.0e-10/15.0; //Tiny perturbation
+         g.d[i].x+=max*(i&15);
+         g.d[i].y+=max*((i+5)&15);
+  }
+
+  int fid=FEM_Create_simple_field(FEM_DOUBLE,2);
+
+  //Timeloop
+  if (myId==0)
+    CkPrintf("Entering timeloop\n");
+  int tSteps=5000;
+  double startTime, totalStart;
+  startTime=totalStart=CkWallTimer();
+  for (int t=0;t<tSteps;t++) {
+    if (1) { //Structural mechanics
+
+    //Compute forces on nodes exerted by elements
+       CST_NL(g.coord,g.conn,g.R_net,g.d,matConst,nnodes,nelems,g.S11,g.S22,g.S12);
+
+    //Communicate net force on shared nodes
+       FEM_Update_field(fid,g.R_net);
+
+    //Advance node positions
+       advanceNodes(dt,nnodes,g.coord,g.R_net,g.a,g.v,g.d,0);
+    }
+
+    //Debugging/perf. output
+    double curTime=CkWallTimer();
+    double total=curTime-startTime;
+    startTime=curTime;
+    if (myId==0 && (t%64==0)) {
+           CkPrintf("%d %.6f sec for loop %d \n",CkNumPes(),total,t);
+           if (0) {
+             CkPrintf("    Triangle 0:\n");
+             for (int j=0;j<3;j++) {
+                   int n=g.conn[0][j];
+                   CkPrintf("    Node %d: coord=(%.4f,%.4f)  d=(%.4g,%.4g)\n",
+                            n,g.coord[n].x,g.coord[n].y,g.d[n].x,g.d[n].y);
+             }
+           }
+    }
+    /* perform migration-based load balancing */
+    if (t%1024==0)
+      FEM_Migrate();
+    
+    if (t%1024==0) { //Publish data to the net
+           NetFEM n=NetFEM_Begin(FEM_My_partition(),t,2,NetFEM_POINTAT);
+           
+           NetFEM_Nodes(n,nnodes,(double *)g.coord,"Position (m)");
+           NetFEM_Vector(n,(double *)g.d,"Displacement (m)");
+           NetFEM_Vector(n,(double *)g.v,"Velocity (m/s)");
+           
+           NetFEM_Elements(n,nelems,3,(int *)g.conn,"Triangles");
+               NetFEM_Scalar(n,g.S11,1,"X Stress (pure)");
+               NetFEM_Scalar(n,g.S22,1,"Y Stress (pure)");
+               NetFEM_Scalar(n,g.S12,1,"Shear Stress (pure)");
+           
+           NetFEM_End(n);
+    }
+  }
+
+  if (myId==0) {
+    double elapsed=CkWallTimer()-totalStart;
+    CkPrintf("Driver finished: average %.6f s/step\n",elapsed/tSteps);
+  }
+}
+
+
+// A PUP function to allow for migration and load balancing of mesh partitions.
+// The PUP function is not needed if no migration or load balancing is desired.
+void pup_myGlobals(pup_er p,myGlobals *g) 
+{
+  FEM_Print("-------- called pup routine -------");
+  pup_int(p,&g->nnodes);
+  pup_int(p,&g->nelems);
+  int nnodes=g->nnodes, nelems=g->nelems;
+  if (pup_isUnpacking(p)) {
+    g->coord=new vector2d[nnodes];
+    g->conn=new connRec[nelems];
+    g->R_net=new vector2d[nnodes]; //Net force
+    g->d=new vector2d[nnodes];//Node displacement
+    g->v=new vector2d[nnodes];//Node velocity
+    g->a=new vector2d[nnodes];
+       g->S11=new double[nelems];
+       g->S22=new double[nelems];
+       g->S12=new double[nelems];
+  }
+  pup_doubles(p,(double *)g->coord,2*nnodes);
+  pup_ints(p,(int *)g->conn,3*nelems);
+  pup_doubles(p,(double *)g->R_net,2*nnodes);
+  pup_doubles(p,(double *)g->d,2*nnodes);
+  pup_doubles(p,(double *)g->v,2*nnodes);
+  pup_doubles(p,(double *)g->a,2*nnodes);
+  pup_doubles(p,(double *)g->S11,nelems);
+  pup_doubles(p,(double *)g->S22,nelems);
+  pup_doubles(p,(double *)g->S12,nelems);
+  if (pup_isDeleting(p)) {
+    delete[] g->coord;
+    delete[] g->conn;
+    delete[] g->R_net;
+    delete[] g->d;
+    delete[] g->v;
+    delete[] g->a;
+       delete[] g->S11;
+       delete[] g->S22;
+       delete[] g->S12;
+  }
+}
+
+
+
+//Update node position, velocity, acceleration based on net force.
+void advanceNodes(const double dt,int nnodes,const vector2d *coord,
+                 vector2d *R_net,vector2d *a,vector2d *v,vector2d *d,bool dampen)
+{
+  const double nodeMass=1.0e-6; //Node mass, kilograms (HACK: hardcoded)
+  const double xm=1.0/nodeMass; //Inverse of node mass
+  const vector2d z(0,0);
+
+  const double shearForce=1.0e-6/(dt*dt*xm);
+
+  bool someNaNs=false;
+  int i;
+  for (i=0;i<nnodes;i++) {
+    vector2d R_n=R_net[i];
+#if NANCHECK
+    if (((R_n.x-R_n.x)!=0)) {
+           CkPrintf("%d (%.4f,%.4f)   ",i,coord[i].x,coord[i].y);
+           someNaNs=true;
+    }
+#endif
+    R_net[i]=z;
+//Apply boundary conditions (HACK: hardcoded!)
+    if (1) {
+       if (coord[i].x<0.00001)
+              continue; //Left edge will NOT move
+       if (coord[i].y>0.02-0.00001)
+              R_n.x+=shearForce; //Top edge pushed hard to right
+    }
+//Update displacement and velocity
+    vector2d aNew=R_n*xm;
+    v[i]+=(dt*0.5)*(aNew+a[i]);
+    d[i]+=dt*v[i]+(dt*dt*0.5)*aNew;
+    a[i]=aNew;   
+    //if (coord[i].y>0.02-0.00001) d[i].y=0.0; //Top edge in horizontal slot
+  }
+  if (dampen)
+    for (i=0;i<nnodes;i++)
+         v[i]*=0.9; //Dampen velocity slightly (prevents eventual blowup)
+  if (someNaNs) {
+         CkPrintf("Nodes all NaN!\n");
+         CkAbort("Node forces NaN!");
+  }
+}
diff --git a/examples/ParFUM/simple2D/pgm.h b/examples/ParFUM/simple2D/pgm.h
new file mode 100644 (file)
index 0000000..fc0ae37
--- /dev/null
@@ -0,0 +1,51 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include "charm++.h"
+#include "ParFUM.h"
+#include "netfem.h"
+#include "vector2d.h"
+
+//One element's connectivity information
+typedef int connRec[3];
+
+// A structure for handling data that may need to be migrated
+struct myGlobals {
+  int nnodes;
+  int nelems;
+  vector2d *coord;
+  connRec *conn;
+
+  vector2d *R_net, *d, *v, *a;
+  
+  double *S11, *S22, *S12;
+};
+
+//Compute forces on constant-strain triangles:
+void CST_NL(const vector2d *coor,const connRec *lm,vector2d *R_net,
+           const vector2d *d,const double *c,
+           int numnp,int numel,
+           double *S11o,double *S22o,double *S12o);
+
+// Prototypes
+void advanceNodes(const double dt,int nnodes,const vector2d *coord,
+                  vector2d *R_net,vector2d *a,vector2d *v,vector2d *d,bool dampen);
+
+void pup_myGlobals(pup_er p,myGlobals *g);
+
+//The material constants c, as computed by fortran mat_const
+// I think the units here are Pascals (N/m^2)
+const double matConst[4]={3.692e9,  1.292e9,  3.692e9,  1.200e9 };
+
+//The timestep, in seconds
+const double dt=1.0e-9;
+
+// A convenient error function
+static void die(const char *str) {
+  CkError("Fatal error: %s\n",str);
+  CkExit();
+}
+
+
+
+#define NANCHECK 1 /*Check for NaNs at each timestep*/
diff --git a/examples/ParFUM/simple2D/vector2d.h b/examples/ParFUM/simple2D/vector2d.h
new file mode 100644 (file)
index 0000000..08cb1e4
--- /dev/null
@@ -0,0 +1,102 @@
+/*
+Orion's Standard Library
+Orion Sky Lawlor, 2/22/2000
+NAME:          vector2d.h
+
+DESCRIPTION:   C++ 2-Dimentional vector library (no templates)
+
+This file provides various utility routines for easily
+manipulating 2-D vectors-- included are arithmetic,
+dot product, magnitude and normalization terms. 
+All routines are provided right in the header file (for inlining).
+
+Converted from vector3d.h.
+
+*/
+
+#ifndef __OSL_VECTOR_2D_H
+#define __OSL_VECTOR_2D_H
+
+#include <math.h>
+
+typedef double real;
+
+//vector2d is a cartesian vector in 2-space-- an x and y.
+class vector2d {
+public:
+       real x,y;
+       vector2d(void) {}//Default consructor
+       //Simple 1-value constructor
+       explicit vector2d(const real init) {x=y=init;}
+       //Simple 1-value constructor
+       explicit vector2d(int init) {x=y=init;}
+       //2-value constructor
+       vector2d(const real Nx,const real Ny) {x=Nx;y=Ny;}
+       //Copy constructor
+       vector2d(const vector2d &copy) {x=copy.x;y=copy.y;}
+       
+       //Cast-to-real * operators (treat vector as array)
+       operator real *() {return &x;}
+       operator const real *() const {return &x;}
+       
+/*Arithmetic operations: these are carefully restricted to just those
+ that make unambiguous sense (to me... now...  ;-)
+Counterexamples: vector*vector makes no sense (use .dot()) because
+real/vector is meaningless (and we'd want a*b/b==a for b!=0), 
+ditto for vector&vector (dot?), vector|vector (projection?), 
+vector^vector (cross?),real+vector, vector+=real, etc.
+*/
+       vector2d &operator=(const vector2d &b) {x=b.x;y=b.y;return *this;}
+       int operator==(const vector2d &b) const {return (x==b.x)&&(y==b.y);}
+       int operator!=(const vector2d &b) const {return (x!=b.x)||(y!=b.y);}
+       vector2d operator+(const vector2d &b) const {return vector2d(x+b.x,y+b.y);}
+       vector2d operator-(const vector2d &b) const {return vector2d(x-b.x,y-b.y);}
+       vector2d operator*(const real scale) const 
+               {return vector2d(x*scale,y*scale);}
+       friend vector2d operator*(const real scale,const vector2d &v)
+               {return vector2d(v.x*scale,v.y*scale);}
+       vector2d operator/(const real &div) const
+               {real scale=1.0/div;return vector2d(x*scale,y*scale);}
+       vector2d operator-(void) const {return vector2d(-x,-y);}
+       void operator+=(const vector2d &b) {x+=b.x;y+=b.y;}
+       void operator-=(const vector2d &b) {x-=b.x;y-=b.y;}
+       void operator*=(const real scale) {x*=scale;y*=scale;}
+       void operator/=(const real div) {real scale=1.0/div;x*=scale;y*=scale;}
+
+//Vector-specific operations
+       //Return the square of the magnitude of this vector
+       real magSqr(void) const {return x*x+y*y;}
+       //Return the magnitude (length) of this vector
+       real mag(void) const {return sqrt(magSqr());}
+       
+       //Return the square of the distance to the vector b
+       real distSqr(const vector2d &b) const 
+               {return (x-b.x)*(x-b.x)+(y-b.y)*(y-b.y);}
+       //Return the distance to the vector b
+       real dist(const vector2d &b) const {return sqrt(distSqr(b));}
+       
+       //Return the dot product of this vector and b
+       real dot(const vector2d &b) const {return x*b.x+y*b.y;}
+       //Return the cosine of the angle between this vector and b
+       real cosAng(const vector2d &b) const {return dot(b)/(mag()*b.mag());}
+       
+       //Return the "direction" (unit vector) of this vector
+       vector2d dir(void) const {return (*this)/mag();}
+
+       //Return the CCW perpendicular vector
+       vector2d perp(void) const {return vector2d(-y,x);}
+
+       //Return this vector scaled by that
+       vector2d &scale(const vector2d &b) {x*=b.x;y*=b.y;return *this;}
+       
+       //Return the largest coordinate in this vector
+       real max(void) {return (x>y)?x:y;}
+       //Make each of this vector's coordinates at least as big
+       // as the given vector's coordinates.
+       void enlarge(const vector2d &by)
+       {if (by.x>x) x=by.x; if (by.y>y) y=by.y;}
+};
+
+#endif //__OSL_VECTOR2D_H
+
+
diff --git a/examples/ParFUM/simple2D/xxx.1.ele b/examples/ParFUM/simple2D/xxx.1.ele
new file mode 100644 (file)
index 0000000..40e33bd
--- /dev/null
@@ -0,0 +1,64 @@
+62  3  1
+   1      16    28    29  101
+   2      34    41    18  101
+   3      40     8    33  101
+   4      10    36    23  101
+   5      21    36    44  101
+   6      30     9    31  101
+   7       3    12    32  101
+   8      32    12    31  101
+   9      23    25    11  101
+  10      22    23    11  101
+  11      26    14    25  101
+  12      39    40    13  101
+  13      30    31    12  101
+  14      38    41    28  101
+  15      14    34    11  101
+  16      26    25     5  101
+  17      37    45    27  101
+  18       8    22    24  101
+  19      27    29    28  101
+  20      35    19    11  101
+  21      42    43    35  101
+  22      17     1    20  101
+  23       5    23    36  101
+  24      11    19    24  101
+  25      10    23    22  101
+  26      25    23     5  101
+  27      11    24    22  101
+  28      26     5    15  101
+  29      11    25    14  101
+  30      27    26    15  101
+  31      28    14    26  101
+  32      17    20    29  101
+  33      27    28    26  101
+  34      38     7    41  101
+  35      17    29    27  101
+  36      16    29    20  101
+  37      12     4    30  101
+  38      44    31     9  101
+  39      22    40    10  101
+  40      39    32    31  101
+  41       3    32    13  101
+  42      13    32    39  101
+  43      41    34    14  101
+  44      35    34    18  101
+  45      34    35    11  101
+  46      35    18    42  101
+  47      39    44    10  101
+  48       5    36    21  101
+  49      27    15    37  101
+  50      28    16    38  101
+  51      44    39    31  101
+  52      10    40    39  101
+  53       8    40    22  101
+  54      40    33    13  101
+  55      28    41    14  101
+  56      18    41     7  101
+  57      43    42     2  101
+  58      21    44     9  101
+  59      43    19    35  101
+  60      10    44    36  101
+  61      45    37     6  101
+  62      45    17    27  101
+# Generated by /home/net/olawlor/bin/triangle -pqVAa xxx.poly
diff --git a/examples/ParFUM/simple2D/xxx.1.node b/examples/ParFUM/simple2D/xxx.1.node
new file mode 100644 (file)
index 0000000..ff48d9b
--- /dev/null
@@ -0,0 +1,47 @@
+45  2  0  1
+   1    0  0    -6
+   2    0.02  0    -2
+   3    0.02  0.02    -2
+   4    0.01  0.02    -2
+   5    0.01  0.01    -1
+   6    0  0.01    -6
+   7    0.01  0    -1
+   8    0.02  0.01    -2
+   9    0.01  0.014999999999999999    -1
+  10    0.015000000000000001  0.012499999999999999    0
+  11    0.015000000000000001  0.006249999999999996    0
+  12    0.014999999999999999  0.02    -2
+  13    0.02  0.014999999999999999    -2
+  14    0.01015625  0.0050000000000000001    0
+  15    0.0050000000000000001  0.01    -1
+  16    0.0050000000000000001  0    -1
+  17    0  0.00390625    -6
+  18    0.014999999999999999  0    -1
+  19    0.02  0.0050000000000000001    -2
+  20    0.001953125  0    -1
+  21    0.01  0.012500000000000001    -1
+  22    0.016562500000000001  0.0093749999999999979    0
+  23    0.013224431818181817  0.009090909090909087    0
+  24    0.02  0.0074999999999999997    -2
+  25    0.011044034090909089  0.0075301846590909091    0
+  26    0.0074518335752407029  0.0073697619682169752    0
+  27    0.0034935834451690091  0.0061378770449895311    0
+  28    0.0065259385819577895  0.0035850672436060298    0
+  29    0.0033651099586732946  0.0029078733795882871    0
+  30    0.01  0.018046875    -1
+  31    0.013476562500000001  0.016523437499999998    0
+  32    0.017035827636718751  0.017035827636718751    0
+  33    0.02  0.012500000000000001    -2
+  34    0.013223286290322579  0.0031249999999999984    0
+  35    0.016875000000000001  0.0031249999999999971    0
+  36    0.012500000000000001  0.011635150331439392    0
+  37    0.0025000000000000001  0.01    -1
+  38    0.0074999999999999997  0    -1
+  39    0.015513188995643435  0.014994450809030037    0
+  40    0.017500000000000002  0.011796874999999998    0
+  41    0.010722244901359637  0.0024798712530825114    0
+  42    0.017500000000000002  0    -1
+  43    0.02  0.0025000000000000001    -2
+  44    0.012426329003210664  0.0141915584926731    0
+  45    0  0.0069531250000000001    -6
+# Generated by /home/net/olawlor/bin/triangle -pqVAa xxx.poly
diff --git a/examples/ParFUM/simple2D/xxx.poly b/examples/ParFUM/simple2D/xxx.poly
new file mode 100644 (file)
index 0000000..2a40e66
--- /dev/null
@@ -0,0 +1,21 @@
+6 2 0 1
+1    0.000    0.000   -6
+2    0.020    0.000   -2
+3    0.020    0.020   -2
+4    0.010    0.020   -2
+5    0.010    0.010   -1
+6    0.000    0.010   -6
+6 1
+1    1    2   -1
+2    2    3   -2
+3    3    4   -2
+4    4    5   -1
+5    5    6   -1
+6    6    1   -6
+0
+1
+1       0.005     0.005       0101    0.000002
+
+
+# use:  triangle -pqVAa xxx.poly
+