Cleaned up code and use new interface FEM_Mesh_data().
authorIsaac Dooley <idooley2@illinois.edu>
Mon, 13 Jun 2005 21:33:49 +0000 (21:33 +0000)
committerIsaac Dooley <idooley2@illinois.edu>
Mon, 13 Jun 2005 21:33:49 +0000 (21:33 +0000)
examples/fem/simple2D/pgm.C
examples/fem/simple2D/pgm.h

index 8d5066f44377f43bedc293ced2e27d2065d22b32..f1e57394a31e42aafdb488ff798e6a35011c7861 100644 (file)
       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, olawlor@acm.org
- */
-#include <stdlib.h>
-#include <stdio.h>
-#include <math.h>
-#include "charm++.h"
-#include "fem.h"
-#include "netfem.h"
-
-#include "vector2d.h"
-#include "pgm.h"
 
-//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;
-
-static void die(const char *str) {
-  CkError("Fatal error: %s\n",str);
-  CkExit();
-}
+ Updated to new FEM interface by
+    Isaac Dooley, 2005
+ */
 
+#include "pgm.h"
 
-#define NANCHECK 1 /*Check for NaNs at each timestep*/
 
 extern "C" void
 init(void)
@@ -93,11 +74,10 @@ init(void)
       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!");
+       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);
   }
@@ -143,107 +123,15 @@ init(void)
  
    
   delete[] ele;
-
   delete[] pts;
 
   CkPrintf("Finished with init (Reading took %.f s)\n",CmiWallTimer()-startTime);
 
 }
 
-//Update node position, velocity, accelleration 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!");
-  }
-}
-
-struct myGlobals {
-  int nnodes;
-  int nelems;
-  vector2d *coord;
-  connRec *conn;
-
-  vector2d *R_net, *d, *v, *a;
-  
-  double *S11, *S22, *S12;
-};
-
-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;
-  }
-}
 
+// A driver() function 
+// driver() is required in all FEM programs
 extern "C" void
 driver(void)
 {
@@ -268,7 +156,7 @@ driver(void)
   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 accelleration
+  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);
 
@@ -297,7 +185,7 @@ driver(void)
        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);//(t%16)==0);
+       advanceNodes(dt,nnodes,g.coord,g.R_net,g.a,g.v,g.d,0);
     }
 
     //Debugging/perf. output
@@ -316,10 +204,10 @@ driver(void)
            }
     }
     /* perform migration-based load balancing */
-    if (t%1024==0)
+    if (0)
       FEM_Migrate();
     
-    if (1) { //Publish data to the net
+    if (t==4608) { //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)");
@@ -342,3 +230,89 @@ driver(void)
 }
 
 
+// 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!");
+  }
+}
index b401793d47ecc288be5da43322848f201560c3de..b8a195181b092bc508fec0f6fa781d1ec8f7bd96 100644 (file)
@@ -1,10 +1,51 @@
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include "charm++.h"
+#include "fem.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*/