← Back to branch summary

~ubuntu-branches/debian/squeeze/gmsh/squeeze

« back to all changes in this revision

Viewing changes to Solver/elasticitySolver.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme
  • Date: 2009-09-27 17:36:40 UTC
  • mfrom: (1.4.2 upstream)
  • mto: This revision was merged to the branch mainline in revision 10.
  • Revision ID: james.westby@ubuntu.com-20090927173640-meoeywl4f5hq5qas
Tags: 2.4.2.dfsg-1
[Christophe Prud'homme]
* New upstream release
  + solver code refactoring
  + better IDE integration.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Solver/elasticitySolver.cpp
 
1
// Gmsh - Copyright (C) 1997-2009 C. Geuzaine, J.-F. Remacle
 
2
//
 
3
// See the LICENSE.txt file for license information. Please report all
 
4
// bugs and problems to <gmsh@geuz.org>.
 
5
 
 
6
#include <string.h>
 
7
#include "elasticityTerm.h"
 
8
#include "MTriangle.h"
 
9
#include "MTetrahedron.h"
 
10
#include "elasticitySolver.h"
 
11
#include "linearSystemTAUCS.h"
 
12
#include "Numeric.h"
 
13
 
 
14
static double vonMises(dofManager<double,double> *a, MElement *e, 
 
15
                       double u, double v, double w, 
 
16
                       double E, double nu, int _tag)
 
17
{
 
18
  double valx[256];
 
19
  double valy[256];
 
20
  double valz[256];
 
21
  for (int k = 0; k < e->getNumVertices(); k++){
 
22
    valx[k] = a->getDofValue(e->getVertex(k), 0, _tag);
 
23
    valy[k] = a->getDofValue(e->getVertex(k), 1, _tag);
 
24
    valz[k] = a->getDofValue(e->getVertex(k), 2, _tag);
 
25
  }
 
26
  double gradux[3];
 
27
  double graduy[3];
 
28
  double graduz[3];
 
29
  e->interpolateGrad(valx, u, v, w, gradux);
 
30
  e->interpolateGrad(valy, u, v, w, graduy);
 
31
  e->interpolateGrad(valz, u, v, w, graduz);
 
32
  
 
33
  double eps[6] = {gradux[0], graduy[1], graduz[2],
 
34
                   0.5 * (gradux[1] + graduy[0]),
 
35
                   0.5 * (gradux[2] + graduz[0]),
 
36
                   0.5 * (graduy[2] + graduz[1])};
 
37
  double A = E / (1. + nu);
 
38
  double B = A * (nu / (1. - 2 * nu));
 
39
  double trace = eps[0] + eps[1] + eps[2] ;
 
40
  double sxx = A * eps[0] + B * trace;
 
41
  double syy = A * eps[1] + B * trace;
 
42
  double szz = A * eps[2] + B * trace;
 
43
  double sxy = A * eps[3];
 
44
  double sxz = A * eps[4];
 
45
  double syz = A * eps[5];
 
46
 
 
47
  double s[9] = {sxx, sxy, sxz,
 
48
                 sxy, syy, syz,
 
49
                 sxz, syz, szz};
 
50
  
 
51
  return ComputeVonMises(s);
 
52
}
 
53
  
 
54
void elasticitySolver::setMesh(const std::string &meshFileName)
 
55
{
 
56
  pModel = new GModel();
 
57
  pModel->readMSH(meshFileName.c_str());
 
58
  _dim = pModel->getNumRegions() ? 3 : 2;  
 
59
}
 
60
 
 
61
void elasticitySolver::readInputFile(const std::string &fn)
 
62
{
 
63
  FILE *f = fopen(fn.c_str(), "r");
 
64
  char what[256];
 
65
  while(!feof(f)){
 
66
    if(fscanf(f, "%s", what) != 1) return;
 
67
    // printf("%s\n", what);
 
68
    if (!strcmp(what,"ElasticMaterial")){
 
69
      double E, nu;
 
70
      int physical;
 
71
      if(fscanf(f, "%d %lf %lf", &physical, &E, &nu) != 3) return;
 
72
      elasticConstants[physical] = std::make_pair(E, nu);    
 
73
    }
 
74
    else if (!strcmp(what, "NodalDisplacement")){
 
75
      double val;
 
76
      int node, comp;
 
77
      if(fscanf(f, "%d %d %lf", &node, &comp, &val) != 3) return;
 
78
      nodalDisplacements[ std::make_pair(node, comp) ] = val;    
 
79
    }
 
80
    else if (!strcmp(what, "EdgeDisplacement")){
 
81
      double val;
 
82
      int edge, comp;
 
83
      if(fscanf(f, "%d %d %lf", &edge, &comp, &val) != 3) return;
 
84
      edgeDisplacements[ std::make_pair(edge, comp) ] = val;    
 
85
    }
 
86
    else if (!strcmp(what, "FaceDisplacement")){
 
87
      double val;
 
88
      int face, comp;
 
89
      if(fscanf(f, "%d %d %lf", &face, &comp, &val) != 3) return;
 
90
      faceDisplacements[ std::make_pair(face, comp) ] = val;    
 
91
    }
 
92
    else if (!strcmp(what, "NodalForce")){
 
93
      double val1, val2, val3;
 
94
      int node;
 
95
      if(fscanf(f, "%d %lf %lf %lf", &node, &val1, &val2, &val3) != 4) return;
 
96
      nodalForces[node] = SVector3(val1, val2, val3);    
 
97
    }
 
98
    else if (!strcmp(what, "LineForce")){
 
99
      double val1, val2, val3;
 
100
      int node;
 
101
      if(fscanf(f, "%d %lf %lf %lf", &node, &val1, &val2, &val3) != 4) return;
 
102
      //printf("%d %lf %lf %lf\n", node, val1, val2, val3);
 
103
      lineForces[node] = SVector3(val1, val2, val3);    
 
104
    }
 
105
    else if (!strcmp(what, "FaceForce")){
 
106
      double val1, val2, val3;
 
107
      int face;
 
108
      if(fscanf(f, "%d %lf %lf %lf", &face, &val1, &val2, &val3) != 4) return;
 
109
      faceForces[face] = SVector3(val1, val2, val3);    
 
110
    }
 
111
    else if (!strcmp(what, "VolumeForce")){
 
112
      double val1, val2, val3;
 
113
      int volume;
 
114
      if(fscanf(f, "%d %lf %lf %lf", &volume, &val1, &val2, &val3) != 4) return;
 
115
      volumeForces[volume] = SVector3(val1, val2, val3);    
 
116
    }
 
117
    else if (!strcmp(what, "MeshFile")){
 
118
      char name[245];
 
119
      if(fscanf(f, "%s", name) != 1) return;
 
120
      setMesh(name);
 
121
    }
 
122
  }
 
123
  fclose(f);
 
124
}
 
125
  
 
126
void elasticitySolver::solve()
 
127
{
 
128
#ifdef HAVE_TAUCS
 
129
  linearSystemCSRTaucs<double> *lsys = new linearSystemCSRTaucs<double>;
 
130
#else
 
131
  linearSystemCSRGmm<double> *lsys = new linearSystemCSRGmm<double>;
 
132
#endif
 
133
  pAssembler = new dofManager<double, double>(lsys);
 
134
  
 
135
  std::map<int, std::vector<GEntity*> > groups[4];
 
136
  pModel->getPhysicalGroups(groups);
 
137
  
 
138
  // we first do all fixations. the behavior of the dofManager is to 
 
139
  // give priority to fixations : when a dof is fixed, it cannot be
 
140
  // numbered afterwards
 
141
  
 
142
  for (std::map<std::pair<int, int>, double>::iterator it = nodalDisplacements.begin();
 
143
       it != nodalDisplacements.end(); ++it){
 
144
    MVertex *v = pModel->getMeshVertexByTag(it->first.first);
 
145
    pAssembler->fixVertex(v , it->first.second, _tag, it->second);    
 
146
    printf("-- Fixing node %3d comp %1d to %8.5f\n", 
 
147
           it->first.first, it->first.second, it->second);
 
148
  }
 
149
  
 
150
  for (std::map<std::pair<int, int>, double>::iterator it = edgeDisplacements.begin();
 
151
       it != edgeDisplacements.end(); ++it){
 
152
    elasticityTerm El(pModel, 1, 1, _tag);
 
153
    El.dirichletNodalBC(it->first.first, 1, it->first.second, _tag, 
 
154
                        simpleFunction<double>(it->second), *pAssembler);
 
155
    printf("-- Fixing edge %3d comp %1d to %8.5f\n", 
 
156
           it->first.first, it->first.second, it->second);
 
157
  }
 
158
  
 
159
  for (std::map<std::pair<int, int>, double>::iterator it = faceDisplacements.begin();
 
160
       it != faceDisplacements.end(); ++it){
 
161
    elasticityTerm El(pModel, 1, 1, _tag);
 
162
    El.dirichletNodalBC(it->first.first, 2, it->first.second, _tag, 
 
163
                        simpleFunction<double>(it->second), *pAssembler);
 
164
    printf("-- Fixing face %3d comp %1d to %8.5f\n", 
 
165
           it->first.first, it->first.second, it->second);
 
166
  }
 
167
  
 
168
  // we number the nodes : when a node is numbered, it cannot be numbered
 
169
  // again with another number. so we loop over elements
 
170
  for (std::map<int, std::pair<double, double> >::iterator it = elasticConstants.begin();
 
171
       it != elasticConstants.end() ; it++){
 
172
    int physical = it->first;
 
173
    std::vector<MVertex *> v;     
 
174
    pModel->getMeshVerticesForPhysicalGroup(_dim, physical, v);
 
175
    printf("Physical %d, dim: %d, nb vert: %d\n", physical, _dim, (int)v.size());
 
176
    for (unsigned int i = 0; i < v.size(); i++){  
 
177
      pAssembler->numberVertex(v[i], 0, _tag);  
 
178
      pAssembler->numberVertex(v[i], 1, _tag);  
 
179
      pAssembler->numberVertex(v[i], 2, _tag);  
 
180
    }
 
181
  }
 
182
  
 
183
  // Now we start the assembly process
 
184
  // First build the force vector
 
185
  // Nodes at geometrical vertices
 
186
  for (std::map<int, SVector3 >::iterator it = nodalForces.begin();
 
187
       it != nodalForces.end(); ++it){
 
188
    int iVertex = it->first;
 
189
    SVector3 f = it->second;
 
190
    std::vector<GEntity*> ent = groups[0][iVertex];
 
191
    for (unsigned int i = 0; i < ent.size(); i++){      
 
192
      pAssembler->assemble(ent[i]->mesh_vertices[0], 0, _tag, f.x());
 
193
      pAssembler->assemble(ent[i]->mesh_vertices[0], 1, _tag, f.y());
 
194
      pAssembler->assemble(ent[i]->mesh_vertices[0], 2, _tag, f.z());
 
195
      printf("-- Force on node %3d(%3d) : %8.5f %8.5f %8.5f\n", 
 
196
             ent[i]->mesh_vertices[0]->getNum(), iVertex, f.x(), f.y(), f.z());
 
197
    }
 
198
  }
 
199
 
 
200
  // line forces
 
201
  for (std::map<int, SVector3 >::iterator it = lineForces.begin();
 
202
       it != lineForces.end(); ++it){
 
203
    elasticityTerm El(pModel, 1, 1, _tag);
 
204
    int iEdge = it->first;
 
205
    SVector3 f = it->second;
 
206
    El.neumannNodalBC(iEdge, 1, 0, _tag, simpleFunction<double>(f.x()), *pAssembler);
 
207
    El.neumannNodalBC(iEdge, 1, 1, _tag, simpleFunction<double>(f.y()), *pAssembler);
 
208
    El.neumannNodalBC(iEdge, 1, 2, _tag, simpleFunction<double>(f.z()), *pAssembler);
 
209
    printf("-- Force on edge %3d : %8.5f %8.5f %8.5f\n", iEdge, f.x(), f.y(), f.z());
 
210
  }
 
211
 
 
212
  // face forces
 
213
  for (std::map<int, SVector3 >::iterator it = faceForces.begin();
 
214
       it != faceForces.end(); ++it){
 
215
    elasticityTerm El(pModel, 1, 1, _tag);
 
216
    int iFace = it->first;
 
217
    SVector3 f = it->second;
 
218
    El.neumannNodalBC(iFace, 2, 0, _tag, simpleFunction<double>(f.x()), *pAssembler);
 
219
    El.neumannNodalBC(iFace, 2, 1, _tag, simpleFunction<double>(f.y()), *pAssembler);
 
220
    El.neumannNodalBC(iFace, 2, 2, _tag, simpleFunction<double>(f.z()), *pAssembler);
 
221
    printf("-- Force on face %3d : %8.5f %8.5f %8.5f\n", iFace, f.x(), f.y(), f.z());
 
222
  }
 
223
  
 
224
  // volume forces
 
225
  for (std::map<int, SVector3 >::iterator it = volumeForces.begin();
 
226
       it != volumeForces.end(); ++it){
 
227
    elasticityTerm El(pModel, 1, 1, _tag);
 
228
    int iVolume = it->first;
 
229
    SVector3 f = it->second;
 
230
    El.setVector(f);
 
231
    printf("-- Force on volume %3d : %8.5f %8.5f %8.5f\n", iVolume, f.x(), f.y(), f.z());
 
232
    std::vector<GEntity*> ent = groups[_dim][iVolume];
 
233
    for (unsigned int i = 0; i < ent.size(); i++){
 
234
      El.addToRightHandSide(*pAssembler, ent[i]);
 
235
    }
 
236
  }
 
237
  
 
238
  // assembling the stifness matrix
 
239
  for (std::map<int, std::pair<double, double> >::iterator it = elasticConstants.begin();
 
240
       it != elasticConstants.end() ; it++){
 
241
    int physical = it->first;
 
242
    double E = it->second.first;
 
243
    double nu = it->second.second;
 
244
    elasticityTerm El(0, E, nu, _tag);
 
245
    std::vector<GEntity*> ent = groups[_dim][physical];
 
246
    for (unsigned int i = 0; i < ent.size(); i++){
 
247
      El.addToMatrix(*pAssembler, ent[i]);
 
248
    }
 
249
  }
 
250
 
 
251
  // solving
 
252
  lsys->systemSolve();
 
253
}