Charm++: Charm Source Code Documentation

00001 /*
00002 Triangle Mesh Refinement (TMR) demo code.
00003 
00004 Reads in a Triangle mesh.
00005 Passes it to the REFINE2D subsystem.
00006 Asks for refinements occasionally.
00007 */
00008 #include <stdlib.h>
00009 #include <unistd.h>
00010 #include <stdio.h>
00011 #include <math.h>
00012 #include "charm++.h"
00013 #include "fem.h"
00014 /*My HEADERS*/
00015 #include "ckvector3d.h"
00016 #include "charm-api.h"
00017 #include "fem_mesh.h"
00018 #include "netfem.h"
00019 #include "refine.h"
00020 #include "femrefine.h"
00021 #include "pgm.h"
00022 
00023 //The material constants c, as computed by fortran mat_const
00024 // I think the units here are Pascals (N/m^2)
00025 const double matConst[4]={3.692e9,  1.292e9,  3.692e9,  1.200e9 };
00026 
00027 //Material density, Kg/m^3
00028 const double density=5.0*1000.0;
00029 //Plate thickness, meters
00030 const double thickness=0.0001;
00031 
00032 //The timestep, in seconds
00033 // This aught to be adjusted when the mesh changes!
00034 const double dt=1.0e-12;
00035 
00036 static void die(const char *str) {
00037   CkError("Fatal error: %s\n",str);
00038   CkExit();
00039 }
00040 
00042 void resize_nodes(void *data,int *len,int *max);
00043 void resize_elems(void *data,int *len,int *max);
00044 
00045 #define NANCHECK 1 /*Check for NaNs at each timestep*/
00046 
00047 extern "C" void
00048 init(void)
00049 {
00050   CkPrintf("init started\n");
00051 
00052   const char *eleName="out.1024.ele";
00053   const char *nodeName="out.1024.node";
00054   int nPts=0; //Number of nodes
00055   vector2d *pts=0; //Node coordinates
00056 
00057   CkPrintf("Reading node coordinates from %s\n",nodeName);
00058   //Open and read the node coordinate file
00059   {
00060     char line[1024];
00061     FILE *f=fopen(nodeName,"r");
00062     if (f==NULL) die("Can't open node file!");
00063     fgets(line,1024,f);
00064     if (1!=sscanf(line,"%d",&nPts)) die("Can't read number of points!");
00065     pts=new vector2d[nPts];
00066     for (int i=0;i<nPts;i++) {
00067       int ptNo;
00068       if (NULL==fgets(line,1024,f)) die("Can't read node input line!");
00069       if (3!=sscanf(line,"%d%lf%lf",&ptNo,&pts[i].x,&pts[i].y)) 
00070     die("Can't parse node input line!");
00071     }
00072     fclose(f);
00073   }
00074   CkPrintf("Passing node coords to framework\n");
00075   
00076   // Register the node entity and its data arrays that will be used later. This
00077   // needs to be done so that the width of the data segments are set correctly 
00078   // and can be used later
00079   FEM_Register_entity(FEM_Mesh_default_write(),FEM_NODE,NULL,nPts,nPts,resize_nodes);
00080   for(int k=0;k<=4;k++){
00081     if(k != 0){
00082       vector2d *t = new vector2d[nPts];
00083       FEM_Register_array(FEM_Mesh_default_write(),FEM_NODE,FEM_DATA+k,t,FEM_DOUBLE,2);
00084     }else{
00085       FEM_Register_array(FEM_Mesh_default_write(),FEM_NODE,FEM_DATA+k,pts,FEM_DOUBLE,2);
00086     }
00087   }
00088   double *td = new double[nPts];
00089   FEM_Register_array(FEM_Mesh_default_write(),FEM_NODE,FEM_DATA+5,td,FEM_DOUBLE,1);
00090   int *validNodes = new int[nPts];
00091   for(int ii=0;ii<nPts;ii++){
00092     validNodes[ii]=1;
00093   }
00094   FEM_Register_array(FEM_Mesh_default_write(),FEM_NODE,FEM_VALID,validNodes,FEM_INT,1);
00095   
00096   int nEle=0;
00097   int *ele=NULL;
00098   CkPrintf("Reading elements from %s\n",eleName);
00099   //Open and read the element connectivity file
00100   {
00101     char line[1024];
00102     FILE *f=fopen(eleName,"r");
00103     if (f==NULL) die("Can't open element file!");
00104     fgets(line,1024,f);
00105     if (1!=sscanf(line,"%d",&nEle)) die("Can't read number of elements!");
00106     ele=new int[3*nEle];
00107     for (int i=0;i<nEle;i++) {
00108       int elNo;
00109       if (NULL==fgets(line,1024,f)) die("Can't read element input line!");
00110       if (4!=sscanf(line,"%d%d%d%d",&elNo,&ele[3*i+0],&ele[3*i+1],&ele[3*i+2])) 
00111     die("Can't parse element input line!");
00112       ele[3*i+0]--; //Fortran to C indexing
00113       ele[3*i+1]--; //Fortran to C indexing
00114       ele[3*i+2]--; //Fortran to C indexing
00115     }
00116     fclose(f);
00117   }
00118   
00119   CkPrintf("Passing elements to framework\n");
00120   // Register the Element entity and its connectivity array. Register the
00121   // data arrays to set up the widths correctly at the beginning
00122   FEM_Register_entity(FEM_Mesh_default_write(),FEM_ELEM,NULL,nEle,nEle,resize_elems);
00123   FEM_Register_array(FEM_Mesh_default_write(),FEM_ELEM,FEM_CONN,ele,FEM_INDEX_0,3);
00124   
00125   for(int k=0;k<3;k++){
00126     void *t = new double[nEle];
00127     FEM_Register_array(FEM_Mesh_default_write(),FEM_ELEM,FEM_DATA+k,t,FEM_DOUBLE,1);
00128   }
00129   int *validElem = new int[nEle];
00130   for(int ii=0;ii<nEle;ii++){
00131     validElem[ii]=1;
00132   }
00133   FEM_Register_array(FEM_Mesh_default_write(),FEM_ELEM,FEM_VALID,validElem,FEM_INT,1);
00134   
00135   /*Build the ghost layer for refinement border*/
00136   FEM_Add_ghost_layer(2,0); /*2 nodes/tuple, do not add ghost nodes*/
00137   //FEM_Add_ghost_layer(2,1); /*2 nodes/tuple, do not add ghost nodes*/
00138   const static int tri2edge[6]={0,1, 1,2, 2,0};
00139   FEM_Add_ghost_elem(0,3,tri2edge);
00140   CkPrintf("Finished with init\n");
00141 }
00142 
00143 struct myGlobals {
00144   int nnodes,maxnodes;
00145   int nelems,maxelems;
00146   int *conn; //Element connectivity table
00147   vector2d *coord; //Undeformed coordinates of each node
00148   vector2d *R_net, *d, *v, *a; //Physical fields of each node
00149   double *m_i; //Inverse of mass at each node
00150   int m_i_fid; //Field ID for m_i
00151   int *validNode,*validElem;
00152   double *S11, *S22, *S12; //Stresses for each element
00153 };
00154 
00155 void pup_myGlobals(pup_er p,myGlobals *g) 
00156 {
00157   FEM_Print("-------- called pup routine -------");
00158   pup_int(p,&g->nnodes);
00159   pup_int(p,&g->nelems);
00160   pup_int(p,&g->maxelems);
00161   pup_int(p,&g->maxnodes);
00162   int nnodes=g->nnodes, nelems=g->nelems;
00163   if (pup_isUnpacking(p)) {
00164     g->coord=new vector2d[g->maxnodes];
00165     g->conn=new int[3*g->maxelems];
00166     g->R_net=new vector2d[g->maxnodes]; //Net force
00167     g->d=new vector2d[g->maxnodes];//Node displacement
00168     g->v=new vector2d[g->maxnodes];//Node velocity
00169     g->a=new vector2d[g->maxnodes];
00170     g->m_i=new double[g->maxnodes];
00171     g->S11=new double[g->maxelems];
00172     g->S22=new double[g->maxelems];
00173     g->S12=new double[g->maxelems];
00174     g->validNode = new int[g->maxnodes];
00175     g->validElem = new int[g->maxelems];
00176   }
00177   pup_doubles(p,(double *)g->coord,2*nnodes);
00178   pup_ints(p,(int *)g->conn,3*nelems);
00179   pup_doubles(p,(double *)g->R_net,2*nnodes);
00180   pup_doubles(p,(double *)g->d,2*nnodes);
00181   pup_doubles(p,(double *)g->v,2*nnodes);
00182   pup_doubles(p,(double *)g->a,2*nnodes);
00183   pup_doubles(p,(double *)g->m_i,nnodes);
00184   pup_doubles(p,(double *)g->S11,nelems);
00185   pup_doubles(p,(double *)g->S22,nelems);
00186   pup_doubles(p,(double *)g->S12,nelems);
00187   pup_ints(p,(int *)g->validNode,nnodes);
00188   pup_ints(p,(int *)g->validElem,nelems);
00189   if (pup_isDeleting(p)) {
00190     delete[] g->coord;
00191     delete[] g->conn;
00192     delete[] g->R_net;
00193     delete[] g->d;
00194     delete[] g->v;
00195     delete[] g->a;
00196     delete[] g->m_i;
00197     delete[] g->S11;
00198     delete[] g->S22;
00199     delete[] g->S12;
00200     delete[] g->validNode;
00201     delete[] g->validElem;
00202   }
00203 }
00204 
00205 //Return the signed area of triangle i
00206 double calcArea(myGlobals &g, int i)
00207 {
00208   int n1=g.conn[3*i+0];
00209   int n2=g.conn[3*i+1];
00210   int n3=g.conn[3*i+2];
00211   vector2d a=g.coord[n1];
00212   vector2d b=g.coord[n2];
00213   vector2d c=g.coord[n3];
00214   c-=a; b-=a;
00215   double area=0.5*(b.x*c.y-c.x*b.y);
00216   return area;
00217 }
00218 
00219 // Check the quality of triangle i
00220 void checkTriangle(myGlobals &g, int i)
00221 {
00222   double area=calcArea(g,i);
00223   if (area<0) {
00224     CkError("Triangle %d of chunk %d is inverted! (area=%g)\n",
00225         i,FEM_My_partition(),area);
00226     CkAbort("Inverted triangle");
00227   }
00228   if (area<1.0e-15) {
00229     CkError("Triangle %d of chunk %d is a sliver!\n",i,FEM_My_partition());
00230     CkAbort("Sliver triangle");
00231   }
00232 }
00233 
00234 //Compute forces on constant-strain triangles:
00235 void CST_NL(const vector2d *coor,const int *lm,vector2d *R_net,
00236         const vector2d *d,const double *c,
00237         int numnp,int numel,
00238         double *S11o,double *S22o,double *S12o);
00239 
00240 //Update node position, velocity, accelleration based on net force.
00241 void advanceNodes(const double dt,int nnodes,const vector2d *coord,
00242           vector2d *R_net,vector2d *a,vector2d *v,vector2d *d,
00243           const double *m_i,bool dampen)
00244 {
00245   const vector2d z(0,0);
00246   const double shearForce=1.0e-11/(dt*dt);
00247   bool someNaNs=false;
00248   int i;
00249   for (i=0;i<nnodes;i++) {
00250     vector2d R_n=R_net[i];
00251 #if NANCHECK
00252     if (((R_n.x-R_n.x)!=0)) {
00253       CkPrintf("R_net[%d]=NaN at (%.4f,%.4f)   ",i,coord[i].x,coord[i].y);
00254       CmiAbort("nan node");
00255       someNaNs=true;
00256     }
00257     if (fabs(d[i].x)>1.0) {
00258       CkPrintf("d[%d] %f large at (%.4f,%.4f)   ",i,d[i].x,coord[i].x,coord[i].y);
00259       someNaNs=true;
00260     }
00261 #endif
00262     R_net[i]=z;
00263     //Apply boundary conditions (HACK: hardcoded!)
00264     if (1) {
00265       if (coord[i].x<0.00001)
00266     R_n.y+=shearForce/m_i[i]; //Bottom edge pushed hard down
00267       if (coord[i].y>0.02-0.00001)
00268     R_n=z; //Top edge held in place
00269     }
00270     //Update displacement and velocity
00271     vector2d aNew=R_n*m_i[i];
00272     v[i]+=(dt*0.5)*(aNew+a[i]);
00273     d[i]+=dt*v[i]+(dt*dt*0.5)*aNew;
00274     a[i]=aNew;   
00275     //if (coord[i].y>0.02-0.00001) d[i].y=0.0; //Top edge in horizontal slot
00276   }
00277   if (dampen)
00278     for (i=0;i<nnodes;i++)
00279       v[i]*=0.9; //Dampen velocity slightly (prevents eventual blowup)
00280   
00281   if (someNaNs) {
00282     CkPrintf("Nodes all NaN!\n");
00283     CkAbort("Node forces NaN!");
00284   }
00285 }
00286 
00287 //Fill out the m_i array with the inverse of the node masses
00288 void calcMasses(myGlobals &g) {
00289   int i;
00290   double *m_i=g.m_i;
00291   //Zero out node masses
00292   for (i=0;i<g.nnodes;i++) m_i[i]=0.0;
00293   //Add mass from surrounding triangles:
00294   for (i=0;i<g.nelems;i++) {
00295     if(g.validElem[i]){
00296       int n1=g.conn[3*i+0];
00297       int n2=g.conn[3*i+1];
00298       int n3=g.conn[3*i+2];
00299       double area=calcArea(g,i);
00300       //        if (1 || i%100==0) CkPrintf("Triangle %d (%d %d %d) has area %.3g\n",i,n1,n2,n3,area);
00301       double mass=0.333*density*(thickness*area);
00302       m_i[n1]+=mass;
00303       m_i[n2]+=mass;
00304       m_i[n3]+=mass;
00305     }
00306   }
00307   //Include mass from other processors
00308   FEM_Update_field(g.m_i_fid,m_i);
00309   //Invert masses to get m_i
00310   for (i=0;i<g.nnodes;i++) {
00311     double mass=m_i[i];
00312     if (mass<1.0e-10) m_i[i]=1.0; //Disconnected node (!)
00313     else m_i[i]=1.0/mass;
00314   }
00315 }
00316 
00317 void init_myGlobal(myGlobals *g){
00318   g->coord = g->R_net = g->d = g->v = g->a = NULL;
00319   g->m_i = NULL;
00320   g->conn = NULL;
00321   g->S11 = g->S22 = g->S12 = NULL;
00322 }
00323 
00324 
00325 void resize_nodes(void *data,int *len,int *max){
00326   printf("[%d] resize nodes called len %d max %d\n",FEM_My_partition(),*len,*max);
00327   FEM_Register_entity(FEM_Mesh_default_read(),FEM_NODE,data,*len,*max,resize_nodes);
00328   myGlobals *g = (myGlobals *)data;
00329   vector2d *coord=g->coord,*R_net=g->R_net,*d=g->d,*v=g->v,*a=g->a;
00330   double *m_i=g->m_i;
00331   int *validNode = g->validNode;
00332   
00333   g->coord=new vector2d[*max];
00334   g->coord[0].x = 0.9;
00335   g->coord[0].y = 0.8;
00336   g->maxnodes = *max;
00337   g->R_net=new vector2d[g->maxnodes]; //Net force
00338   g->d=new vector2d[g->maxnodes];//Node displacement
00339   g->v=new vector2d[g->maxnodes];//Node velocity
00340   g->a=new vector2d[g->maxnodes];//Node accelleration
00341   g->m_i=new double[g->maxnodes];//Node mass
00342   g->validNode = new int[g->maxnodes]; //is the node valid
00343   
00344   if(coord != NULL){
00345     for(int k=0;k<*len;k++){
00346       printf("before resize node %d ( %.6f %.6f ) \n",k,coord[k].x,coord[k].y);
00347     }
00348   } 
00349   
00350   FEM_Register_array(FEM_Mesh_default_read(),FEM_NODE,FEM_DATA,(void *)g->coord,FEM_DOUBLE,2);
00351   FEM_Register_array(FEM_Mesh_default_read(),FEM_NODE,FEM_DATA+1,(void *)g->R_net,FEM_DOUBLE,2);
00352   FEM_Register_array(FEM_Mesh_default_read(),FEM_NODE,FEM_DATA+2,(void *)g->d,FEM_DOUBLE,2);
00353   FEM_Register_array(FEM_Mesh_default_read(),FEM_NODE,FEM_DATA+3,(void *)g->v,FEM_DOUBLE,2);
00354   FEM_Register_array(FEM_Mesh_default_read(),FEM_NODE,FEM_DATA+4,(void *)g->a,FEM_DOUBLE,2);
00355   FEM_Register_array_layout(FEM_Mesh_default_read(),FEM_NODE,FEM_DATA+5,(void *)g->m_i,g->m_i_fid);
00356   FEM_Register_array(FEM_Mesh_default_read(),FEM_NODE,FEM_VALID,(void *)g->validNode,FEM_INT,1);
00357   
00358   for(int k=0;k<*len;k++){
00359     printf("after resize node %d ( %.6f %.6f )\n",k,g->coord[k].x,g->coord[k].y);
00360   }
00361   
00362   if(coord != NULL){
00363     delete [] coord;
00364     delete [] R_net;
00365     delete [] d;
00366     delete [] v;
00367     delete [] a;
00368     delete [] m_i;
00369     delete [] validNode;
00370   }
00371 };
00372 
00373 void resize_elems(void *data,int *len,int *max){
00374   printf("[%d] resize elems called len %d max %d\n",FEM_My_partition(),*len,*max);
00375   FEM_Register_entity(FEM_Mesh_default_read(),FEM_ELEM,data,*len,*max,resize_elems);
00376   myGlobals *g = (myGlobals *)data;
00377   int *conn=g->conn;
00378   double *S11 = g->S11,*S22 = g->S22,*S12 = g->S12;
00379   int *validElem = g->validElem;
00380   
00381   g->conn = new int[3*(*max)];
00382   g->maxelems = *max;
00383   g->S11=new double[g->maxelems];
00384   g->S22=new double[g->maxelems];
00385   g->S12=new double[g->maxelems];
00386   g->validElem = new int[g->maxelems];
00387   
00388   FEM_Register_array(FEM_Mesh_default_read(),FEM_ELEM,FEM_CONN,(void *)g->conn,FEM_INDEX_0,3);  
00389   CkPrintf("Connectivity array starts at %p \n",g->conn);
00390   FEM_Register_array(FEM_Mesh_default_read(),FEM_ELEM,FEM_DATA,(void *)g->S11,FEM_DOUBLE,1);    
00391   FEM_Register_array(FEM_Mesh_default_read(),FEM_ELEM,FEM_DATA+1,(void *)g->S22,FEM_DOUBLE,1);  
00392   FEM_Register_array(FEM_Mesh_default_read(),FEM_ELEM,FEM_DATA+2,(void *)g->S12,FEM_DOUBLE,1);  
00393   FEM_Register_array(FEM_Mesh_default_read(),FEM_ELEM,FEM_VALID,(void *)g->validElem,FEM_INT,1);    
00394   
00395   if(conn != NULL){
00396     delete [] conn;
00397     delete [] S11;
00398     delete [] S22;
00399     delete [] S12;
00400     delete [] validElem;
00401   }
00402 };
00403 
00404 void repeat_after_split(void *data){
00405   myGlobals *g = (myGlobals *)data;
00406   g->nelems = FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_ELEM);
00407   g->nnodes = FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_NODE);
00408   for(int k=0;k<g->nnodes;k++){
00409     if(g->validNode[k]){
00410       printf(" node %d ( %.6f %.6f )\n",k,g->coord[k].x,g->coord[k].y);
00411     }   
00412   }
00413   calcMasses(*g);
00414 };
00415 
00416 
00417 extern "C" void
00418 driver(void)
00419 {
00420   int ignored;
00421   int i;  
00422   int myChunk=FEM_My_partition();
00423   
00424   /*Add a refinement object to FEM array*/
00425   CkPrintf("[%d] begin init\n",myChunk);
00426   FEM_REFINE2D_Init();
00427   CkPrintf("[%d] end init\n",myChunk);
00428   
00429   myGlobals g;
00430   FEM_Register(&g,(FEM_PupFn)pup_myGlobals);
00431   init_myGlobal(&g);
00432   
00433   g.nnodes = FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_NODE);
00434   int maxNodes = g.nnodes;
00435   g.maxnodes=2*maxNodes;
00436   g.m_i_fid=FEM_Create_field(FEM_DOUBLE,1,0,sizeof(double));
00437   resize_nodes((void *)&g,&g.nnodes,&maxNodes);
00438   int nghost=0;
00439   g.nelems=FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_ELEM);
00440   g.maxelems=g.nelems;
00441   resize_elems((void *)&g,&g.nelems,&g.maxelems);
00442 
00443   FEM_REFINE2D_Newmesh(FEM_Mesh_default_read(),FEM_NODE,FEM_ELEM);
00444   
00445   //Initialize associated data
00446   for (i=0;i<g.maxnodes;i++) {
00447     g.R_net[i]=g.d[i]=g.v[i]=g.a[i]=vector2d(0.0);
00448   }
00449   
00450   //Apply a small initial perturbation to positions
00451   for (i=0;i<g.nnodes;i++) {
00452     const double max=1.0e-15/15.0; //Tiny perturbation
00453     g.d[i].x+=max*(i&15);
00454     g.d[i].y+=max*((i+5)&15);
00455   }
00456   
00457   int fid=FEM_Create_field(FEM_DOUBLE,2,0,sizeof(vector2d));
00458   
00459   for (i=0;i<g.nelems;i++){
00460     checkTriangle(g,i);
00461   } 
00462   sleep(5);
00463   //Timeloop
00464   if (CkMyPe()==0){
00465     CkPrintf("Entering timeloop\n");
00466   } 
00467   //  int tSteps=0x70FF00FF;
00468   int tSteps=10;
00469   int z=13;
00470   calcMasses(g);
00471   double startTime=CkWallTimer();
00472   double curArea=2.5e-5/1024;
00473   int t = 0;
00474   if (1) { //Publish data to the net
00475     NetFEM n=NetFEM_Begin(myChunk,t,2,NetFEM_WRITE);
00476     int count=0;
00477     double *vcoord = new double[2*g.nnodes];
00478     double *vnodeid = new double[g.nnodes];
00479     int *maptovalid = new int[g.nnodes];
00480     for(int i=0;i<g.nnodes;i++){
00481       if(g.validNode[i]){
00482     vcoord[2*count] = ((double *)g.coord)[2*i];
00483     vcoord[2*count+1] = ((double *)g.coord)[2*i+1];
00484     maptovalid[i] = count;
00485     printf("~~~~~~~ %d %d %.6lf %.6lf \n",count,i,vcoord[2*count],vcoord[2*count+1]);
00486     vnodeid[count] = i;
00487     count++;    
00488       }
00489     }
00490     NetFEM_Nodes(n,count,(double *)vcoord,"Position (m)");
00491     NetFEM_Scalar(n,vnodeid,1,"Node ID");
00492     /*    NetFEM_Vector(n,(double *)g.d,"Displacement (m)");
00493       NetFEM_Vector(n,(double *)g.v,"Velocity (m/s)");*/
00494     count=0;
00495     int *vconn = new int[3*g.nelems];
00496     double *vid = new double[3*g.nelems];
00497     for(int i=0;i<g.nelems;i++){
00498       if(g.validElem[i]){
00499     vconn[3*count] = maptovalid[g.conn[3*i]];
00500     vconn[3*count+1] = maptovalid[g.conn[3*i+1]];
00501     vconn[3*count+2] = maptovalid[g.conn[3*i+2]];
00502     printf("~~~~~~~ %d %d < %d,%d %d,%d %d,%d >\n",count,i,vconn[3*count],g.conn[3*i],vconn[3*count+1],g.conn[3*i+1],vconn[3*count+2],g.conn[3*i+2]);
00503     vid[count]=count;
00504     count++;    
00505       }
00506     }
00507     NetFEM_Elements(n,count,3,(int *)vconn,"Triangles");
00508     NetFEM_Scalar(n,vid,1,"Element ID");
00509     /*  NetFEM_Scalar(n,g.S22,1,"Y Stress (pure)");
00510     NetFEM_Scalar(n,g.S12,1,"Shear Stress (pure)");*/
00511     NetFEM_End(n);
00512     delete [] vcoord;
00513     delete [] vconn;
00514     delete [] maptovalid;
00515     delete [] vid;
00516     delete [] vnodeid;
00517   }
00518   for (t=1;t<=tSteps;t++) {
00519     /*    if (1) { //Structural mechanics
00520     //Compute forces on nodes exerted by elements
00521     CST_NL(g.coord,g.conn,g.R_net,g.d,matConst,g.nnodes,g.nelems,g.S11,g.S22,g.S12);
00522     //Communicate net force on shared nodes
00523     FEM_Update_field(fid,g.R_net);
00524     //Advance node positions
00525     advanceNodes(dt,g.nnodes,g.coord,g.R_net,g.a,g.v,g.d,g.m_i,(t%4)==0);
00526     }*/
00527     
00528     //Debugging/perf. output
00529     double curTime=CkWallTimer();
00530     double total=curTime-startTime;
00531     startTime=curTime;
00532     /*    if (CkMyPe()==0 && (t%64==0))
00533       CkPrintf("%d %.6f sec for loop %d \n",CkNumPes(),total,t);*/
00534     /*   if (0 && t%16==0) {
00535      CkPrintf("    Triangle 0:\n");
00536      for (int j=0;j<3;j++) {
00537      int n=g.conn[0][j];
00538      CkPrintf("    Node %d: coord=(%.4f,%.4f)  d=(%.4g,%.4g)\n",
00539      n,g.coord[n].x,g.coord[n].y,g.d[n].x,g.d[n].y);
00540      }
00541      }*/
00542     //    if (t%512==0)
00543     //      FEM_Migrate();
00544     
00545     double *areas=new double[g.nelems];
00546     for (i=0;i<g.nelems;i++) {
00547       areas[i]=1.5*calcArea(g,i);
00548 //      areas[i]=calcArea(g,i);
00549     }
00550     
00551     //coarsen all steps
00552 //    areas[z] *= 2.0;
00553     z += 3;
00554     CkPrintf("[%d] Starting coarsening step: %d nodes, %d elements to %.3g\n", myChunk,g.nnodes,g.nelems,curArea);
00555     FEM_REFINE2D_Coarsen(FEM_Mesh_default_read(),FEM_NODE,(double *)g.coord,FEM_ELEM,areas);
00556     repeat_after_split((void *)&g);
00557     g.nelems = FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_ELEM);
00558     g.nnodes = FEM_Mesh_get_length(FEM_Mesh_default_read(),FEM_NODE);
00559     CkPrintf("[%d] Done with coarsening step: %d nodes, %d elements\n",
00560          myChunk,g.nnodes,g.nelems);
00561     if (1) { //Publish data to the net
00562       NetFEM n=NetFEM_Begin(myChunk,t,2,NetFEM_WRITE);
00563       int count=0;
00564       double *vcoord = new double[2*g.nnodes];
00565       double *vnodeid = new double[g.nnodes];
00566       int *maptovalid = new int[g.nnodes];
00567       for(int i=0;i<g.nnodes;i++){
00568     maptovalid[i] = -1;
00569     if(g.validNode[i]){
00570       vcoord[2*count] = ((double *)g.coord)[2*i];
00571       vcoord[2*count+1] = ((double *)g.coord)[2*i+1];
00572       maptovalid[i] = count;
00573       printf("node~~~~~~~ %d %d %.6lf %.6lf \n",count,i,vcoord[2*count],vcoord[2*count+1]);
00574       vnodeid[count] = i;
00575       count++;  
00576     }
00577       }
00578       NetFEM_Nodes(n,count,(double *)vcoord,"Position (m)");
00579       NetFEM_Scalar(n,vnodeid,1,"Node ID");
00580       /*    NetFEM_Vector(n,(double *)g.d,"Displacement (m)");
00581         NetFEM_Vector(n,(double *)g.v,"Velocity (m/s)");*/
00582       count=0;
00583       int *vconn = new int[3*g.nelems];
00584       double *vid = new double[3*g.nelems];
00585       for(int i=0;i<g.nelems;i++){
00586     if(g.validElem[i]){
00587       vconn[3*count] = maptovalid[g.conn[3*i]];
00588       vconn[3*count+1] = maptovalid[g.conn[3*i+1]];
00589       vconn[3*count+2] = maptovalid[g.conn[3*i+2]];
00590       printf("~~~~~~~ %d %d < %d,%d %d,%d %d,%d >\n",count,i,vconn[3*count],g.conn[3*i],vconn[3*count+1],g.conn[3*i+1],vconn[3*count+2],g.conn[3*i+2]);
00591       vid[count]=count;
00592       count++;  
00593     }
00594       }
00595       NetFEM_Elements(n,count,3,(int *)vconn,"Triangles");
00596       NetFEM_Scalar(n,vid,1,"Element ID");
00597       /*    NetFEM_Scalar(n,g.S22,1,"Y Stress (pure)");
00598         NetFEM_Scalar(n,g.S12,1,"Shear Stress (pure)");*/
00599       NetFEM_End(n);
00600       delete [] vcoord;
00601       delete [] vconn;
00602       delete [] maptovalid;
00603       delete [] vid;
00604       delete [] vnodeid;
00605     }
00606   }
00607   if (CkMyPe()==0)
00608     CkPrintf("Driver finished\n");
00609 }
00610 
00611
libs/ck-libs/TMRC2D/pgm-coarsen.C
