bbeb9bba47c43050c377a2a8b5cd9082266983c6
[charm.git] / examples / fem / crack2D / mesh.C
1 /**
2  * Read and maintain basic mesh for crack propagation code.
3  */
4
5 #include "converse.h"
6 #if CMK_STL_USE_DOT_H  /* Pre-standard C++ */
7 #  include <fstream.h>
8 #else /* ISO C++ */
9 #  include <fstream>
10    using namespace std;
11 #endif
12
13 #include <stddef.h>
14 #include "crack.h"
15
16 #define MAXLINE 1024
17 // Silly "skip n lines" macro.  Should actually only skip 
18 //  comment lines, but this will have to do...
19 #define cl(fd, buffer, n) do {\
20                                for(int _i=0; _i<(n); _i++) \
21                                  fd.getline(buffer, MAXLINE);\
22                              } while(0)
23
24 /// Read this mesh from this file:
25 void
26 readMesh(MeshData *mesh, const char *meshFile)
27 {
28   ifstream fm; // mesh file
29   char buf[MAXLINE];
30   int i, itmp;
31   double dtmp;
32   
33   fm.open(meshFile);
34   if (!fm)
35   {
36     crack_abort("Cannot open mesh file for reading.\n");
37     return;
38   }
39   
40   //Skip ramp-up for boundary conditions (already read by readConfig)
41   int numProp;
42   fm >> itmp >> itmp >> dtmp >> numProp;
43   for (i=0; i< numProp; i++) {
44     fm >> itmp >> dtmp;
45   }
46   
47   //read nodal co-ordinates
48   fm >> mesh->nn;
49   cl(fm,buf,1);
50   mesh->nodes=new Node[mesh->nn];
51   for(i=0; i<mesh->nn; i++)
52   {
53     Node *np = &(mesh->nodes[i]);
54     fm >> itmp >> np->pos.x >> np->pos.y;
55     np->isbnd = 0;
56   }
57   
58   //read nodal boundary conditions
59   fm >> mesh->numBound;
60   for (i=0; i<mesh->numBound; i++)
61   {
62     int j; // Boundary condition i applies to node j
63     fm >> j;
64     j--; //Fortran to C indexing
65     Node *np = &(mesh->nodes[j]);
66     np->isbnd = 1;
67     fm >> np->id1 >> np->id2 >> np->r.x >> np->r.y;
68   }
69   
70   //read cohesive elements
71   fm >> itmp >> mesh->nc >> itmp >> itmp >> itmp;
72   cl(fm,buf,1);
73   mesh->cohs=new Coh[mesh->nc];
74   for(i=0; i<mesh->nc; i++)
75   {
76     Coh *coh = &(mesh->cohs[i]);
77     fm >> coh->material; coh->material--;
78     int k;
79     for(k=0;k<6;k++) {
80       fm >> coh->conn[k]; coh->conn[k]--;
81     }
82     fm >> itmp; // What *is* this?
83   }
84   
85   // Read volumetric elements
86   fm >> itmp >> mesh->ne >> itmp;
87   cl(fm,buf,1);
88   mesh->vols=new Vol[mesh->ne];
89   for(i=0;i<mesh->ne;i++)
90   {
91     Vol *vol = &(mesh->vols[i]);
92     fm >> vol->material; vol->material--;
93     int k;
94     for(k=0;k<6;k++) {
95       fm >> vol->conn[k]; vol->conn[k]--;
96     }
97   }
98   fm.close();
99 }
100
101 void setupMesh(MeshData *mesh) {
102   int i;
103   //Zero out nodes
104   for(i=0; i<mesh->nn; i++)
105   {
106     Node *np = &(mesh->nodes[i]);
107     // np->pos is already set
108     np->xM.x = np->xM.y = 0;
109     np->disp.x = np->disp.y = 0.0;
110     np->vel.x = np->vel.y = 0.0;
111     np->accel.x = np->accel.y = 0.0;
112   }
113   
114   // Init cohesive elements, determine the length and angle
115   for(i=0; i<mesh->nc; i++)
116   {
117     Coh *coh = &(mesh->cohs[i]);
118     // coh->material and coh->conn[k] are already set
119     
120     Node *np1 = &(mesh->nodes[coh->conn[1]]);
121     Node *np2 = &(mesh->nodes[coh->conn[0]]);
122     coh->Sthresh[2] = coh->Sthresh[1] =
123       coh->Sthresh[0] = config.cohm[coh->material].Sinit;
124     double x = np1->pos.x - np2->pos.x;
125     double y = np1->pos.y - np2->pos.y;
126     coh->sidel[0] = sqrt(x*x+y*y);
127     coh->sidel[1] = x/coh->sidel[0];
128     coh->sidel[2] = y/coh->sidel[0];
129   }
130   
131   // Set up volumetric elements
132   for(i=0;i<mesh->ne;i++)
133   {
134     Vol *vol = &(mesh->vols[i]);
135     // vol->material and vol->conn[k] are already set
136     for(int k=0;k<3;k++)
137     {
138       vol->s11l[k] = 0.0;
139       vol->s22l[k] = 0.0;
140       vol->s12l[k] = 0.0;
141     }
142   }
143   
144   // Compute the mass of local elements:
145   for (i=0;i<mesh->ne;i++)
146   {
147     Vol *v=&(mesh->vols[i]);
148     Node *n[6];              //Pointers to each of the triangle's nodes
149     int k;                  //Loop index
150     for (k=0;k<6;k++)
151       n[k]=&(mesh->nodes[v->conn[k]]);
152     //Compute the mass of this element
153     double area=((n[1]->pos.x-n[0]->pos.x)*(n[2]->pos.y-n[0]->pos.y)-
154                  (n[2]->pos.x-n[0]->pos.x)*(n[1]->pos.y-n[0]->pos.y));
155     double mass=config.volm[v->material].rho*area/114.0;
156     //Divide the element's mass among the element's nodes
157     for (k=0;k<3;k++) {
158       n[k]->xM.x+=mass*3.0;
159       n[k]->xM.y+=mass*3.0;
160     }
161     for (k=3;k<6;k++) {
162       n[k]->xM.x+=mass*16.0;
163       n[k]->xM.y+=mass*16.0;
164     }
165   }  
166 }
167
168 void deleteMesh(MeshData *mesh)
169 {
170   delete[] mesh->nodes; mesh->nodes=NULL;
171   delete[] mesh->vols;  mesh->vols=NULL;
172   delete[] mesh->cohs;  mesh->cohs=NULL;
173 }
174