configure: Drop archaic and incompletely used test
[charm.git] / examples / fem / crack2D / config.C
1 /**
2  * Read configuration file for crack propagation code.
3  */
4
5 #include "converse.h"
6
7 #include <fstream>
8 using namespace std;
9
10 #include <stddef.h>
11 #include "crack.h"
12
13 /// Globally accessible (write-once) configuration data.
14 ConfigurationData config;
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 /// This routine reads the config data from the config file.
25 void readConfig(const char *configFile,const char *meshFile)
26 {
27   ifstream fg; // global parameter file
28   ifstream fm; // mesh file, for reading loading data
29   ConfigurationData *cfg=&config; //Always write to the global version
30   char buf[MAXLINE];
31   int i, itmp;
32   double dtmp;
33   
34   fg.open(configFile);
35   fm.open(meshFile);
36   if (!fg || !fm)
37   {
38     crack_abort("Cannot open configuration file for reading.\n");
39     return;
40   }
41   cl(fg,buf,3); // ignore first three lines
42   fg >> cfg->nTime >> cfg->steps >> cfg->tsintFull 
43      >> cfg->tsintEnergy >> cfg->restart;
44   cl(fg, buf, 2);
45   fg >> cfg->nplane >> cfg->ncoh >> cfg->lin;
46   cl(fg, buf, 2);
47   
48   //read volumetric material properties
49   fg >> cfg->numMatVol;
50   cfg->volm = new VolMaterial[cfg->numMatVol];
51   VolMaterial *v;
52   cl(fg, buf, 2);
53   for (i=0; i<cfg->numMatVol; i++)
54   {
55     v = &(cfg->volm[i]);
56     fg >> v->e1 >> v->e2 >> v->g12 >> v->g23;
57     fg >> v->xnu12 >> v->xnu23 >> v->rho;
58     fg >> v->alpha1 >> v->alpha2 >> v->theta;
59     cl(fg, buf, 2);
60   }
61   
62   //Compute the elastic stiffness constants for each material type
63   for (int matNo=0;matNo<cfg->numMatVol;matNo++)
64   {
65     double x, x1, x2, x3;
66     VolMaterial *vm=&(cfg->volm[matNo]);
67     switch (cfg->nplane)
68     {
69       case 1://Orthotropic plane strain
70         double sT,cT,xnu21;
71         cT = cos(vm->theta*1.74532925199e-2);
72         sT = sin(vm->theta*1.74532925199e-2);
73         xnu21 = vm->xnu12*vm->e2/vm->e1;
74         x = 1.0 - vm->xnu23*vm->xnu23 - 
75           2.0*vm->xnu12*xnu21*(1.0 + vm->xnu23);
76         x1 = vm->e1*(1.0-vm->xnu23*vm->xnu23) / x;
77         x2 = xnu21*vm->e1*(1.0+vm->xnu23) / x;
78         x3 = vm->e2*(vm->xnu23+vm->xnu12*xnu21) / x;
79         vm->c[2] = vm->e2*(1.0-vm->xnu12*xnu21) / x;
80         vm->c[0] = x1*cT*cT*cT*cT + 2.0*(x2+2.0*vm->g12)*cT*cT*sT*sT +
81           vm->c[2]*sT*sT*sT*sT;
82         vm->c[1] = x2*cT*cT + x3*sT*sT;
83         vm->c[3] = vm->g12*cT*cT + vm->g23*sT*sT;
84         break;
85       case 0: //Plane stress (isotropic)
86         vm->c[0] = vm->e1 / (1.0 - vm->xnu12*vm->xnu12);
87         vm->c[1] = vm->e1*vm->xnu12 / (1.0 - vm->xnu12*vm->xnu12);
88         vm->c[2] = vm->c[0];
89         vm->c[3] = vm->e1/ (2.0 * (1.0 + vm->xnu12));
90         break;
91       case 2: //Axisymmetric (isotropic)
92         vm->c[0] = vm->e1 * (1.0 - vm->xnu12) / ((1.0 + vm->xnu12)*
93                                                  (1.0 - 2.0*vm->xnu12));
94         vm->c[1] = vm->e1 * vm->xnu12 / ((1.0 + vm->xnu12)*
95                                          (1.0 - 2.0*vm->xnu12));
96         vm->c[2] = vm->e1 / (2.0*(1.0 + vm->xnu12));
97         break;
98       default:
99         crack_abort("Unknown planar analysis type in config file");
100     }
101   }
102
103   //read cohesive material properties
104   fg >> cfg->numMatCoh;
105   cfg->cohm = new CohMaterial[cfg->numMatCoh];
106   CohMaterial *c;
107   cl(fg, buf, 2);
108   for (i=0; i<cfg->numMatCoh; i++)
109   {
110     c = &(cfg->cohm[i]);
111     fg >> c->deltan >> c->deltat >> c->sigmax
112        >> c->taumax >> c->mu;
113     if (cfg->ncoh)
114       fg >> c->Sinit;
115     cl(fg, buf, 2);
116   }
117   
118   //read impact data
119   fg >> cfg->imp >> cfg->voImp >> cfg->massImp >> cfg->radiusImp;
120   cl(fg, buf, 2);
121   fg >> cfg->eImp >> cfg->xnuImp;
122   cl(fg, buf, 2);
123   cfg->voImp = 0; cfg->massImp = 1.0; cfg->radiusImp = 0; 
124   cfg->eImp = 1.0; cfg->xnuImp = 0.3;
125   
126   //read (& ignore) thermal load
127   fg >> dtmp;
128   fg.close();
129   
130   //read proportional ramp-up for boundary conditions
131   fm >> itmp >> itmp >> cfg->delta >> cfg->numProp;
132   cfg->delta /= (double) cfg->steps;
133   cfg->delta2 = cfg->delta*cfg->delta*0.5;
134   cfg->ts_proportion = new int[cfg->numProp];
135   cfg->proportion = new double[cfg->numProp];
136   for (i=0; i< cfg->numProp; i++) {
137     fm >> cfg->ts_proportion[i] >> cfg->proportion[i];
138   }
139   fm.close();
140 }
141
142
143