Major reorganization and cleaning--separated FEM-specific
authorOrion Lawlor <olawlor@acm.org>
Wed, 22 Jan 2003 20:34:53 +0000 (20:34 +0000)
committerOrion Lawlor <olawlor@acm.org>
Wed, 22 Jan 2003 20:34:53 +0000 (20:34 +0000)
code (fem_*) from rest of code, to show how to write a
serial-or-FEM program.  Also updated for new FEM framework
"mesh" interface.

Still needed:
  -Example problem with cohesive elements
  -NetFEM output (and/or a real output file)

16 files changed:
examples/fem/crack2D/Makefile
examples/fem/crack2D/README
examples/fem/crack2D/cohesive.inp
examples/fem/crack2D/config.C [new file with mode: 0644]
examples/fem/crack2D/crack.h
examples/fem/crack2D/driver.C [deleted file]
examples/fem/crack2D/fem_main.C [new file with mode: 0644]
examples/fem/crack2D/fem_mesh.C [new file with mode: 0644]
examples/fem/crack2D/getmesh.C [deleted file]
examples/fem/crack2D/lst_NL.C
examples/fem/crack2D/lst_coh2.C
examples/fem/crack2D/mesh.C [new file with mode: 0644]
examples/fem/crack2D/node.C [new file with mode: 0644]
examples/fem/crack2D/readfile.C [deleted file]
examples/fem/crack2D/serial_main.C [new file with mode: 0644]
examples/fem/crack2D/volume.C [deleted file]

index 7a3c75e13d67a0ff495674c715f79f3b6ea15e4c..671871c0eacdc470a347499ef8341234ee90006e 100644 (file)
@@ -1,22 +1,20 @@
 CHARMC=../../../../bin/charmc $(OPTS)
+COMMON_OBJS=   config.o mesh.o node.o \
+       lst_NL.o lst_coh2.o
+FEM_OBJS=fem_main.o fem_mesh.o $(COMMON_OBJS)
+SERIAL_OBJS=serial_main.o $(COMMON_OBJS)
 
-all: getmesh pgm
+all: pgm
 
-getmesh: getmesh.C
-       $(CHARMC) -seq -o getmesh getmesh.C -language c++
+# Common files:
+config.o: config.C crack.h
+       $(CHARMC) -c config.C
 
-pgm: driver.o readfile.o volume.o lst_NL.o lst_coh2.o
-       $(CHARMC) -o pgm readfile.o volume.o lst_NL.o lst_coh2.o driver.o \
-         -language fem 
+mesh.o: mesh.C crack.h
+       $(CHARMC) -c mesh.C
 
-readfile.o: readfile.C crack.h
-       $(CHARMC) -c readfile.C
-
-driver.o: driver.C crack.h
-       $(CHARMC) -c driver.C
-
-volume.o: volume.C crack.h
-       $(CHARMC) -c volume.C
+node.o: node.C crack.h
+       $(CHARMC) -c node.C
 
 lst_NL.o: lst_NL.C crack.h
        $(CHARMC) -c lst_NL.C
@@ -24,6 +22,25 @@ lst_NL.o: lst_NL.C crack.h
 lst_coh2.o: lst_coh2.C crack.h
        $(CHARMC) -c lst_coh2.C
 
+
+# FEM parallel program
+fem_main.o: fem_main.C crack.h
+       $(CHARMC) -c fem_main.C
+
+fem_mesh.o: fem_mesh.C crack.h
+       $(CHARMC) -c fem_mesh.C
+
+pgm: $(FEM_OBJS)
+       $(CHARMC) -o pgm $(FEM_OBJS) -language fem 
+
+# Simple serial program
+serial_main.o: serial_main.C crack.h
+       $(CHARMC) -c serial_main.C
+
+serial: $(SERIAL_OBJS)
+       $(CHARMC) -o serial $(SERIAL_OBJS) -seq -language c++
+
+
 clean:
        rm -rf *.o pgm *.sts *.log *.bak *~ SunWS_cache ii_files ti_files
        rm -f conv-host getmesh charmrun
index 0a456e040ce3d282b336ccbd46a9495a41a9b650..348aedbfe05f1986180bac111af3789582043ce8 100644 (file)
@@ -1,39 +1,58 @@
+INTRODUCTION
+
 This is a 2D crack propagation application built on top of the
-FEM framework.  This program is now OBSOLETE, as it uses an older
-version of the FEM framework.
+FEM framework. It is a fully functional analysis application,
+and includes input, partitioning, and processing.  There is no
+output for now, however.
 
 
-In order to compile the crack2d program, you need to just type "make"
-in this directory. This wll make two programs: pgm and getmesh
+BUILD
 
-getmesh is needed to prepare inputs for the program. Where as pgm is the
+In order to compile the crack2d program, you need to just type "make"
+in this directory. This will make one program: pgm, the parallel
 crack2d application.
 
-After compilation, run:
+There is also a serial version of the exact same program, which
+can be made by typing "make serial".
 
-# this will extract the mesh description from the input file crck_bar.inp
-# and will put it in crck_bar.mesh
+The files in this directory include the input files:
+    cohesive.inp: Main configuration file--timestep and materials
+    crck_bar.inp: Mesh input file--lists nodes and elements
 
-./getmesh
+The common I/O and physics files:
+    config.C: Read configuration file
+    mesh.C: Read and set up nodes and elements 
+    node.C: Physics for node timestepping
+    lst_coh2.C: Physics for cohesive elements
+    lst_NL.C: Physics for volumetric elements
+    crack.h: Main header file
 
-# this will convert the mesh in crck_bar.mesh to a graph in crck_bar.graph
-# needed to do partitioning with metis
+The FEM framework version files:
+    fem_main.C: FEM version's main routine
+    fem_mesh.C: Interfaces mesh with FEM framework
 
-../../../../bin/mesh2graph crck_bar.mesh crck_bar.graph
+The serial version files:
+    serial_main.C: Serial version of fem_main.C
 
-# this will partition the graph in crck_bar.graph into 4 partitions
-# and will create individual files meshdata.Pe* for each partition.
 
-../../../../bin/gmap crck_bar.graph 4
+RUN
 
-# this will run the program pgm on 2 processors, by mapping 4 partitions
+# This will run the program pgm on 2 processors, by mapping 4 partitions
 # onto these two processors. This program runs for 2000 iterations, and
 # prints the time taken for each iteration at the end. Migration is
 # performed after every 25 iterations. Note that the number of iterations can
 # be changed in the cohesive.inp file.
 
-./conv-host +p2 pgm +vp4
+./charmrun +p2 pgm +vp4
 
 I have bigger data files available for this program. SO, if you are interested
 in using those, you should contact me. (They are so big, I dont want to
-check them in this directory.
+check them in this directory.)
+
+
+HISTORY
+
+Originally written by Scott Breitenfeld (1999)
+Converted to C and the FEM framework by members of PPL (1999)
+Used to develop FEM framework by Milind Bhandarkar (2000)
+Updated by Orion Lawlor (2003)
index 5bc76b192ebbb209f43da13dca2aafb54b34316a..5f351f03ae913e54f845f7d45f1224674a7ca601 100644 (file)
@@ -1,7 +1,7 @@
 prefix of output files
 crck_bar
 Execution: ntime,nsteps,ntsint_full,ntsint_energy,restart(0=no,1=yes)
-        2000 20 500 15003 0
+        10000 20 500 15003 0
 nplane(0=stress,1=strain,2=axi) ncoh(0=exp,1=lin) lin(0=nonlin,1=linear)
         0 1 0
 number of volumetric material sets : numat_vol
diff --git a/examples/fem/crack2D/config.C b/examples/fem/crack2D/config.C
new file mode 100644 (file)
index 0000000..7e47add
--- /dev/null
@@ -0,0 +1,138 @@
+/**
+ * Read configuration file for crack propagation code.
+ */
+#include <fstream.h>
+#include <stddef.h>
+#include "crack.h"
+
+/// Globally accessible (write-once) configuration data.
+ConfigurationData config;
+
+#define MAXLINE 1024
+// Silly "skip n lines" macro.  Should actually only skip 
+//  comment lines, but this will have to do...
+#define cl(fd, buffer, n) do {\
+                               for(int _i=0; _i<(n); _i++) \
+                                 fd.getline(buffer, MAXLINE);\
+                             } while(0)
+
+/// This routine reads the config data from the config file.
+void readConfig(const char *configFile,const char *meshFile)
+{
+  ifstream fg; // global parameter file
+  ifstream fm; // mesh file, for reading loading data
+  ConfigurationData *cfg=&config; //Always write to the global version
+  char buf[MAXLINE];
+  int i, itmp;
+  double dtmp;
+  
+  fg.open(configFile);
+  fm.open(meshFile);
+  if (!fg || !fm)
+  {
+    crack_abort("Cannot open configuration file for reading.\n");
+    return;
+  }
+  cl(fg,buf,3); // ignore first three lines
+  fg >> cfg->nTime >> cfg->steps >> cfg->tsintFull 
+     >> cfg->tsintEnergy >> cfg->restart;
+  cl(fg, buf, 2);
+  fg >> cfg->nplane >> cfg->ncoh >> cfg->lin;
+  cl(fg, buf, 2);
+  
+  //read volumetric material properties
+  fg >> cfg->numMatVol;
+  cfg->volm = new VolMaterial[cfg->numMatVol];
+  VolMaterial *v;
+  cl(fg, buf, 2);
+  for (i=0; i<cfg->numMatVol; i++)
+  {
+    v = &(cfg->volm[i]);
+    fg >> v->e1 >> v->e2 >> v->g12 >> v->g23;
+    fg >> v->xnu12 >> v->xnu23 >> v->rho;
+    fg >> v->alpha1 >> v->alpha2 >> v->theta;
+    cl(fg, buf, 2);
+  }
+  
+  //Compute the elastic stiffness constants for each material type
+  for (int matNo=0;matNo<cfg->numMatVol;matNo++)
+  {
+    double x, x1, x2, x3;
+    VolMaterial *vm=&(cfg->volm[matNo]);
+    switch (cfg->nplane)
+    {
+      case 1://Orthotropic plane strain
+        double sT,cT,xnu21;
+        cT = cos(vm->theta*1.74532925199e-2);
+        sT = sin(vm->theta*1.74532925199e-2);
+        xnu21 = vm->xnu12*vm->e2/vm->e1;
+        x = 1.0 - vm->xnu23*vm->xnu23 - 
+          2.0*vm->xnu12*xnu21*(1.0 + vm->xnu23);
+        x1 = vm->e1*(1.0-vm->xnu23*vm->xnu23) / x;
+        x2 = xnu21*vm->e1*(1.0+vm->xnu23) / x;
+        x3 = vm->e2*(vm->xnu23+vm->xnu12*xnu21) / x;
+        vm->c[2] = vm->e2*(1.0-vm->xnu12*xnu21) / x;
+        vm->c[0] = x1*cT*cT*cT*cT + 2.0*(x2+2.0*vm->g12)*cT*cT*sT*sT +
+          vm->c[2]*sT*sT*sT*sT;
+        vm->c[1] = x2*cT*cT + x3*sT*sT;
+        vm->c[3] = vm->g12*cT*cT + vm->g23*sT*sT;
+        break;
+      case 0: //Plane stress (isotropic)
+        vm->c[0] = vm->e1 / (1.0 - vm->xnu12*vm->xnu12);
+        vm->c[1] = vm->e1*vm->xnu12 / (1.0 - vm->xnu12*vm->xnu12);
+        vm->c[2] = vm->c[0];
+        vm->c[3] = vm->e1/ (2.0 * (1.0 + vm->xnu12));
+        break;
+      case 2: //Axisymmetric (isotropic)
+        vm->c[0] = vm->e1 * (1.0 - vm->xnu12) / ((1.0 + vm->xnu12)*
+                                                 (1.0 - 2.0*vm->xnu12));
+        vm->c[1] = vm->e1 * vm->xnu12 / ((1.0 + vm->xnu12)*
+                                         (1.0 - 2.0*vm->xnu12));
+        vm->c[2] = vm->e1 / (2.0*(1.0 + vm->xnu12));
+        break;
+      default:
+        crack_abort("Unknown planar analysis type in config file");
+    }
+  }
+
+  //read cohesive material properties
+  fg >> cfg->numMatCoh;
+  cfg->cohm = new CohMaterial[cfg->numMatCoh];
+  CohMaterial *c;
+  cl(fg, buf, 2);
+  for (i=0; i<cfg->numMatCoh; i++)
+  {
+    c = &(cfg->cohm[i]);
+    fg >> c->deltan >> c->deltat >> c->sigmax
+       >> c->taumax >> c->mu;
+    if (cfg->ncoh)
+      fg >> c->Sinit;
+    cl(fg, buf, 2);
+  }
+  
+  //read impact data
+  fg >> cfg->imp >> cfg->voImp >> cfg->massImp >> cfg->radiusImp;
+  cl(fg, buf, 2);
+  fg >> cfg->eImp >> cfg->xnuImp;
+  cl(fg, buf, 2);
+  cfg->voImp = 0; cfg->massImp = 1.0; cfg->radiusImp = 0; 
+  cfg->eImp = 1.0; cfg->xnuImp = 0.3;
+  
+  //read (& ignore) thermal load
+  fg >> dtmp;
+  fg.close();
+  
+  //read proportional ramp-up for boundary conditions
+  fm >> itmp >> itmp >> cfg->delta >> cfg->numProp;
+  cfg->delta /= (double) cfg->steps;
+  cfg->delta2 = cfg->delta*cfg->delta*0.5;
+  cfg->ts_proportion = new int[cfg->numProp];
+  cfg->proportion = new double[cfg->numProp];
+  for (i=0; i< cfg->numProp; i++) {
+    fm >> cfg->ts_proportion[i] >> cfg->proportion[i];
+  }
+  fm.close();
+}
+
+
+
index 99f58bb6c6e0aa52230c7be80b83a61e356c3eae..d31e18c435cc696dfd0d5299f4be3c660bb276ae 100644 (file)
@@ -1,5 +1,9 @@
+/**
+ * Main header file for crack propagation code.
+ */
 #include <string>
 #include <stdlib.h>
+#include <math.h>
 #include "fem.h"
 
 //Coord describes a property with X and Y components
@@ -13,6 +17,8 @@ struct Coord
 // These are 6-vertex triangles (each side has an extra vertex)
 struct Vol
 {
+  int material; //VolMaterial type
+  int conn[6]; // 6 nodes around this element
   double s11l[3], s22l[3], s12l[3];//Stress coefficients
 };
 
@@ -21,23 +27,14 @@ struct Vol
 //  These are 6-vertex rectangles (long sides have an extra vertex)
 struct Coh
 {
+  int material; //CohMaterial type
+  int conn[6]; // 6 nodes around this element
   Coord T[3];//Tractions at each sample point
   double sidel[3];//[0]->length of element; 
                 //[1],[2] give cosine and sine of orientation
   double Sthresh[3];//The threshold, and damage for this edge
 };
 
-// Declaration for class Element (generic elements)
-struct Element
-{
-  int material;          // matlst, matclst, matc
-  union {
-    Vol v;
-    Coh c;
-  };
-};
-
-  
 // Declaration for class Material
 struct Material
 {
@@ -66,20 +63,20 @@ struct CohMaterial : public Material
   
 struct Node
 {
-  Coord Rin;             //Internal force
-  Coord Rco;             //Cohesive traction load
-  Coord xM;              //Mass at this node (xM.x==xM.y, too)
+  Coord Rin;   //Internal force, from volumetric elements
+  Coord Rco;   //Cohesive traction load, from cohesive elements
+  Coord xM;    //Mass at this node (xM.x==xM.y)
   Coord vel;
   Coord accel;
-  Coord disp;
-  Coord pos;
-  Coord r;
-  unsigned char id1, id2;
-  unsigned char isbnd;
+  Coord disp;  //Distance this node has moved
+  Coord pos;   //Undeformed position of this node
+  //FIXME: boundary conditions should be stored in a separate array
+  Coord r; //Boundary condition vector
+  unsigned char isbnd; //Flag: is this a boundary node?
+  unsigned char id1, id2; //Boundary condition flags
 };
 
 //Global constants
-const int    numBoundMax = 1123;   // Maximum number of boundary nodes
 const double   g1          = -0.774596669241483;
 const double   g3          = 0.774596669241483;
 const double   w1          = 0.555555555555555;
@@ -87,25 +84,35 @@ const double   w2          = 0.888888888888888;
 const double   w3          = 0.555555555555555;
 const double   pi          = 3.14159265358979;
 
-struct GlobalData {
-  int myid;
-  int numNP;               //number of nodal points (numnp)
-  int numLST;              //number of LST elements (numlst)
-  int numCLST;             //number of LST cohesive elements (numclst)
-  int numBound;            //number of boundary nodes w/loads
+/// Contains data read from the configuration file that
+/// never changes during the run.
+struct ConfigurationData {
+  /// Simulation control
   int nTime;               //total number of time steps
+  int tsintFull;          //output interval in time steps
+  int tsintEnergy;        //partial output interval
+  int restart;
+  
+  /// Timestep control
   int steps;               //ratio of delta and Courant cond
   double delta;            //timestep
   double delta2;           //delta squared times 2????
-  int *ts_proportion;      //time step for given constant
-  double *proportion;      //load proportionality constant
+  int numProp;             //number of proportionality constants
+  int *ts_proportion;      //time step to apply proportionality constant at
+  double *proportion;      //boundary load proportionality constant
 
+  /// Material formulation
   int lin;                 //elastic formulation, 0=nl, 1=linear
   int ncoh;                //type of cohesive law, 0=exp., 1=linear
   int nplane;              //type of plane analysis, 0=stress, 1=strain
-  int tsintFull;          //output interval in time steps
-  int tsintEnergy;        //partial output interval
-  int restart;
+  
+  /// Material properties
+  int numMatVol;           //number of volumetric materials (numat_vol)
+  int numMatCoh;           //number of cohesive materials (numat_coh)
+  VolMaterial *volm;  // Properties of volumetric materials
+  CohMaterial *cohm;  // Properties of cohesive materials
+  
+  /// "Impact": special boundary condition on just one node
   int imp;                //Is there impact?
   double voImp;         //Velocity of impactor
   double dImp;             //Displacement
@@ -113,26 +120,50 @@ struct GlobalData {
   double vImp;             //Velocity
   double xnuImp,eImp,eTop,radiusImp;   //Impact parameters
   double indent;           //indentation of impactor
-  int nImp;            //node hit
   double fImp;             //contact force
   double massImp;          //mass
+};
 
-  int numMatVol;           //number of volumetric materials (numat_vol)
-  int numMatCoh;           //number of cohesive materials (numat_coh)
-  int numProp;             //number of proportionality constants
-  VolMaterial *volm;
-  CohMaterial *cohm;
+extern ConfigurationData config;
+void readConfig(const char *configFile,const char *meshFile);
 
-  double *itimes;
+/// This structure describes the nodes and elements of a mesh:
+struct MeshData {
+  int nn;               //number of nodal points (numnp)
+  int ne;              //number of LST elements (numlst)
+  int nc;             //number of LST cohesive elements (numclst)
+  int numBound;            //number of boundary nodes w/loads
+  int nImp;            //node hit by impactor
 
-  int nn, ne, npere;
-  int *nnums, *enums, *conn;
   Node *nodes;
-  Element *elements;
-  int scoh, ecoh, svol, evol;
+  Vol *vols;
+  Coh *cohs;
 };
 
-extern void readFile(GlobalData *gd);
-extern void vol_elem(GlobalData *gd);
-extern void lst_NL(GlobalData *gd);
-extern void lst_coh2(GlobalData *gd);
+//Serial mesh routines: in mesh.C
+void readMesh(MeshData *mesh,const char *meshFile);
+void setupMesh(MeshData *mesh); //Initialize newly-read mesh
+void deleteMesh(MeshData *mesh); //Free storage allocated in mesh
+
+//Parallel mesh routines: in fem_mesh.C
+void sendMesh(MeshData *mesh,int fem_mesh);
+void recvMesh(MeshData *mesh,int fem_mesh);
+extern "C" void pupMesh(pup_er p,MeshData *mesh); //For migration
+
+
+//Node physics: in node.C
+struct NodeSlope {
+  int kk; //Index into config.ts_proportion array
+  double prop; //Proportion of boundary conditions to apply
+  double slope; //Slope of prop
+};
+void nodeSetup(NodeSlope *sl);
+void nodeBeginStep(MeshData *mesh);
+void nodeFinishStep(MeshData *mesh, NodeSlope *sl,int tstep);
+
+//Element physics: in lst_NL.C, lst_coh2.C
+extern void lst_NL(MeshData *mesh);
+extern void lst_coh2(MeshData *mesh);
+
+
+void crack_abort(const char *why);
diff --git a/examples/fem/crack2D/driver.C b/examples/fem/crack2D/driver.C
deleted file mode 100644 (file)
index f2f223c..0000000
+++ /dev/null
@@ -1,198 +0,0 @@
-#include <fstream.h>
-#include <stddef.h>
-#include "crack.h"
-
-extern "C" void
-init(void)
-{
-}
-
-static GlobalData *
-mypup(pup_er p, GlobalData *gd)
-{
-  if(pup_isUnpacking(p))
-    gd = new GlobalData;
-  pup_bytes(p, (void*)gd, sizeof(GlobalData));
-  if(pup_isUnpacking(p))
-  {
-    gd->ts_proportion = new int[gd->numProp];
-    gd->proportion = new double[gd->numProp];
-    gd->volm = new VolMaterial[gd->numMatVol];
-    gd->cohm = new CohMaterial[gd->numMatCoh];
-    gd->itimes = new double[gd->nTime];
-    gd->nodes = new Node[gd->nn];
-    gd->elements = new Element[gd->ne];
-  }
-  pup_ints(p, gd->ts_proportion, gd->numProp);
-  pup_doubles(p, gd->proportion, gd->numProp);
-  pup_bytes(p, (void*)gd->volm, gd->numMatVol*sizeof(VolMaterial));
-  pup_bytes(p, (void*) gd->cohm, gd->numMatCoh*sizeof(CohMaterial));
-  pup_doubles(p, gd->itimes, gd->nTime);
-  pup_bytes(p, (void*) gd->nodes, gd->nn*sizeof(Node));
-  pup_bytes(p, (void*) gd->elements, gd->ne*sizeof(Element));
-  if(pup_isPacking(p))
-  {
-    delete[] gd->ts_proportion;
-    delete[] gd->proportion;
-    delete[] gd->volm;
-    delete[] gd->cohm;
-    delete[] gd->itimes;
-    delete[] gd->nodes;
-    delete[] gd->elements;
-    delete gd;
-    gd = 0;
-  }
-  return gd;
-}
-
-static void
-initNodes(GlobalData *gd)
-{
-  int idx;
-  for(idx=0; idx<gd->nn; idx++) {
-    Node *np = &(gd->nodes[idx]);
-    // Zero out node accumulators, update node positions
-    np->Rco.y = np->Rco.x = np->Rin.y =  np->Rin.x = 0.0;
-    np->disp.x += gd->delta2*np->accel.x + gd->delta*np->vel.x;
-    np->disp.y += gd->delta2*np->accel.y + gd->delta*np->vel.y;
-  }
-}
-
-static void 
-updateNodes(GlobalData *gd, double prop, double slope)
-{
-  for(int idx=0; idx<gd->nn; idx++) {
-    Node *np = &(gd->nodes[idx]);
-    if(!np->isbnd) {
-      double aX, aY;
-      aX = (np->Rco.x-np->Rin.x)*np->xM.x;
-      aY = (np->Rco.y-np->Rin.y)*np->xM.y;
-      np->vel.x += (gd->delta*(np->accel.x+aX)*(double) 0.5);
-      np->vel.y += (gd->delta*(np->accel.y+aY)*(double) 0.5);
-      np->accel.x = aX;
-      np->accel.y = aY;
-    } else {
-      double acc;
-      if (!(np->id1)) {
-        np->vel.x = (np->r.x)*prop;
-        np->accel.x = (np->r.x)*slope;
-      } else {
-        acc = (np->r.x*prop+ np->Rco.x - np->Rin.x)*np->xM.x;
-        np->vel.x += (gd->delta*(np->accel.x+acc)*0.5);
-        np->accel.x = acc;
-      }
-      if (!(np->id2)) {
-        np->vel.y = np->r.y*prop;
-        np->accel.y = np->r.y*slope;
-      } else {
-        acc = (np->r.y*prop+ np->Rco.y - np->Rin.y)*np->xM.y;
-        np->vel.y = (np->vel.y + gd->delta*(np->accel.y+acc)*0.5);
-        np->accel.y = acc;
-      }
-    }
-  }
-}
-
-extern "C" double CmiCpuTimer(void);
-static void
-_DELAY_(int microsecs)
-{
-  double upto = CmiCpuTimer() + 1.e-6 * microsecs;
-  while(upto > CmiCpuTimer());
-}
-
-extern "C" void
-driver(int nn, int *nnums, int ne, int *enums, int npere, int *conn)
-{
-  int myid = FEM_My_Partition();
-  int numparts = FEM_Num_Partitions();
-  // CkPrintf("[%d] starting driver\n", myid);
-  GlobalData *gd = new GlobalData;
-  int uidx;
-  uidx = FEM_Register((void*)gd, (FEM_PupFn)mypup);
-  Node *nodes = new Node[nn];
-  Element *elements = new Element[ne];
-  gd->myid = myid;
-  gd->nn = nn;
-  gd->nnums = nnums;
-  gd->ne = ne;
-  gd->enums = enums;
-  gd->npere = npere;
-  gd->conn = conn;
-  gd->nodes = nodes;
-  gd->elements = elements;
-  readFile(gd);
-  gd->itimes = new double[gd->nTime];
-  int i;
-  for(i=0;i<gd->nTime;i++)
-    gd->itimes[i] = 0.0;
-  // CkPrintf("[%d] read file\n", myid);
-  int massfield = FEM_Create_Field(FEM_DOUBLE, 2, offsetof(Node, xM), 
-                                   sizeof(Node));
-  int rfield = FEM_Create_Field(FEM_DOUBLE, 4, offsetof(Node, Rin), 
-                                sizeof(Node));
-  int tfield = FEM_Create_Field(FEM_DOUBLE, 1, 0, 0);
-  FEM_Update_Field(massfield, gd->nodes);
-  int kk = -1;
-  double prop, slope;
-  //double stime, etime;
-  int phase = 0;
-  //stime = CkWallTimer();
-  for(i=0;i<gd->nTime;i++)
-  {
-    gd->itimes[i] = CkWallTimer();
-    // CkPrintf("[%d] iteration %d at %lf secs\n", gd->myid, i, CkTimer());
-    if (gd->ts_proportion[kk+1] == i)
-    {
-      kk++;
-      prop = gd->proportion[kk];
-      slope = (gd->proportion[kk+1]-prop)/gd->delta;
-      slope /= (double) (gd->ts_proportion[kk+1]- gd->ts_proportion[kk]);
-    }
-    else
-    {
-      prop = (double)(i - gd->ts_proportion[kk])*
-                    slope*gd->delta+gd->proportion[kk];
-    }
-    initNodes(gd);
-    lst_NL(gd);
-    lst_coh2(gd);
-    if(myid==79 && i>35)
-    {
-      int biter = (i < 40 ) ? i : 40;
-      _DELAY_((biter-35)*19000);
-    }
-    updateNodes(gd, prop, slope);
-    FEM_Update_Field(rfield, gd->nodes);
-    if(i%25==24)
-    {
-      // if(gd->myid == 0)
-      // {
-        // etime = CkWallTimer();
-        // CkPrintf("Phase=%d\tTime=%lf seconds\n", i, (etime-stime));
-        phase++;
-      // }
-      FEM_Migrate();
-      gd = (GlobalData*) FEM_Get_Userdata(uidx);
-      gd->nnums = FEM_Get_Node_Nums();
-      gd->enums = FEM_Get_Elem_Nums();
-      gd->conn = FEM_Get_Conn();
-      //stime = CkWallTimer();
-    }
-    gd->itimes[i] = CkWallTimer()-(gd->itimes[i]);
-    FEM_Reduce(tfield, &gd->itimes[i], &gd->itimes[i], FEM_SUM);
-  }
-  if(gd->myid==0)
-  {
-    for(i=0;i<gd->nTime;i++)
-    {
-      CkPrintf("Iter=%d\tTime=%lf seconds\n",i,gd->itimes[i]/numparts);
-    }
-  }
-  // CkPrintf("[%d] exiting driver\n", myid);
-}
-
-extern "C" void
-finalize()
-{
-}
diff --git a/examples/fem/crack2D/fem_main.C b/examples/fem/crack2D/fem_main.C
new file mode 100644 (file)
index 0000000..ed25e8c
--- /dev/null
@@ -0,0 +1,147 @@
+/**
+ * Main implementation file for FEM version of crack propagation code.
+ */
+#include <fstream.h>
+#include <stddef.h>
+#include "crack.h"
+#include "charm++.h" // for CkWallTimer, CkPrintf, etc.
+
+void crack_abort(const char *why)
+{
+  CkAbort(why);
+}
+
+extern "C" void
+init(void)
+{
+  //Read and send off mesh data:
+  CkPrintf("Reading mesh...\n");
+  MeshData mesh;
+  readMesh(&mesh,"crck_bar.inp");
+  CkPrintf("Uploading mesh to FEM...\n");
+  sendMesh(&mesh,FEM_Mesh_default_write());
+  deleteMesh(&mesh);
+  CkPrintf("Partitioning mesh...\n");
+}
+
+/// This structure contains all the data held by one chunk of the mesh.
+struct GlobalData {
+  int myid;      // my chunk of the mesh
+  
+  MeshData mesh;
+  
+  double *itimes; // Length of each timestep, for performance tuning
+};
+
+static void mypup(pup_er p, GlobalData *gd)
+{
+  pup_int(p,&gd->myid);
+  pupMesh(p,&gd->mesh);
+  if(pup_isUnpacking(p))
+  {
+    gd->itimes = new double[config.nTime];
+  }
+  pup_doubles(p, gd->itimes, config.nTime);
+  if(pup_isDeleting(p))
+  {
+    delete[] gd->itimes;
+  }
+}
+
+extern "C" double CmiCpuTimer(void);
+static void
+_DELAY_(int microsecs)
+{
+  double upto = CmiCpuTimer() + 1.e-6 * microsecs;
+  while(upto > CmiCpuTimer());
+}
+
+extern "C" void
+driver(void)
+{
+  //Read our configuration data.  We should 
+  //  really use TCHARM_Readonly_globals here, rather than
+  //  re-reading the input file for each chunk.
+  readConfig("cohesive.inp","crck_bar.inp");
+  
+// Set up my chunk's global data:
+  GlobalData gd_storage;
+  GlobalData *gd=&gd_storage;
+  int myid = FEM_My_partition();
+  int numparts = FEM_Num_partitions();
+  FEM_Register((void*)gd, (FEM_PupFn)mypup); //allows migration
+  
+  gd->myid = myid;
+  gd->myid = myid;
+  gd->itimes = new double[config.nTime];
+  for(int i=0;i<config.nTime;i++) gd->itimes[i] = 0.0;
+  
+// Get the FEM mesh:
+  if (myid==0)
+    CkPrintf("Extracting partitioned mesh...\n");
+  int fem_mesh=FEM_Mesh_default_read();
+  recvMesh(&gd->mesh,fem_mesh);
+  
+  // Prepare for communication:
+  int rfield = FEM_Create_field(FEM_DOUBLE, 4, offsetof(Node, Rin), 
+                                sizeof(Node));
+
+// Prepare for the timeloop:
+  if (myid==0)
+    CkPrintf("Running timeloop...\n");
+  NodeSlope sl;
+  nodeSetup(&sl);
+  int t;
+  for(t=0;t<config.nTime;t++) //Timeloop:
+  {
+    double startTime = CkWallTimer(); //Start timing:
+    
+    nodeBeginStep(&gd->mesh);
+    lst_NL(&gd->mesh);
+    lst_coh2(&gd->mesh);
+    
+    //We have the node forces from local elements, but not remote:
+    //  ask FEM to add in the remote node forces.
+    FEM_Update_field(rfield, gd->mesh.nodes);
+    
+    nodeFinishStep(&gd->mesh, &sl, t);
+    
+    if (myid==0) { //For debugging, print the status of my node 0:
+      Node *n=&gd->mesh.nodes[0];
+      int g=-1; //Global number of my node 0
+      FEM_Mesh_get_data(fem_mesh,FEM_NODE,FEM_GLOBALNO, &g, 
+              0,1, FEM_INDEX_0,1);
+      CkPrintf("t=%d  node=%d  d=(%g,%g)  v=(%g,%g)\n",
+           t, g, n->disp.x,n->disp.y, n->vel.x,n->vel.y);
+    }
+    
+    if(0 && myid==79 && t>35) // Add fake load imbalance
+    {
+      int biter = (t < 40 ) ? t : 40;
+      _DELAY_((biter-35)*19000);
+    }
+    
+    if(t%1000==999) // Migrate, for load balance
+    {
+       FEM_Migrate();
+    }
+    
+    // Keep track of how long that timestep took, across the whole machine:
+    gd->itimes[t] = CkWallTimer()-startTime;
+  }
+  
+  if (0) 
+  {
+    // Do a debugging printout of how long each step took:
+    int tfield = FEM_Create_field(FEM_DOUBLE, 1, 0, 0);
+    for(t=0;t<config.nTime;t++)
+    {
+        double thisStep;
+        // Sum across the whole machine:
+        FEM_Reduce(tfield, &gd->itimes[t], &thisStep, FEM_SUM);
+        if(gd->myid==0)
+          CkPrintf("Iteration\t%d\t%.9f\n",t,thisStep/numparts);
+    }
+  }
+}
+
diff --git a/examples/fem/crack2D/fem_mesh.C b/examples/fem/crack2D/fem_mesh.C
new file mode 100644 (file)
index 0000000..a19b821
--- /dev/null
@@ -0,0 +1,120 @@
+/**
+ * Connect mesh to FEM framework for crack propagation code.
+ */
+#include <fstream.h>
+#include <stddef.h>
+#include "crack.h"
+
+
+/// Send/recv this element's material type and element connectivity:
+void sendRecvElement(int fem_mesh,int fem_entity,
+       int matOffset,int connOffset,int connLen,
+       int elSize,int nEl, void *elData)
+{
+  // Element connectivity:
+  IDXL_Layout_t eConn=IDXL_Layout_offset(IDXL_INDEX_0,connLen, //connLen indices
+       connOffset, elSize, 0);
+  FEM_Mesh_data_layout(fem_mesh,fem_entity,FEM_CONN,elData,
+       0,nEl, eConn);
+  IDXL_Layout_destroy(eConn);
+  
+  // Element material type:
+  IDXL_Layout_t eMat=IDXL_Layout_offset(IDXL_INT,1, //1 int
+       matOffset, elSize, 0);
+  FEM_Mesh_data_layout(fem_mesh,fem_entity,FEM_DATA+0,elData,
+       0,nEl, eMat);
+  IDXL_Layout_destroy(eMat);
+}
+
+/// Send/recv this node field
+void sendRecvNode(MeshData *mesh,int fem_mesh,int attrib,
+       int offset,int datatype,int dataLen)
+{
+  IDXL_Layout_t l=IDXL_Layout_offset(datatype,dataLen,
+       offset, sizeof(Node), 0);
+  FEM_Mesh_data_layout(fem_mesh,FEM_NODE,attrib,mesh->nodes,
+       0,mesh->nn, l);
+  IDXL_Layout_destroy(l);
+}
+
+/// Send/recv this mesh's data to/from the FEM framework:
+void sendRecvMesh(MeshData *mesh,int fem_mesh)
+{
+// Nodes:
+  sendRecvNode(mesh,fem_mesh,FEM_DATA+10,
+       offsetof(Node,pos), FEM_DOUBLE, 2);
+  
+  // FIXME: should a list of boundary nodes, instead of
+  // (normally zero!) boundary data for every node
+  sendRecvNode(mesh,fem_mesh,FEM_DATA+20,
+       offsetof(Node,r), FEM_DOUBLE, 2);
+  sendRecvNode(mesh,fem_mesh,FEM_DATA+30,
+       offsetof(Node,isbnd), FEM_BYTE, 3);
+  
+// Volumetric element:
+  sendRecvElement(fem_mesh,FEM_ELEM+0,
+       offsetof(Vol,material),offsetof(Vol,conn),6,
+       sizeof(Vol),mesh->ne,mesh->vols);
+  
+// Cohesive element:
+  sendRecvElement(fem_mesh,FEM_ELEM+1,
+       offsetof(Coh,material),offsetof(Coh,conn),6,
+       sizeof(Coh),mesh->nc,mesh->cohs);
+  
+}
+
+/// Send this mesh out to the FEM framework:
+void sendMesh(MeshData *mesh,int fem_mesh)
+{
+   sendRecvMesh(mesh,fem_mesh);
+}
+
+/// Recv this mesh from the FEM framework:
+void recvMesh(MeshData *mesh,int fem_mesh)
+{
+  // Allocate storage for the mesh
+  mesh->nn=FEM_Mesh_get_length(fem_mesh,FEM_NODE);
+  mesh->ne=FEM_Mesh_get_length(fem_mesh,FEM_ELEM+0);
+  mesh->nc=FEM_Mesh_get_length(fem_mesh,FEM_ELEM+1);
+  mesh->nodes=new Node[mesh->nn];
+  mesh->vols=new Vol[mesh->ne];
+  mesh->cohs=new Coh[mesh->nc];
+   
+  // Grab the data from the FEM framework
+  sendRecvMesh(mesh,fem_mesh);
+  
+  // Initialize the mesh, which fills out the fields we haven't copied
+  setupMesh(mesh);
+  
+  // Sum up element masses (nodes[i].xM) from other processors:
+  int massfield = IDXL_Layout_offset(FEM_DOUBLE, 2, 
+       offsetof(Node, xM), sizeof(Node), 0);
+  int nodeComm = FEM_Comm_shared(fem_mesh,FEM_NODE);
+  IDXL_Comm_sendsum(0,nodeComm,massfield, mesh->nodes);
+  IDXL_Layout_destroy(massfield);
+}
+
+extern "C" void pupMesh(pup_er p,MeshData *mesh)
+{
+  pup_int(p,&mesh->nn);
+  pup_int(p,&mesh->ne);
+  pup_int(p,&mesh->nc);
+  pup_int(p,&mesh->numBound);
+  pup_int(p,&mesh->nImp);
+  if(pup_isUnpacking(p)) {
+    mesh->nodes=new Node[mesh->nn];
+    mesh->vols=new Vol[mesh->ne];
+    mesh->cohs=new Coh[mesh->nc];
+  }
+  //EVIL: instead of packing each node/vol/coh field separately,
+  //  just pack them as flat runs of bytes.
+  pup_bytes(p, (void*)mesh->nodes, mesh->nn*sizeof(Node));
+  pup_bytes(p, (void*)mesh->vols, mesh->ne*sizeof(Vol));
+  pup_bytes(p, (void*)mesh->cohs, mesh->nc*sizeof(Coh));
+  
+  if (pup_isDeleting(p)) {
+    deleteMesh(mesh);
+  }
+}
+
+
diff --git a/examples/fem/crack2D/getmesh.C b/examples/fem/crack2D/getmesh.C
deleted file mode 100644 (file)
index 13134a6..0000000
+++ /dev/null
@@ -1,79 +0,0 @@
-#include <fstream.h>
-#include <stdio.h>
-
-int 
-main(int, char**)
-{
-  ifstream fi;
-  ofstream fo;
-  fi.open("crck_bar.inp");
-  if (!fi)
-  {
-    fprintf(stderr, "Cannot open crck_bar.inp for reading\n");
-    return 1;
-  }
-  fo.open("crck_bar.mesh");
-  if (!fo)
-  {
-    fprintf(stderr, "Cannot open crck_bar.mesh for writing\n");
-    return 1;
-  }
-  int itmp, i, j, k;
-  double dtmp;
-  fi >> itmp >> itmp;
-  fi >> dtmp;
-  int nump;
-  fi >> nump; // number of props
-  for (i=0;i<nump;i++)
-    fi >> itmp >> dtmp;
-  int numnp;
-  fi >> numnp; // number of nodes
-  for (i=0;i<numnp;i++)
-    fi >> itmp >> dtmp >> dtmp;
-  int numb; // num boundary nodes
-  fi >> numb;
-  for(i=0;i<numb;i++)
-    fi >> itmp >> itmp >> itmp >> dtmp >> dtmp;
-  int numclst; // number of cohesive elements
-  fi >> itmp >> numclst >> itmp >> itmp >> itmp;
-  int *cnodes = new int[numclst*6];
-  k = 0;
-  for (i=0; i< numclst; i++)
-  {
-    fi >> itmp;
-    for(j=0;j<6;j++)
-      fi >> cnodes[k++];
-    fi >> itmp;
-  }
-  int numlst; // number of vol elements
-  fi >> itmp >> numlst >> itmp;
-  int *vnodes = new int[numlst*6];
-  k = 0;
-  for (i=0; i< numlst; i++)
-  {
-    fi >> itmp;
-    for(j=0;j<6;j++)
-      fi >> vnodes[k++];
-  }
-  fo << numlst+numclst << ' ' << numnp << ' ' << 6 << endl;
-  k = 0;
-  for(i=0;i<numclst;i++)
-  {
-    for(j=0;j<6;j++)
-      fo << cnodes[k++]-1 << ' ';
-    fo << endl;
-  }
-  k = 0;
-  for(i=0;i<numlst;i++)
-  {
-    for(j=0;j<6;j++)
-      fo << vnodes[k++]-1 << ' ';
-    fo << endl;
-  }
-  cout << "Total Nodes: " << numnp << endl;
-  cout << "Total Cohesive Elements: " << numclst << endl;
-  cout << "Total Volumetric Elements: " << numlst << endl;
-  fi.close();
-  fo.close();
-  return 0;
-}
index 545b4595739a1234c1af65ab66e336a8a0024b88..8e37b2601cb63e72975792fd1e8f996515512964 100644 (file)
@@ -1,8 +1,7 @@
-/*Subroutine: 
-  lst_NL:  converted from Fortran on 9/7/99 by Orion Sky Lawlor
-           Computes each node's Rin.
-*/
-
+/**
+ * Volumetric element physics for crack propagation code.
+ *  converted from Fortran on 9/7/99 by Orion Sky Lawlor
+ */
 #include "crack.h"
 
 //BtoR: sum the "B" spatial derivative array into the
@@ -10,7 +9,7 @@
 
 static void BtoR(const Coord *B/*6-array*/,
     Coord *R/*6-array*/,Node **n/*6-array*/,
-    Element *v,const VolMaterial *vm,
+    Vol *v,const VolMaterial *vm,
     int pkCount,double iaa)
 
 {
@@ -36,9 +35,9 @@ static void BtoR(const Coord *B/*6-array*/,
 //C-----Calculate the pkCount'th Piola-Kirchhoff stress (S)
   
   double s11c,s12c,s22c;
-  s11c=v->v.s11l[pkCount] = E11*vm->c[0] + E22*vm->c[1];
-  s22c=v->v.s22l[pkCount] = E11*vm->c[1] + E22*vm->c[2];
-  s12c=v->v.s12l[pkCount] = E12*vm->c[3];
+  s11c=v->s11l[pkCount] = E11*vm->c[0] + E22*vm->c[1];
+  s22c=v->s22l[pkCount] = E11*vm->c[1] + E22*vm->c[2];
+  s12c=v->s12l[pkCount] = E12*vm->c[3];
   
 //Update R
   for (k=0;k<6;k++) {
@@ -52,15 +51,15 @@ static void BtoR(const Coord *B/*6-array*/,
 
 
 void
-lst_NL(GlobalData *gd)
+lst_NL(MeshData *mesh)
 {
   int idx;
-  for(idx=gd->svol; idx<gd->evol; idx++) {
-    Element *v = &(gd->elements[idx]);
+  for(idx=0; idx<mesh->ne; idx++) {
+    Vol *v = &(mesh->vols[idx]);
     Node *n[6];
     int k;
     for(k=0;k<6;k++)
-      n[k] = &(gd->nodes[gd->conn[idx*6+k]]);
+      n[k] = &(mesh->nodes[v->conn[k]]);
 
     const double c5v3 = 1.66666666666667;
     const double c1v3 = 0.33333333333333;
@@ -76,7 +75,7 @@ lst_NL(GlobalData *gd)
     //Nodes 0, 1, and 2 are the outer corners;
     //Nodes 4, 5, and 6 are on the midpoints of the sides.
       
-    VolMaterial *vm = &(gd->volm[v->material]); //Current material
+    VolMaterial *vm = &(config.volm[v->material]); //Current material
       
     for (k=0;k<6;k++) {
         R[k].x=R[k].y=0.0;
index 752ad9b81d68400502a5ff732a5940516cebcd9f..788026899990786aaf219c944cfb89d8f9fcba71 100644 (file)
@@ -1,40 +1,43 @@
+/**
+ * Cohesive element physics for crack propagation code.
+ */
 #include "crack.h"
 #include <math.h>
 
 
 //Update cohesive element traction itNo given Dn, Dt
 static void 
-updateTraction(double Dn, double Dt,Element *coh,CohMaterial *m,int itNo)
+updateTraction(double Dn, double Dt,Coh *coh,CohMaterial *m,int itNo)
 {
   double Tn, Tt;
-  double c=coh->c.sidel[1],s=coh->c.sidel[2];
+  double c=coh->sidel[1],s=coh->sidel[2];
   double x = 1.0 - sqrt(Dn*Dn+Dt*Dt);
   
   if (x <= 0.0) {
-    coh->c.Sthresh[itNo] = 0.0;
-  } else if (x <= coh->c.Sthresh[itNo])
-    coh->c.Sthresh[itNo] = x;
+    coh->Sthresh[itNo] = 0.0;
+  } else if (x <= coh->Sthresh[itNo])
+    coh->Sthresh[itNo] = x;
   
   if (Dn > 0.0)
-    Tn = coh->c.Sthresh[itNo]/(1.0-coh->c.Sthresh[itNo])*m->sigmax*Dn;
+    Tn = coh->Sthresh[itNo]/(1.0-coh->Sthresh[itNo])*m->sigmax*Dn;
   else
     Tn = m->Sinit/(1.0-m->Sinit)*m->sigmax*Dn;
   
-  Tt = coh->c.Sthresh[itNo]/(1.0-coh->c.Sthresh[itNo])*m->taumax*Dt;
-  coh->c.T[itNo].x = c*Tt - s*Tn;
-  coh->c.T[itNo].y = s*Tt + c*Tn;
-}       
+  Tt = coh->Sthresh[itNo]/(1.0-coh->Sthresh[itNo])*m->taumax*Dt;
+  coh->T[itNo].x = c*Tt - s*Tn;
+  coh->T[itNo].y = s*Tt + c*Tn;
+}
 
 void
-lst_coh2(GlobalData *gd
+lst_coh2(MeshData *mesh
 {
   int idx;
-  for(idx=gd->scoh; idx<gd->ecoh; idx++) {
-    Element *coh = &(gd->elements[idx]);
+  for(idx=0; idx<mesh->nc; idx++) {
+    Coh *coh = &(mesh->cohs[idx]);
     Node *n[6];
     int k;
     for(k=0;k<6;k++)
-      n[k] = &(gd->nodes[gd->conn[idx*6+k]]);
+      n[k] = &(mesh->nodes[coh->conn[k]]);
 
     // local variables:
     double deltn, deltt;  // the char length of current element
@@ -51,13 +54,13 @@ lst_coh2(GlobalData *gd)
     // g1,g3: gauss quadrature points
     // w1,w2,w3: weights; w1=w2, w3=w4
     
-    CohMaterial *m = &(gd->cohm[coh->material]);
+    CohMaterial *m = &(config.cohm[coh->material]);
       
     deltn = (double)1.0 / m->deltan;
     deltt = (double)1.0 / m->deltat;
-    length = coh->c.sidel[0] * (double)0.5;
-    c = coh->c.sidel[1];
-    s = coh->c.sidel[2];
+    length = coh->sidel[0] * (double)0.5;
+    c = coh->sidel[1];
+    s = coh->sidel[2];
     
     Dx1 = n[3]->disp.x - n[0]->disp.x;
     Dy1 = n[3]->disp.y - n[0]->disp.y;
@@ -87,12 +90,12 @@ lst_coh2(GlobalData *gd)
   
     // cohesive forces
     x = length*w1*g1*0.5;
-    Rx1 = (coh->c.T[0].x*(g1-1.0)+coh->c.T[2].x*(g1+1.0))*x;
-    Ry1 = (coh->c.T[0].y*(g1-1.0)+coh->c.T[2].y*(g1+1.0))*x;
-    Rx2 = (coh->c.T[0].x*(g1+1.0)+coh->c.T[2].x*(g1-1.0))*x;
-    Ry2 = (coh->c.T[0].y*(g1+1.0)+coh->c.T[2].y*(g1-1.0))*x;
-    Rx3 = ((coh->c.T[0].x+coh->c.T[2].x)*w1*(1.0-g1*g1)+coh->c.T[1].x*w2)*length;
-    Ry3 = ((coh->c.T[0].y+coh->c.T[2].y)*w1*(1.0-g1*g1)+coh->c.T[1].y*w2)*length;
+    Rx1 = (coh->T[0].x*(g1-1.0)+coh->T[2].x*(g1+1.0))*x;
+    Ry1 = (coh->T[0].y*(g1-1.0)+coh->T[2].y*(g1+1.0))*x;
+    Rx2 = (coh->T[0].x*(g1+1.0)+coh->T[2].x*(g1-1.0))*x;
+    Ry2 = (coh->T[0].y*(g1+1.0)+coh->T[2].y*(g1-1.0))*x;
+    Rx3 = ((coh->T[0].x+coh->T[2].x)*w1*(1.0-g1*g1)+coh->T[1].x*w2)*length;
+    Ry3 = ((coh->T[0].y+coh->T[2].y)*w1*(1.0-g1*g1)+coh->T[1].y*w2)*length;
     n[0]->Rco.x += Rx1;
     n[0]->Rco.y += Ry1;
     n[3]->Rco.x -= Rx1;
diff --git a/examples/fem/crack2D/mesh.C b/examples/fem/crack2D/mesh.C
new file mode 100644 (file)
index 0000000..62e8d61
--- /dev/null
@@ -0,0 +1,166 @@
+/**
+ * Read and maintain basic mesh for crack propagation code.
+ */
+#include <fstream.h>
+#include <stddef.h>
+#include "crack.h"
+
+#define MAXLINE 1024
+// Silly "skip n lines" macro.  Should actually only skip 
+//  comment lines, but this will have to do...
+#define cl(fd, buffer, n) do {\
+                               for(int _i=0; _i<(n); _i++) \
+                                 fd.getline(buffer, MAXLINE);\
+                             } while(0)
+
+/// Read this mesh from this file:
+void
+readMesh(MeshData *mesh, const char *meshFile)
+{
+  ifstream fm; // mesh file
+  char buf[MAXLINE];
+  int i, itmp;
+  double dtmp;
+  
+  fm.open(meshFile);
+  if (!fm)
+  {
+    crack_abort("Cannot open mesh file for reading.\n");
+    return;
+  }
+  
+  //Skip ramp-up for boundary conditions (already read by readConfig)
+  int numProp;
+  fm >> itmp >> itmp >> dtmp >> numProp;
+  for (i=0; i< numProp; i++) {
+    fm >> itmp >> dtmp;
+  }
+  
+  //read nodal co-ordinates
+  fm >> mesh->nn;
+  cl(fm,buf,1);
+  mesh->nodes=new Node[mesh->nn];
+  for(i=0; i<mesh->nn; i++)
+  {
+    Node *np = &(mesh->nodes[i]);
+    fm >> itmp >> np->pos.x >> np->pos.y;
+    np->isbnd = 0;
+  }
+  
+  //read nodal boundary conditions
+  fm >> mesh->numBound;
+  for (i=0; i<mesh->numBound; i++)
+  {
+    int j; // Boundary condition i applies to node j
+    fm >> j;
+    j--; //Fortran to C indexing
+    Node *np = &(mesh->nodes[j]);
+    np->isbnd = 1;
+    fm >> np->id1 >> np->id2 >> np->r.x >> np->r.y;
+  }
+  
+  //read cohesive elements
+  fm >> itmp >> mesh->nc >> itmp >> itmp >> itmp;
+  cl(fm,buf,1);
+  mesh->cohs=new Coh[mesh->nc];
+  for(i=0; i<mesh->nc; i++)
+  {
+    Coh *coh = &(mesh->cohs[i]);
+    fm >> coh->material; coh->material--;
+    int k;
+    for(k=0;k<6;k++) {
+      fm >> coh->conn[k]; coh->conn[k]--;
+    }
+    fm >> itmp; // What *is* this?
+  }
+  
+  // Read volumetric elements
+  fm >> itmp >> mesh->ne >> itmp;
+  cl(fm,buf,1);
+  mesh->vols=new Vol[mesh->ne];
+  for(i=0;i<mesh->ne;i++)
+  {
+    Vol *vol = &(mesh->vols[i]);
+    fm >> vol->material; vol->material--;
+    int k;
+    for(k=0;k<6;k++) {
+      fm >> vol->conn[k]; vol->conn[k]--;
+    }
+  }
+  fm.close();
+}
+
+void setupMesh(MeshData *mesh) {
+  int i;
+  //Zero out nodes
+  for(i=0; i<mesh->nn; i++)
+  {
+    Node *np = &(mesh->nodes[i]);
+    // np->pos is already set
+    np->xM.x = np->xM.y = 0;
+    np->disp.x = np->disp.y = 0.0;
+    np->vel.x = np->vel.y = 0.0;
+    np->accel.x = np->accel.y = 0.0;
+  }
+  
+  // Init cohesive elements, determine the length and angle
+  for(i=0; i<mesh->nc; i++)
+  {
+    Coh *coh = &(mesh->cohs[i]);
+    // coh->material and coh->conn[k] are already set
+    
+    Node *np1 = &(mesh->nodes[coh->conn[1]]);
+    Node *np2 = &(mesh->nodes[coh->conn[0]]);
+    coh->Sthresh[2] = coh->Sthresh[1] =
+      coh->Sthresh[0] = config.cohm[coh->material].Sinit;
+    double x = np1->pos.x - np2->pos.x;
+    double y = np1->pos.y - np2->pos.y;
+    coh->sidel[0] = sqrt(x*x+y*y);
+    coh->sidel[1] = x/coh->sidel[0];
+    coh->sidel[2] = y/coh->sidel[0];
+  }
+  
+  // Set up volumetric elements
+  for(i=0;i<mesh->ne;i++)
+  {
+    Vol *vol = &(mesh->vols[i]);
+    // vol->material and vol->conn[k] are already set
+    for(int k=0;k<3;k++)
+    {
+      vol->s11l[k] = 0.0;
+      vol->s22l[k] = 0.0;
+      vol->s12l[k] = 0.0;
+    }
+  }
+  
+  // Compute the mass of local elements:
+  for (i=0;i<mesh->ne;i++)
+  {
+    Vol *v=&(mesh->vols[i]);
+    Node *n[6];              //Pointers to each of the triangle's nodes
+    int k;                  //Loop index
+    for (k=0;k<6;k++)
+      n[k]=&(mesh->nodes[v->conn[k]]);
+    //Compute the mass of this element
+    double area=((n[1]->pos.x-n[0]->pos.x)*(n[2]->pos.y-n[0]->pos.y)-
+                 (n[2]->pos.x-n[0]->pos.x)*(n[1]->pos.y-n[0]->pos.y));
+    double mass=config.volm[v->material].rho*area/114.0;
+    //Divide the element's mass among the element's nodes
+    for (k=0;k<3;k++) {
+      n[k]->xM.x+=mass*3.0;
+      n[k]->xM.y+=mass*3.0;
+    }
+    for (k=3;k<6;k++) {
+      n[k]->xM.x+=mass*16.0;
+      n[k]->xM.y+=mass*16.0;
+    }
+  }  
+}
+
+void deleteMesh(MeshData *mesh)
+{
+  delete[] mesh->nodes; mesh->nodes=NULL;
+  delete[] mesh->vols;  mesh->vols=NULL;
+  delete[] mesh->cohs;  mesh->cohs=NULL;
+}
+
diff --git a/examples/fem/crack2D/node.C b/examples/fem/crack2D/node.C
new file mode 100644 (file)
index 0000000..03f724c
--- /dev/null
@@ -0,0 +1,75 @@
+/**
+ * Node manipulation for crack propagation code.
+ */
+#include <fstream.h>
+#include <stddef.h>
+#include "crack.h"
+
+void nodeSetup(NodeSlope *sl) {
+  sl->kk=-1;
+  sl->slope=sl->prop=0;
+}
+
+void nodeBeginStep(MeshData *mesh)
+{
+  for(int idx=0; idx<mesh->nn; idx++) {
+    Node *np = &(mesh->nodes[idx]);
+    // Zero out node accumulators, update node positions
+    np->Rco.y = np->Rco.x = np->Rin.y =  np->Rin.x = 0.0;
+    np->disp.x += config.delta2*np->accel.x + config.delta*np->vel.x;
+    np->disp.y += config.delta2*np->accel.y + config.delta*np->vel.y;
+  }
+}
+
+void nodeFinishStep(MeshData *mesh, NodeSlope *sl,int tstep)
+{
+  // Slowly ramp up boundary conditions:
+  if (config.ts_proportion[sl->kk+1] == tstep)
+  {
+      sl->kk++;
+      sl->prop = config.proportion[sl->kk];
+      sl->slope = (config.proportion[sl->kk+1]-sl->prop)/config.delta;
+      sl->slope /= (double) (config.ts_proportion[sl->kk+1]- config.ts_proportion[sl->kk]);
+  }
+  else
+  {
+      sl->prop = (double)(tstep - config.ts_proportion[sl->kk])*
+                    sl->slope*config.delta+config.proportion[sl->kk];
+  }
+  double slope=sl->slope;
+  double prop=sl->prop;
+
+  // Update each node:
+  for(int idx=0; idx<mesh->nn; idx++) {
+    Node *np = &(mesh->nodes[idx]);
+    if(!np->isbnd) { //Apply normal timestep to this node:
+      double aX, aY;
+      aX = (np->Rco.x-np->Rin.x)*np->xM.x;
+      aY = (np->Rco.y-np->Rin.y)*np->xM.y;
+      np->vel.x += (config.delta*(np->accel.x+aX)*(double) 0.5);
+      np->vel.y += (config.delta*(np->accel.y+aY)*(double) 0.5);
+      np->accel.x = aX;
+      np->accel.y = aY;
+    } else { //Apply boundary conditions to this node:
+      double acc;
+      if (!(np->id1)) {
+        np->vel.x = (np->r.x)*prop;
+        np->accel.x = (np->r.x)*slope;
+      } else {
+        acc = (np->r.x*prop+ np->Rco.x - np->Rin.x)*np->xM.x;
+        np->vel.x += (config.delta*(np->accel.x+acc)*0.5);
+        np->accel.x = acc;
+      }
+      if (!(np->id2)) {
+        np->vel.y = np->r.y*prop;
+        np->accel.y = np->r.y*slope;
+      } else {
+        acc = (np->r.y*prop+ np->Rco.y - np->Rin.y)*np->xM.y;
+        np->vel.y = (np->vel.y + config.delta*(np->accel.y+acc)*0.5);
+        np->accel.y = acc;
+      }
+    }
+  }
+}
+
+
diff --git a/examples/fem/crack2D/readfile.C b/examples/fem/crack2D/readfile.C
deleted file mode 100644 (file)
index 60251bc..0000000
+++ /dev/null
@@ -1,169 +0,0 @@
-#include "crack.h"
-#include <fstream.h>
-#include <math.h>
-
-#define MAXLINE 1024
-#define cl(fd, buffer, n) do {\
-                               for(int _i=0; _i<(n); _i++) \
-                                 fd.getline(buffer, MAXLINE);\
-                             } while(0)
-void
-readFile(GlobalData *gd)
-{
-  ifstream fg; // global parameter file
-  ifstream fm; // mesh file
-  char buf[MAXLINE];
-  int i, itmp;
-  double dtmp;
-
-  fg.open("cohesive.inp");
-  fm.open("crck_bar.inp");
-  if (!fg || !fm)
-  {
-    CkAbort("Cannot open input files for reading.\n");
-    return;
-  }
-  cl(fg,buf,3); // ignore first three lines
-  fg >> gd->nTime >> gd->steps >> gd->tsintFull 
-     >> gd->tsintEnergy >> gd->restart;
-  cl(fg, buf, 2);
-  fg >> gd->nplane >> gd->ncoh >> gd->lin;
-  cl(fg, buf, 2);
-  //read volumetric material properties
-  fg >> gd->numMatVol;
-  gd->volm = new VolMaterial[gd->numMatVol];
-  VolMaterial *v;
-  cl(fg, buf, 2);
-  for (i=0; i<gd->numMatVol; i++)
-  {
-    v = &(gd->volm[i]);
-    fg >> v->e1 >> v->e2 >> v->g12 >> v->g23;
-    fg >> v->xnu12 >> v->xnu23 >> v->rho;
-    fg >> v->alpha1 >> v->alpha2 >> v->theta;
-    cl(fg, buf, 2);
-  }
-  //read cohesive material properties
-  fg >> gd->numMatCoh;
-  gd->cohm = new CohMaterial[gd->numMatCoh];
-  CohMaterial *c;
-  cl(fg, buf, 2);
-  for (i=0; i<gd->numMatCoh; i++)
-  {
-    c = &(gd->cohm[i]);
-    fg >> c->deltan >> c->deltat >> c->sigmax
-       >> c->taumax >> c->mu;
-    if (gd->ncoh)
-      fg >> c->Sinit;
-    cl(fg, buf, 2);
-  }
-  //read impact data
-  fg >> gd->imp >> gd->voImp >> gd->massImp >> gd->radiusImp;
-  cl(fg, buf, 2);
-  fg >> gd->eImp >> gd->xnuImp;
-  cl(fg, buf, 2);
-  gd->voImp = 0; gd->massImp = 1.0; gd->radiusImp = 0; 
-  gd->eImp = 1.0; gd->xnuImp = 0.3;
-  //read thermal load
-  fg >> dtmp;
-  fg.close();
-  fm >> itmp >> itmp >> gd->delta >> gd->numProp;
-  gd->delta /= (double) gd->steps;
-  gd->delta2 = gd->delta*gd->delta*0.5;
-  gd->ts_proportion = new int[gd->numProp];
-  gd->proportion = new double[gd->numProp];
-  for (i=0; i< gd->numProp; i++) {
-    fm >> gd->ts_proportion[i] >> gd->proportion[i];
-  }
-  //read nodal co-ordinates
-  fm >> gd->numNP;
-  cl(fm,buf,1);
-  int curline = 0;
-  for(i=0; i<gd->nn; i++)
-  {
-    cl(fm, buf, gd->nnums[i]-curline);
-    curline = gd->nnums[i];
-    Node *np = &(gd->nodes[i]);
-    fm >> itmp >> np->pos.x >> np->pos.y;
-    if(itmp != gd->nnums[i]+1)
-    {
-      CkError("[%d] Expected to read node %d got %d\n", gd->myid, 
-              gd->nnums[i]+1, itmp);
-      CkAbort("");
-    }
-    np->xM.x = np->xM.y = 0;
-    np->disp.x = np->disp.y = 0.0;
-    np->vel.x = np->vel.y = 0.0;
-    np->accel.x = np->accel.y = 0.0;
-  }
-  cl(fm, buf, gd->numNP-curline);
-  //read nodal boundary conditions
-  fm >> gd->numBound;
-  for (i=0; i<gd->numBound; i++)
-  {
-    int j, id0;
-    fm >> id0;
-    id0--;
-    for (j=0; j<gd->nn && id0!=gd->nnums[j]; j++);
-    if(j==gd->nn)
-    {
-      cl(fm, buf, 1);
-      continue;
-    }
-    Node *np = &(gd->nodes[j]);
-    np->isbnd = 1;
-    fm >> np->id1 >> np->id2 >> np->r.x >> np->r.y;
-  }
-  //read cohesive element data, determine the length and angle
-  fm >> itmp >> gd->numCLST >> itmp >> itmp >> itmp;
-  cl(fm,buf,1);
-  curline = 0;
-  for(i=0;i<gd->ne && gd->enums[i]<gd->numCLST;i++);
-  gd->scoh = 0;
-  gd->ecoh = i;
-  gd->svol = i;
-  gd->evol = gd->ne;
-  for(i=0; i<gd->ne && gd->enums[i]<gd->numCLST; i++)
-  {
-    cl(fm, buf, gd->enums[i]-curline);
-    curline = gd->enums[i];
-    Element *ep = &(gd->elements[i]);
-    fm >> ep->material; ep->material--;
-    int k;
-    for(k=0;k<6;k++) {
-      fm >> itmp;
-    }
-    fm >> itmp;
-    Node *np1 = &(gd->nodes[gd->conn[6*i+1]]);
-    Node *np2 = &(gd->nodes[gd->conn[6*i]]);
-    Coh *coh = &(ep->c);
-    coh->Sthresh[2] = coh->Sthresh[1] =
-      coh->Sthresh[0] = gd->cohm[ep->material].Sinit;
-    double x = np1->pos.x - np2->pos.x;
-    double y = np1->pos.y - np2->pos.y;
-    coh->sidel[0] = sqrt(x*x+y*y);
-    coh->sidel[1] = x/coh->sidel[0];
-    coh->sidel[2] = y/coh->sidel[0];
-  }
-  cl(fm, buf, gd->numCLST-curline);
-  fm >> itmp >> gd->numLST >> itmp;
-  cl(fm,buf,1);
-  curline = gd->numCLST;
-  for(gd->svol;i<gd->ne;i++)
-  {
-    cl(fm, buf, gd->enums[i]-curline);
-    curline = gd->enums[i];
-    Element *ep = &(gd->elements[i]);
-    fm >> ep->material; ep->material--;
-    int k;
-    for(k=0;k<6;k++) fm >> itmp;
-    for(k=0;k<3;k++)
-    {
-      ep->v.s11l[k] = 0.0;
-      ep->v.s22l[k] = 0.0;
-      ep->v.s12l[k] = 0.0;
-    }
-  }
-  cl(fm, buf, gd->numCLST-curline);
-  fm.close();
-  vol_elem(gd);
-}
diff --git a/examples/fem/crack2D/serial_main.C b/examples/fem/crack2D/serial_main.C
new file mode 100644 (file)
index 0000000..a9123bd
--- /dev/null
@@ -0,0 +1,47 @@
+/**
+ * Main implementation file for serial version of
+ * crack propagation code.
+ */
+#include <stdio.h>
+#include <stddef.h>
+#include "crack.h"
+
+void crack_abort(const char *why)
+{
+  fprintf(stderr,"Fatal error> %s\n",why);
+  abort();
+}
+
+int main(int argc,char *argv[])
+{
+  //Read mesh data:
+  printf("Reading mesh...\n");
+  MeshData mesh;
+  readMesh(&mesh,"crck_bar.inp");
+  
+  //Read our configuration data.
+  readConfig("cohesive.inp","crck_bar.inp");
+  
+  //Set up the mesh
+  setupMesh(&mesh);
+  
+// Prepare for the timeloop:
+  NodeSlope sl;
+  nodeSetup(&sl);
+  int t;
+  for(t=0;t<config.nTime;t++) //Timeloop:
+  {
+    nodeBeginStep(&mesh);
+    lst_NL(&mesh);
+    lst_coh2(&mesh);
+    
+    nodeFinishStep(&mesh, &sl, t);
+    
+    //For debugging, print the status of node 0:
+    Node *n=&mesh.nodes[0];
+    int g=0; //In a serial program, my node 0 is the true node 0.
+    printf("t=%d  node=%d  d=(%g,%g)  v=(%g,%g)\n",
+           t, g, n->disp.x,n->disp.y, n->vel.x,n->vel.y);
+  }
+}
+
diff --git a/examples/fem/crack2D/volume.C b/examples/fem/crack2D/volume.C
deleted file mode 100644 (file)
index 69df905..0000000
+++ /dev/null
@@ -1,72 +0,0 @@
-#include "crack.h"
-#include <math.h>
-
-void 
-vol_elem(GlobalData *gd)
-{
-  double x, x1, x2, x3;
-
-  //Compute the elastic stiffness constants for each material type
-  int matNo;
-  for (matNo=0;matNo<gd->numMatVol;matNo++)
-  {
-    VolMaterial *vm=&(gd->volm[matNo]);
-    switch (gd->nplane)
-    {
-      case 1://Orthotropic plane strain
-        double sT,cT,xnu21;
-        cT = cos(vm->theta*1.74532925199e-2);
-        sT = sin(vm->theta*1.74532925199e-2);
-        xnu21 = vm->xnu12*vm->e2/vm->e1;
-        x = 1.0 - vm->xnu23*vm->xnu23 - 
-          2.0*vm->xnu12*xnu21*(1.0 + vm->xnu23);
-        x1 = vm->e1*(1.0-vm->xnu23*vm->xnu23) / x;
-        x2 = xnu21*vm->e1*(1.0+vm->xnu23) / x;
-        x3 = vm->e2*(vm->xnu23+vm->xnu12*xnu21) / x;
-        vm->c[2] = vm->e2*(1.0-vm->xnu12*xnu21) / x;
-        vm->c[0] = x1*cT*cT*cT*cT + 2.0*(x2+2.0*vm->g12)*cT*cT*sT*sT +
-          vm->c[2]*sT*sT*sT*sT;
-        vm->c[1] = x2*cT*cT + x3*sT*sT;
-        vm->c[3] = vm->g12*cT*cT + vm->g23*sT*sT;
-        break;
-      case 0: //Plane stress (isotropic)
-        vm->c[0] = vm->e1 / (1.0 - vm->xnu12*vm->xnu12);
-        vm->c[1] = vm->e1*vm->xnu12 / (1.0 - vm->xnu12*vm->xnu12);
-        vm->c[2] = vm->c[0];
-        vm->c[3] = vm->e1/ (2.0 * (1.0 + vm->xnu12));
-        break;
-      case 2: //Axisymmetric (isotropic)
-        vm->c[0] = vm->e1 * (1.0 - vm->xnu12) / ((1.0 + vm->xnu12)*
-                                                 (1.0 - 2.0*vm->xnu12));
-        vm->c[1] = vm->e1 * vm->xnu12 / ((1.0 + vm->xnu12)*
-                                         (1.0 - 2.0*vm->xnu12));
-        vm->c[2] = vm->e1 / (2.0*(1.0 + vm->xnu12));
-        break;
-      default:
-        CkAbort("Unknown planar analysis type passed to vol_elem!\n");
-    }
-  }
-
-  //Update the node-by-node mass of each element
-  int volNo;
-  for (volNo=gd->svol;volNo<gd->ne;volNo++)
-  {
-    Node *n[6];              //Pointers to each of the triangle's nodes
-    int k;                  //Loop index
-    for (k=0;k<6;k++)
-      n[k]=&(gd->nodes[gd->conn[volNo*6+k]]);
-    //Compute the mass of this element
-    double area=((n[1]->pos.x-n[0]->pos.x)*(n[2]->pos.y-n[0]->pos.y)-
-                 (n[2]->pos.x-n[0]->pos.x)*(n[1]->pos.y-n[0]->pos.y));
-    double mass=gd->volm[gd->elements[volNo].material].rho*area/114.0;
-    //Divide the mass among the element's nodes
-    for (k=0;k<3;k++) {
-      n[k]->xM.x+=mass*3.0;
-      n[k]->xM.y+=mass*3.0;
-    }
-    for (k=3;k<6;k++) {
-      n[k]->xM.x+=mass*16.0;
-      n[k]->xM.y+=mass*16.0;
-    }
-  }  
-}