retreiving the prev version
authorNilesh Choudhury <nchoudh2@uiuc.edu>
Sun, 25 Sep 2005 20:08:36 +0000 (20:08 +0000)
committerNilesh Choudhury <nchoudh2@uiuc.edu>
Sun, 25 Sep 2005 20:08:36 +0000 (20:08 +0000)
examples/fem/simple2D/pgm.C
examples/fem/simple2D/pgm.h

index 73bac73089b56cdb2d185a88829f7fedf5dec764..2a78f1b1dcc218d067023ff742f4274804073edb 100644 (file)
-#include "pgm.h"
-#include "mpi.h"
-
-#include "fem_adapt_new.h"
-#include "fem_mesh_modify.h"
-
-#define SUMMARY_ON
-
-extern void _registerFEMMeshModify(void);
-
-double getArea(double *n1_coord, double *n2_coord, double *n3_coord) {
-  double area=0.0;
-  double aLen, bLen, cLen, sLen, d, ds_sum;
-
-  ds_sum = 0.0;
-  for (int i=0; i<2; i++) {
-    d = n1_coord[i] - n2_coord[i];
-    ds_sum += d*d;
-  }
-  aLen = sqrt(ds_sum);
-  ds_sum = 0.0;
-  for (int i=0; i<2; i++) {
-    d = n2_coord[i] - n3_coord[i];
-    ds_sum += d*d;
-  }
-  bLen = sqrt(ds_sum);
-  ds_sum = 0.0;
-  for (int i=0; i<2; i++) {
-    d = n3_coord[i] - n1_coord[i];
-    ds_sum += d*d;
-  }
-  cLen = sqrt(ds_sum);
-  sLen = (aLen+bLen+cLen)/2;
-  if(sLen-aLen < 0) return 0.0;
-  else if(sLen-bLen < 0) return 0.0;
-  else if(sLen-cLen < 0) return 0.0; //area too small to note
-  return (sqrt(sLen*(sLen-aLen)*(sLen-bLen)*(sLen-cLen)));
-}
-
-void publish_data_netfem(int i,  myGlobals g) {
-  MPI_Barrier(MPI_COMM_WORLD);
-  if (1) { //Publish data to the net
-    int mesh=FEM_Mesh_default_read(); // Tell framework we are reading data from the mesh
-    int rank = FEM_My_partition();
-    g.nnodes=FEM_Mesh_get_length(mesh,FEM_NODE); // Get number of nodes
-    g.coord=new vector2d[g.nnodes];
-    int count = 0;
-    vector2d *coord = new vector2d[g.nnodes];
-    int *maptovalid = new int[g.nnodes];
-    double *nodeid = new double[g.nnodes];
-    // Get node positions
-    FEM_Mesh_data(mesh, FEM_NODE, FEM_DATA+0, (double*)g.coord, 0, g.nnodes, FEM_DOUBLE, 2);
-    for(int j=0; j<g.nnodes; j++) {
-      if(FEM_is_valid(mesh,FEM_NODE+0,j)) {
-       coord[count].x = g.coord[j].x;
-       coord[count].y = g.coord[j].y;
-       maptovalid[j] = count;
-       nodeid[count] = j;
-       count++;
-      }
-    }
-    NetFEM n=NetFEM_Begin(rank,i,2,NetFEM_WRITE);
-    NetFEM_Nodes(n,count,(double *)coord,"Position (m)");
-    //NetFEM_Nodes(n,g.nnodes,(double *)g.coord,"Position (m)");
-    NetFEM_Scalar(n,nodeid,1,"Node ID");
-
-    // Get element data
-    g.nelems=FEM_Mesh_get_length(mesh,FEM_ELEM+0); // Get number of elements
-    g.conn=new connRec[g.nelems];
-    connRec *conn = new connRec[g.nelems];
-    double *elid = new double[g.nelems];
-    count = 0;
-    // Get connectivity for elements
-    FEM_Mesh_data(mesh, FEM_ELEM+0, FEM_CONN, (int *)g.conn, 0, g.nelems, FEM_INDEX_0, 3);
-    double totalArea = 0.0;
-    for(int j=0; j<g.nelems; j++) {
-      if(FEM_is_valid(mesh,FEM_ELEM+0,j)) {
-       conn[count][0] = maptovalid[g.conn[j][0]];
-       conn[count][1] = maptovalid[g.conn[j][1]];
-       conn[count][2] = maptovalid[g.conn[j][2]];
-       elid[count] = j;
-       totalArea += getArea(coord[conn[count][0]],coord[conn[count][1]],coord[conn[count][2]]);
-       if(totalArea != totalArea) {
-         CkPrintf("NAN\n");
-       }
-       count++;
-      }
-    }
-    NetFEM_Elements(n,count,3,(int *)conn,"Triangles");
-    //NetFEM_Elements(n,g.nelems,3,(int *)g.conn,"Triangles");
-    NetFEM_Scalar(n,elid,1,"Element ID");
-    NetFEM_End(n);
+/*
+ Charm++ Finite-Element Framework Program:
+   Performs simple 2D structural simulation on Triangle-style inputs.
+   
+ init: 
+    read node input file
+    pass nodes into FEM framework
+    read element input file
+    pass elements into FEM 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, olawlor@acm.org
 
-    double finalArea;
-    CkPrintf("Chunk[%d]: local area: %.12f\n",rank,totalArea);
-    MPI_Reduce((void*)&totalArea,(void*)&finalArea,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
-    if(rank == 0) CkPrintf("Chunk[%d]: total area: %.12f\n",rank,finalArea);
+ Updated to new FEM interface by
+    Isaac Dooley, 2005
+ */
 
-    delete [] g.coord;
-    delete [] g.conn;
-    delete [] coord;
-    delete [] conn;
-    delete [] maptovalid;
-    delete [] nodeid;
-    delete [] elid;
-  }
-}
+#include "pgm.h"
 
 
 extern "C" void
@@ -112,11 +36,10 @@ init(void)
 {
   CkPrintf("init started\n");
   double startTime=CmiWallTimer();
-  const char *eleName="mesh1.tri";//"adpmm/xxx.1.ele";//"88mesh/mesh1.tri";
-  const char *nodeName="mesh1.node";//"adpmm/xxx.1.node";//"88mesh/mesh1.node";
+  const char *eleName="xxx.1.ele";
+  const char *nodeName="xxx.1.node";
   int nPts=0; //Number of nodes
   vector2d *pts=0; //Node coordinates
-  int *bounds;
 
   CkPrintf("Reading node coordinates from %s\n",nodeName);
   //Open and read the node coordinate file
@@ -127,11 +50,10 @@ init(void)
     fgets(line,1024,f);
     if (1!=sscanf(line,"%d",&nPts)) die("Can't read number of points!");
     pts=new vector2d[nPts];
-    bounds = new int[nPts];
     for (int i=0;i<nPts;i++) {
       int ptNo;
       if (NULL==fgets(line,1024,f)) die("Can't read node input line!");
-      if (4!=sscanf(line,"%d%lf%lf%d",&ptNo,&pts[i].x,&pts[i].y,&bounds[i])) 
+      if (3!=sscanf(line,"%d%lf%lf",&ptNo,&pts[i].x,&pts[i].y)) 
        die("Can't parse node input line!");
     }
     fclose(f);
@@ -159,12 +81,16 @@ init(void)
     }
     fclose(f);
   }
 
 
   int fem_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(fem_mesh,        // Add nodes to the current mesh
                 FEM_NODE,        // We are registering nodes
@@ -176,22 +102,13 @@ init(void)
                 FEM_DOUBLE,      // Coordinates are doubles
                 2);              // Points have dimension 2 (x,y)
  
-  CkPrintf("Passing node bounds to framework\n");
-
-
-  FEM_Mesh_data(fem_mesh,        // Add nodes to the current mesh
-                FEM_NODE,        // We are registering nodes
-                FEM_BOUNDARY,      // Register the point bound info 
-                                 // the first data elements for an FEM_NODE
-                (int *)bounds,  // The array of point bound info
-                0,               // 0 based indexing
-                nPts,            // The number of points
-                FEM_INT,        // bounds are ints
-                1);              // Points have dimension 1
 
   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(fem_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
@@ -203,52 +120,10 @@ init(void)
                 FEM_INDEX_0,     // We use zero based node numbering
                 3);              // Elements have degree 3, since triangles are defined 
                                  // by three nodes
-
-
-
-  // Register values to the elements so we can keep track of them after partitioning
-  double values[1024];
-  for(int i=0;i<1024;i++)values[i]=i;
-
-  // triangles
-  FEM_Mesh_data(fem_mesh,      // Add nodes to the current mesh
-                FEM_ELEM,      // We are registering elements with type 1
-                FEM_DATA,   
-                (int *)values,   // The array of point locations
-                0,               // 0 based indexing
-                nEle,            // The number of elements
-                FEM_DOUBLE,         // We use zero based node numbering
-                1);
  
-
-  //boundary conditions
-  FEM_Mesh_data(fem_mesh,      // Add nodes to the current mesh
-                FEM_NODE,      // We are registering elements with type 1
-                FEM_BOUNDARY,   
-                (int *)bounds,   // The array of point locations
-                0,               // 0 based indexing
-                nPts,            // The number of elements
-                FEM_INT,         // We use zero based node numbering
-                1);
-
-
-
+   
   delete[] ele;
   delete[] pts;
-  delete[] bounds;
-  
-  double *sizingData = new double[nEle];
-  for (int i=0; i<nEle; i++) sizingData[i]=-1.0;
-  FEM_Mesh_data(fem_mesh, FEM_ELEM+0, FEM_MESH_SIZING, sizingData, 0, nEle,
-                  FEM_DOUBLE, 1);
-  delete [] sizingData;
-
-  // add ghost layers
-
-  const int triangleFaces[6] = {0,1,2};
-  CkPrintf("Adding Ghost layers\n");
-  FEM_Add_ghost_layer(1,1);
-  FEM_Add_ghost_elem(0,3,triangleFaces);
 
   CkPrintf("Finished with init (Reading took %.f s)\n",CmiWallTimer()-startTime);
 
@@ -260,104 +135,25 @@ init(void)
 extern "C" void
 driver(void)
 {
-  int nnodes,nelems,nelems2,ignored;
+  int nnodes,nelems,ignored;
   int i, myId=FEM_My_partition();
   myGlobals g;
   FEM_Register(&g,(FEM_PupFn)pup_myGlobals);
   
-  _registerFEMMeshModify();
-
-  printf("partition %d is in driver\n", myId);
-
+  
   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.nnodes=nnodes;
   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);
+  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.nelems=nelems;
+  g.nnodes=nnodes; g.nelems=nelems;
   g.conn=new connRec[nelems];
-  // Get connectivity for elements
-  FEM_Mesh_data(mesh, FEM_ELEM+0, FEM_CONN, (int *)g.conn, 0, nelems, FEM_INDEX_0, 3);
-
-
-  double* tridata =new double[nelems];
-  FEM_Mesh_data(mesh, FEM_ELEM+0, FEM_DATA, (int *)tridata, 0, nelems, FEM_DOUBLE, 1);  
-
-  int nelemsghost=FEM_Mesh_get_length(mesh,FEM_ELEM+0+FEM_GHOST); 
-  double* trighostdata =new double[nelemsghost];
-  FEM_Mesh_data(mesh, FEM_ELEM+0+FEM_GHOST, FEM_DATA, (int *)trighostdata, 0, nelemsghost, FEM_DOUBLE, 1);  
-  int nnodesghost=FEM_Mesh_get_length(mesh,FEM_NODE+0+FEM_GHOST); 
-  double* nodeghostdata =new double[2*nnodesghost];
-  FEM_Mesh_data(mesh, FEM_NODE+0+FEM_GHOST, FEM_DATA, (int *)nodeghostdata, 0, nnodesghost, FEM_DOUBLE, 2);  
-
-
-  {
-       const int triangleFaces[6] = {0,1,1,2,2,0};
-       FEM_Add_elem2face_tuples(mesh, 0, 2, 3, triangleFaces);
-       FEM_Mesh_create_elem_elem_adjacency(mesh);
-       FEM_Mesh_allocate_valid_attr(mesh, FEM_ELEM+0);
-       FEM_Mesh_allocate_valid_attr(mesh, FEM_NODE);
-       FEM_Mesh_create_node_elem_adjacency(mesh);
-       FEM_Mesh_create_node_node_adjacency(mesh);
-
-       int netIndex = 0;
-       publish_data_netfem(netIndex,g); netIndex++;
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       int rank = 0;
-       MPI_Comm comm = MPI_COMM_WORLD;
-       MPI_Comm_rank(comm,&rank);
-       
-       MPI_Barrier(comm);
-       FEM_REF_INIT(mesh);
-       
-       FEM_Mesh *meshP = FEM_Mesh_lookup(FEM_Mesh_default_read(),"driver");
-       FEM_AdaptL *ada = meshP->getfmMM()->getfmAdaptL();
-       int ret_op = -1;
-
-       FEM_Adapt_Algs *adaptAlgs= meshP->getfmMM()->getfmAdaptAlgs();
-       adaptAlgs->FEM_Adapt_Algs_Init(FEM_DATA+0,FEM_DATA+4);
-       FEM_Interpolate *interp = meshP->getfmMM()->getfmInp();
-       //interp->FEM_SetInterpolateNodeEdgeFnPtr(interpolate);
-
-       MPI_Barrier(comm);
-
-       //CkPrintf("Shadow arrays have been bound\n");
-       /*
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-               
-       CkPrintf("Marking 5 nodes and one element as invalid\n");
-       FEM_set_entity_invalid(mesh, FEM_NODE, 5);
-       FEM_set_entity_invalid(mesh, FEM_NODE, 6);
-       FEM_set_entity_invalid(mesh, FEM_NODE, 7);
-       FEM_set_entity_invalid(mesh, FEM_NODE, 8);
-       FEM_set_entity_invalid(mesh, FEM_NODE, 9);      
-       FEM_set_entity_invalid(mesh, FEM_ELEM, 9);
-
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       
-       CkPrintf("Marking 5 nodes and one element as valid again\n");
-       FEM_set_entity_valid(mesh, FEM_NODE, 5);
-       FEM_set_entity_valid(mesh, FEM_NODE, 6);
-       FEM_set_entity_valid(mesh, FEM_NODE, 7);
-       FEM_set_entity_valid(mesh, FEM_NODE, 8);
-       FEM_set_entity_valid(mesh, FEM_NODE, 9);        
-       FEM_set_entity_valid(mesh, FEM_ELEM, 9);
-
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-=======
   g.S11=new double[nelems];
   g.S22=new double[nelems];
   g.S12=new double[nelems];
@@ -390,594 +186,17 @@ driver(void)
   startTime=totalStart=CkWallTimer();
   for (int t=0;t<tSteps;t++) {
     if (1) { //Structural mechanics
->>>>>>> 1.4
 
-       int adjs[3];
-       int elemid;
-       if(rank == 0) {
-         adjs[0] = 15;
-         adjs[1] = 16;
-         adjs[2] = 21; // -5;
-         elemid = 28;
-       } else if(rank == 1) {
-         adjs[0] = 19;
-         adjs[1] = 5;
-         adjs[2] = 7;
-         elemid = 21;
-       } else if(rank == 2) {
-         adjs[0] = 8;
-         adjs[1] = 11;
-         adjs[2] = 6;
-         elemid = 7;
-       } else {
-         adjs[0] = 0;
-         adjs[1] = 1;
-         adjs[2] = 2;
-         elemid = 0;
-       }
-       int newel1 = 0;
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-       //FEM_Print_n2e(mesh,adjs[0]);
-       //FEM_Print_n2e(mesh,adjs[1]);
-       //FEM_Print_n2e(mesh,adjs[2]);
-       //FEM_Print_n2n(mesh,adjs[0]);
-       //FEM_Print_n2n(mesh,adjs[1]);
-       //FEM_Print_n2n(mesh,adjs[2]);
-       //FEM_Print_e2n(mesh,newel1);
-       //FEM_Print_e2e(mesh,newel1);
-#endif
-
-       //FEM_Modify_Lock(mesh, adjs, 3, adjs, 0);
-       if(rank == 0) {
-         FEM_remove_element(mesh, elemid, 0, 1);
-       }
-       //FEM_Modify_Unlock(mesh);
-       MPI_Barrier(comm);
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-
-       //FEM_Modify_Lock(mesh, adjs, 3, adjs, 0);
-       if(rank == 0) {
-         newel1 = FEM_add_element(mesh,adjs,3,0,0);
-         CkPrintf("New Element\n");
-       }
-       //FEM_Modify_Unlock(mesh);
-       publish_data_netfem(netIndex,g); netIndex++;
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-       //FEM_Print_n2e(mesh,adjs[0]);
-       //FEM_Print_n2e(mesh,adjs[1]);
-       //FEM_Print_n2e(mesh,adjs[2]);
-       //FEM_Print_n2n(mesh,adjs[0]);
-       //FEM_Print_n2n(mesh,adjs[1]);
-       //FEM_Print_n2n(mesh,adjs[2]);
-       //FEM_Print_e2n(mesh,newel1);
-       //FEM_Print_e2e(mesh,newel1);
-#endif
-
-       if(rank==0){
-         FEM_Print_Mesh_Summary(mesh);
-         CkPrintf("%d: Removing element \n", rank);
-         
-         int nelemsghost   =FEM_Mesh_get_length(mesh,FEM_ELEM+0+FEM_GHOST); 
-         int numvalidghost =FEM_count_valid(mesh,FEM_ELEM+0+FEM_GHOST);
-         CkPrintf("nelemsghost=%d numvalidghost=%d\n", nelemsghost, numvalidghost);
-       
-         for(int i=1;i<20;i++){
-               if(FEM_is_valid(mesh, FEM_ELEM+FEM_GHOST, i)){
-                 double data[1];
-                 FEM_Mesh_data(mesh, FEM_ELEM+FEM_GHOST, FEM_DATA, (int *)data, i, 1, FEM_DOUBLE, 1);  
-
-                 CkPrintf("%d: Eating ghost element %d with value %f\n", rank, i, data[1]);
-                 int conn[3];
-                 
-                 FEM_Mesh_data(mesh, FEM_ELEM+FEM_GHOST, FEM_CONN, (int *)conn, i, 1, FEM_INDEX_0, 3);
-                 CkPrintf("conn for element is: %d %d %d\n", conn[0], conn[1], conn[2]);
-                 FEM_Modify_Lock(mesh, conn, 3, conn, 0);
-                 FEM_remove_element(mesh, FEM_From_ghost_index(i), 0, 1);
-                 FEM_Modify_Unlock(mesh);
-
-                 MPI_Barrier(comm);
-                 FEM_Print_Mesh_Summary(mesh);
-
-                 FEM_Modify_Lock(mesh, conn, 3, conn, 0);
-                 FEM_add_element(mesh, conn, 3, 0, rank); // add locally
-                 FEM_Modify_Unlock(mesh);
-                 CkPrintf("New conn for element is: %d %d %d\n", conn[0], conn[1], conn[2]);
-                 
-                 publish_data_netfem(netIndex,g); netIndex++;
-                 FEM_Print_Mesh_Summary(mesh);
-               }
-               else{
-                 //  CkPrintf("invalid element %d\n", i);
-               }
-         }
-       }
-       else {
-         CkPrintf("Rank %d\n", rank);
-         for(int i=1;i<20;i++){
-           MPI_Barrier(comm);
-           FEM_Print_Mesh_Summary(mesh);
-
-           publish_data_netfem(netIndex,g); netIndex++;
-           FEM_Print_Mesh_Summary(mesh);
-         }
-       }
-       
-       publish_data_netfem(netIndex,g); netIndex++;
-       */      
-
-       /*
-       CkPrintf("Starting Local edge flips on individual chunks\n");
-       int flip[4];
-       if(rank == 0) {
-         flip[0] = 1;
-         flip[1] = 2;
-         flip[2] = 0;
-         flip[3] = 3;
-       }
-       else if(rank == 1) {
-         flip[0] = 1;
-         flip[1] = 2;
-         flip[2] = 0;
-         flip[3] = 4;
-       }
-       else if(rank == 2) {
-         flip[0] = 13;
-         flip[1] = 14;
-         flip[2] = 15;
-         flip[3] = 7;
-       }
-       else {
-         flip[0] = 0;
-         flip[1] = 1;
-         flip[2] = 2;
-         flip[3] = 3;
-       }
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       ret_op = ada->edge_flip(flip[0],flip[1]);
-       publish_data_netfem(netIndex,g); netIndex++;
-       ret_op = ada->edge_flip(flip[2],flip[3]);
-       publish_data_netfem(netIndex,g); netIndex++;
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       
-       CkPrintf("Starting shared edge flips on individual chunks\n");
-       int sflip[4];
-       if(rank == 0) {
-         sflip[0] = 19;
-         sflip[1] = 18;
-         sflip[2] = 1;
-         sflip[3] = -4;
-       }
-       else if(rank == 1) {
-         sflip[0] = 5;
-         sflip[1] = 6;
-         sflip[2] = 7;
-         sflip[3] = -5;
-       }
-       else if(rank == 2) {
-         sflip[0] = 11;
-         sflip[1] = 2;
-         sflip[2] = 0;
-         sflip[3] = -2;
-       }
-       else {
-         sflip[0] = 0;
-         sflip[1] = 1;
-         sflip[2] = 2;
-         sflip[3] = 3;
-       }
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       ret_op = ada->edge_flip(sflip[0],sflip[1]);
-       publish_data_netfem(netIndex,g); netIndex++;
-       if(ret_op > 0) {
-         if(sflip[2]<0) sflip[2] = ret_op;
-         else if(sflip[3]<0) sflip[3] = ret_op;
-       }
-       ret_op = ada->edge_flip(sflip[2],sflip[3]);
-       publish_data_netfem(netIndex,g); netIndex++;
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       
-       CkPrintf("Starting Local edge bisect on individual chunks\n");
-       int bisect[2];
-       if(rank == 0) {
-         bisect[0] = 16;
-         bisect[1] = 21;
-       }
-       else if(rank == 1) {
-         bisect[0] = 5;
-         bisect[1] = 6;
-       }
-       else if(rank == 2) {
-         bisect[0] = 8;
-         bisect[1] = 11;
-       }
-       else {
-         bisect[0] = 0;
-         bisect[1] = 1;
-       }
-       if(rank==2) ret_op = ada->edge_bisect(bisect[0],bisect[1]);
-       publish_data_netfem(netIndex,g); netIndex++;
-       adaptAlgs->tests();
-       MPI_Barrier(comm);
-       
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       
-       CkPrintf("Starting Local vertex remove on individual chunks\n");
-       int vr[2];
-       if(rank == 0) {
-         vr[0] = ret_op;
-         vr[1] = 6;
-       }
-       else if(rank == 1) {
-         vr[0] = ret_op;
-         vr[1] = 13;
-       }
-       else if(rank == 2) {
-         vr[0] = ret_op;
-         vr[1] = 21;
-       }
-       else {
-         vr[0] = ret_op;
-         vr[1] = 1;
-       }
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       ret_op = ada->vertex_remove(vr[0],vr[1]);
-       publish_data_netfem(netIndex,g); netIndex++;
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       
-       
-       CkPrintf("Starting shared edge bisect on individual chunks\n");
-       int sbisect[2];
-       if(rank == 0) {
-         sbisect[0] = 1;//4;//21;
-         sbisect[1] = 19;//20;
-       }
-       else if(rank == 1) {
-         sbisect[0] = 0;
-         sbisect[1] = 21;
-       }
-       else if(rank == 2) {
-         sbisect[0] = 20;//9;//3;
-         sbisect[1] = 8;//2;//19;
-       }
-       else {
-         sbisect[0] = 0;
-         sbisect[1] = 1;
-       }
-       
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       ret_op = ada->edge_bisect(sbisect[0],sbisect[1]);
-       publish_data_netfem(netIndex,g); netIndex++;
-       CkPrintf("Starting shared vertex remove on individual chunks\n");
-
-       int svr[2];
-       if(rank == 0) {
-         svr[0] = ret_op;
-         svr[1] = 19;
-       }
-       else if(rank == 1) {
-         svr[0] = ret_op;
-         svr[1] = 21;
-       }
-       else if(rank == 2) {
-         svr[0] = ret_op;
-         svr[1] = 20;
-       }
-       else {
-         svr[0] = ret_op;
-         svr[1] = 1;
-       }
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       ret_op = ada->vertex_remove(svr[0],svr[1]);
-       publish_data_netfem(netIndex,g); netIndex++;
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-
-       CkPrintf("Starting Local edge contract on individual chunks\n");
-       int contract[2];
-       if(rank == 0) {
-         contract[1] = 16;
-         contract[0] = 21;
-       }
-       else if(rank == 1) {
-         contract[0] = 5;
-         contract[1] = 6;
-       }
-       else if(rank == 2) {
-         contract[0] = 19;
-         contract[1] = 17;
-       }
-       else {
-         contract[0] = 0;
-         contract[1] = 1;
-       }
-       
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       if(rank==0) ret_op = ada->edge_contraction(contract[0],contract[1]);
-       publish_data_netfem(netIndex,g); netIndex++;
-
-       //CkPrintf("Starting Local vertex split on individual chunks\n");
-       /*
-       int vs[3];
-       if(rank == 0) {
-         vs[0] = ret_op;
-         vs[1] = 9;
-         vs[2] = 13;
-       }
-       else if(rank == 1) {
-         vs[0] = ret_op;
-         vs[1] = 8;
-         vs[2] = 7;
-       }
-       else if(rank == 2) {
-         vs[0] = ret_op;
-         vs[1] = 14;
-         vs[2] = 23;
-       }
-       else {
-         vs[0] = ret_op;
-         vs[1] = 2;
-         vs[2] = 3;
-       }
-#ifdef SUMMARY_ON
-       //FEM_Print_Mesh_Summary(mesh);
-#endif
-       //ret_op = ada->vertex_split(vs[0],vs[1],vs[2]);
-       //publish_data_netfem(netIndex,g); netIndex++;
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       */
-       /*
-       CkPrintf("Starting shared edge contract on individual chunks\n");
-       int scontract[2];
-       if(rank == 0) {
-         scontract[0] = 9;
-         scontract[1] = 10;
-       }
-       else if(rank == 1) {
-         scontract[0] = 5;
-         scontract[1] = 6;
-       }
-       else if(rank == 2) {
-         scontract[0] = 11;
-         scontract[1] = 2;
-       }
-       else {
-         scontract[0] = 0;
-         scontract[1] = 1;
-       }
-       
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       ret_op = ada->edge_contraction(scontract[0],scontract[1]);
-       publish_data_netfem(netIndex,g); netIndex++;
-       /*
-       //CkPrintf("Starting shared vertex split on individual chunks\n");
-       int svs[3];
-       if(rank == 0) {
-         svs[0] = ret_op;
-         svs[1] = 1;
-         svs[2] = -6;
-       }
-       else if(rank == 1) {
-         svs[0] = ret_op;
-         svs[1] = 7;
-         svs[2] = 7;
-       }
-       else if(rank == 2) {
-         svs[0] = ret_op;
-         svs[1] = 0;
-         svs[2] = -2;
-       }
-       else {
-         svs[0] = ret_op;
-         svs[1] = 2;
-         svs[2] = 3;
-       }
-#ifdef SUMMARY_ON
-       //FEM_Print_Mesh_Summary(mesh);
-#endif
-       //ret_op = ada->vertex_split(svs[0],svs[1],svs[2]);
-       //publish_data_netfem(netIndex,g); netIndex++;
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       */
-
-       /*
-       CkPrintf("Starting LEB on individual chunks\n");
-       int *leb_elem = (int*)malloc(1*sizeof(int));
-       if(rank==0) {
-         leb_elem[0] = 2;
-       }
-       else if(rank==1) {
-         leb_elem[0] = 13; //4;
-       }
-       else if(rank==2) {
-         leb_elem[0] = 20; //26;
-       }
-       else if (rank == 3){
-         leb_elem[0] = 14;
-       }
-       else {
-         leb_elem[0] = 0;
-       }
-
-       adaptAlgs->refine_element_leb(leb_elem[0]);
-       publish_data_netfem(netIndex,g); netIndex++;
-       */
-       /*
-         int nEle;
-         //for(int tstep = 0; tstep < 2; tstep++) {
-         nEle = FEM_Mesh_get_length(mesh, FEM_ELEM);   
-         for (int i=0; i<nEle; i++)
-         if (FEM_is_valid(mesh, FEM_ELEM, i))
-         adaptAlgs->refine_element_leb(i);
-         publish_data_netfem(netIndex,g); netIndex++;
-         FEM_Print_Mesh_Summary(mesh);
-         //}
-         */
-
-      
-      double targetArea = 0.00004;
-      
-      for(int tstep = 0; tstep < 1; tstep++) {
-       adaptAlgs->simple_refine(targetArea);
-       publish_data_netfem(netIndex,g); netIndex++;
-       adaptAlgs->tests();
-       MPI_Barrier(comm);
-       
-       //int *nodes = new int[g.nnodes];
-       //for (int i=0; i<g.nnodes; i++) nodes[i]=i;    
-       //FEM_mesh_smooth(mesh, nodes, g.nnodes, FEM_DATA+0);
-       //publish_data_netfem(netIndex,g); netIndex++;
-       //delete [] nodes;
-
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       adaptAlgs->simple_coarsen(targetArea);
-       publish_data_netfem(netIndex,g); netIndex++;
-       adaptAlgs->tests();
-       MPI_Barrier(comm);
-       
-       //int *nodes = new int[g.nnodes];
-       //for (int i=0; i<g.nnodes; i++) nodes[i]=i;
-       //FEM_mesh_smooth(mesh, nodes, g.nnodes, FEM_DATA+0);
-       //publish_data_netfem(netIndex,g); netIndex++;
-       //delete [] nodes;
-
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-      }
+    //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);
 
-      targetArea = 0.00004;
-      
-      for(int tstep = 0; tstep < 4; tstep++) {
-       adaptAlgs->simple_refine(targetArea);
-       publish_data_netfem(netIndex,g); netIndex++;
-       adaptAlgs->tests();
-       MPI_Barrier(comm);
-       
-       //int *nodes = new int[g.nnodes];
-       //for (int i=0; i<g.nnodes; i++) nodes[i]=i;    
-       //FEM_mesh_smooth(mesh, nodes, g.nnodes, FEM_DATA+0);
-       //publish_data_netfem(netIndex,g); netIndex++;
-       //delete [] nodes;
+    //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);
+    }
 
-       targetArea *= 0.5;
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-      }
-      targetArea /= 0.4;
-      
-      for(int tstep = 0; tstep < 4; tstep++) {
-       adaptAlgs->simple_coarsen(targetArea);
-       publish_data_netfem(netIndex,g); netIndex++;
-       adaptAlgs->tests();
-       MPI_Barrier(comm);
-       
-       //int *nodes = new int[g.nnodes];
-       //for (int i=0; i<g.nnodes; i++) nodes[i]=i;
-       //FEM_mesh_smooth(mesh, nodes, g.nnodes, FEM_DATA+0);
-       //publish_data_netfem(netIndex,g); netIndex++;
-       //delete [] nodes;
-
-       targetArea *= 1.5;
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-      }
-
-      //wave propagation on a bar
-      targetArea = 0.00004;
-      double xmin = 0.00;
-      double xmax = 0.1;
-      double ymin = 0.00;
-      double ymax = 0.01;
-      for(int tstep = 0; tstep < 0; tstep++) {
-       targetArea = 0.000002;
-       adaptAlgs->simple_refine(targetArea, xmin, ymin, xmax, ymax);
-       publish_data_netfem(netIndex,g); netIndex++;
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       targetArea = 0.0000014;
-       adaptAlgs->simple_coarsen(targetArea, xmin, ymin, xmax, ymax);
-       publish_data_netfem(netIndex,g); netIndex++;
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       ymin += 0.01;
-       ymax += 0.01;
-      }
-
-      //crack propagation on a block
-      targetArea = 0.00004;
-      xmin = 0.00;
-      xmax = 0.2;
-      double xcrackmin = 0.09;
-      double xcrackmax = 0.10;
-      ymin = 0.00;
-      ymax = 0.02;
-      for(int tstep = 0; tstep < 0; tstep++) {
-       targetArea = 0.000025;
-       adaptAlgs->simple_refine(targetArea, xmin, ymin, xmax, ymax);
-       publish_data_netfem(netIndex,g); netIndex++;
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       targetArea = 0.00005;
-       adaptAlgs->simple_coarsen(targetArea, xmin, ymin, xmax, ymax);
-       //publish_data_netfem(netIndex,g); netIndex++;
-#ifdef SUMMARY_ON
-       FEM_Print_Mesh_Summary(mesh);
-#endif
-       /*if(tstep > 2) {
-         targetArea = 0.000025;
-         adaptAlgs->simple_refine(targetArea, xcrackmin, ymin, xcrackmax, ymax);
-         //publish_data_netfem(netIndex,g); netIndex++;
-#ifdef SUMMARY_ON
-         FEM_Print_Mesh_Summary(mesh);
-#endif
-         xcrackmin -= 0.004;
-         xcrackmax += 0.004;
-       }
-       */
-
-       ymin += 0.02;
-       ymax += 0.02;
-      }
-      /*
     //Debugging/perf. output
     double curTime=CkWallTimer();
     double total=curTime-startTime;
@@ -993,7 +212,7 @@ driver(void)
              }
            }
     }
-    // perform migration-based load balancing 
+    /* perform migration-based load balancing */
     if (t%1024==0)
       FEM_Migrate();
     
@@ -1012,19 +231,11 @@ driver(void)
            NetFEM_End(n);
     }
   }
-*/
 
-      CkPrintf("chunk %d Waiting for Synchronization\n",rank);
-      MPI_Barrier(comm);
-      CkPrintf("Synchronized\n");
-#ifdef SUMMARY_ON
-      FEM_Print_Mesh_Summary(mesh);
-#endif
-      publish_data_netfem(netIndex,g); netIndex++;
-      
-      CkExit();
+  if (myId==0) {
+    double elapsed=CkWallTimer()-totalStart;
+    CkPrintf("Driver finished: average %.6f s/step\n",elapsed/tSteps);
   }
-  
 }
 
 
@@ -1039,149 +250,78 @@ void pup_myGlobals(pup_er p,myGlobals *g)
   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];
+    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);
+  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;
+    delete[] g->R_net;
+    delete[] g->d;
+    delete[] g->v;
+    delete[] g->a;
+       delete[] g->S11;
+       delete[] g->S22;
+       delete[] g->S12;
   }
 }
 
 
-void FEM_mesh_smooth(int mesh, int *nodes, int nNodes, int attrNo)
-{
-  vector2d *centroids, newPos, *coords, *ghostCoords, *vGcoords;
-  int nEle, nGn, *boundVals, nodesInChunk, nVg;
-  int neighbors[3], *adjelems;
-  int gIdxN;
-  int j=0;
-  double x[3], y[3];
-  FEM_Mesh *meshP = FEM_Mesh_lookup(mesh, "driver");
 
-  nodesInChunk = FEM_Mesh_get_length(mesh,FEM_NODE);
-  boundVals = new int[nodesInChunk];
-  nGn = FEM_Mesh_get_length(mesh, FEM_GHOST + FEM_NODE);
-  coords = new vector2d[nodesInChunk+nGn];
+//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);
 
-  FEM_Mesh_data(mesh, FEM_NODE, FEM_BOUNDARY, (int*) boundVals, 0, nodesInChunk, FEM_INT, 1);    
+  const double shearForce=1.0e-6/(dt*dt*xm);
 
-  FEM_Mesh_data(mesh, FEM_NODE, attrNo, (double*)coords, 0, nodesInChunk, FEM_DOUBLE, 2);
-  for (int i=0; i<(nodesInChunk); i++) {
-    //CkPrintf(" coords[%d]: (%f, %f)\n", i, coords[i].x, coords[i].y);
-  }
-  IDXL_Layout_t coord_layout = IDXL_Layout_create(IDXL_DOUBLE, 2);
-  FEM_Update_ghost_field(coord_layout,-1, coords); 
-  ghostCoords = &(coords[nodesInChunk]);
-  /*
-  for (int i=0; i<nGn;i++) {
-    if (FEM_is_valid(mesh, FEM_GHOST+FEM_NODE, i)) {
-      CkPrintf("ghost %d is valid \n", i);       
-      // vGcoords[j]=ghostCoords[i];
-      //j++;
+  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;
     }
-    else
-      CkPrintf("ghost %d is invalid \n", i);
-  }
-  */
-  for (int i=0; i<(nodesInChunk+nGn); i++) {
-    //CkPrintf(" coords[%d]: (%f, %f)\n", i, coords[i].x, coords[i].y);
-  }
-//  FEM_Mesh_data(FEM_Mesh_default_write(), FEM_GHOST+FEM_NODE, attrNo, (double*)ghostCoords, 0, nGn, FEM_DOUBLE, 2);
-  for (int i=0; i<nNodes; i++)
-  {
-    newPos.x=0;
-    newPos.y=0;
-    CkAssert(nodes[i]<nodesInChunk);    
-    if (FEM_is_valid(mesh, FEM_NODE, i) && boundVals[i]>-1) //node must be internal
-    {
-      meshP->n2e_getAll(i, &adjelems, &nEle);
-      centroids = new vector2d[nEle];
-      
-      for (int j=0; j<nEle; j++) { //for all adjacent elements, find centroids
-       meshP->e2n_getAll(adjelems[j], neighbors);
-       for (int k=0; k<3; k++) {
-         if (neighbors[k]<-1) {
-           gIdxN = FEM_From_ghost_index(neighbors[k]);
-           x[k] = ghostCoords[gIdxN].x;
-           y[k] = ghostCoords[gIdxN].y;
-         }
-         else {
-           x[k] = coords[neighbors[k]].x;
-           y[k] = coords[neighbors[k]].y;
-         }
-       }     
-       centroids[j].x=(x[0]+x[1]+x[2])/3.0;
-       centroids[j].y=(y[0]+y[1]+y[2])/3.0;
-       newPos.x += centroids[j].x;
-       newPos.y += centroids[j].y;
-      }
-      newPos.x/=nEle;
-      newPos.y/=nEle;
-      FEM_set_entity_coord2(mesh, FEM_NODE, nodes[i], newPos.x, newPos.y);
-      delete [] centroids;
-      delete [] adjelems;
+#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
   }
-  delete [] coords;
-  delete [] boundVals;
-}
-
-void interpolate(FEM_Interpolate::NodalArgs args, FEM_Mesh *meshP)
-{
-  //CkPrintf("INTERPOLATOR!!!!!!!!!!!\n");
-  int length = meshP->node.realsize();
-  int *boundVals= new int[length];
-
-  FEM_Mesh_dataP(meshP, FEM_NODE, FEM_BOUNDARY, (int*) boundVals, 0, length , FEM_INT, 1);   
-  CkVec<FEM_Attribute *>*attrs = (meshP->node).getAttrVec();
-  for (int i=0; i<attrs->size(); i++) {
-    FEM_Attribute *a = (FEM_Attribute *)(*attrs)[i];
-    if (a->getAttr() < FEM_ATTRIB_TAG_MAX || a->getAttr()==FEM_BOUNDARY) {
-      if (a->getAttr()==FEM_BOUNDARY) {
-       
-       int n1_bound =boundVals[args.nodes[0]]; 
-       int n2_bound =boundVals[args.nodes[1]]; 
-       if (n1_bound == n2_bound && n1_bound < 0) {
-         a->copyEntity(args.n, *a, args.nodes[0]);
-       } else if (n1_bound != n2_bound && n1_bound<0 && n2_bound < 0){
-         a->copyEntity(args.n, *a, args.nodes[0]); //a node which is not on the boundary
-       } else  if(n1_bound < 0){
-         a->copyEntity(args.n, *a, args.nodes[1]); 
-       } else {
-         a->copyEntity(args.n, *a, args.nodes[0]); 
-       }
-      }
-      else {
-       FEM_DataAttribute *d = (FEM_DataAttribute *)a;
-       d->interpolate(args.nodes[0], args.nodes[1], args.n, args.frac);
-      }
-    }
+  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!");
   }
-  delete boundVals;
 }
-
-
-
index 9caa60dc3459f7572527dfc282a3d7a60c62064e..b8a195181b092bc508fec0f6fa781d1ec8f7bd96 100644 (file)
@@ -5,12 +5,9 @@
 #include "fem.h"
 #include "netfem.h"
 #include "vector2d.h"
-#include "fem_interpolate.h"
-
 
 //One element's connectivity information
 typedef int connRec[3];
-typedef int connRec2[4];
 
 // A structure for handling data that may need to be migrated
 struct myGlobals {
@@ -19,17 +16,11 @@ struct myGlobals {
   vector2d *coord;
   connRec *conn;
 
-  int nelems2;
-  connRec2 *conn2;
-
   vector2d *R_net, *d, *v, *a;
   
   double *S11, *S22, *S12;
 };
 
-void FEM_mesh_smooth(int mesh, int *nodes, int nNodes, int attrNo);
-
-void interpolate(FEM_Interpolate::NodalArgs args,FEM_Mesh *meshP);
 //Compute forces on constant-strain triangles:
 void CST_NL(const vector2d *coor,const connRec *lm,vector2d *R_net,
            const vector2d *d,const double *c,