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