← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Plugin/Extract.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme, Christophe Prud'homme
  • Date: 2009-09-02 18:12:15 UTC
  • mfrom: (1.2.8 upstream)
  • Revision ID: james.westby@ubuntu.com-20090902181215-yla8zvcas2ucvkm9
[Christophe Prud'homme]
* New upstream release
  + fixed surface mesh orientation bug introduced in 2.4.0;
  + mesh and graphics code refactoring;
  + small usability enhancements and bug fixes.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Plugin/Extract.cpp
4
4
// bugs and problems to <gmsh@geuz.org>.
5
5
 
6
6
#include "GmshConfig.h"
 
7
#include "GmshDefines.h"
7
8
#include "Extract.h"
8
 
 
9
 
#if defined(HAVE_MATH_EVAL)
10
 
#include "matheval.h"
11
 
#endif
 
9
#include "mathEvaluator.h"
12
10
 
13
11
StringXNumber ExtractOptions_Number[] = {
14
12
  {GMSH_FULLRC, "TimeStep", NULL, -1.},
16
14
};
17
15
 
18
16
StringXString ExtractOptions_String[] = {
19
 
#if defined(HAVE_MATH_EVAL)
20
17
  {GMSH_FULLRC, "Expression0", NULL, "Sqrt(v0^2+v1^2+v2^2)"},
21
 
#else
22
 
  {GMSH_FULLRC, "Expression0", NULL, "v0"},
23
 
#endif
24
18
  {GMSH_FULLRC, "Expression1", NULL, ""},
25
19
  {GMSH_FULLRC, "Expression2", NULL, ""},
26
20
  {GMSH_FULLRC, "Expression3", NULL, ""},
39
33
  }
40
34
}
41
35
 
42
 
void GMSH_ExtractPlugin::getName(char *name) const
43
 
{
44
 
  strcpy(name, "Extract");
45
 
}
46
 
 
47
 
void GMSH_ExtractPlugin::getInfos(char *author, char *copyright, char *help_text) const
48
 
{
49
 
  strcpy(author, "C. Geuzaine, J.-F. Remacle");
50
 
  strcpy(copyright, "C. Geuzaine, J.-F. Remacle");
51
 
  strcpy(help_text,
52
 
         "Plugin(Extract) extracts a combination of\n"
 
36
std::string GMSH_ExtractPlugin::getHelp() const
 
37
{
 
38
  return "Plugin(Extract) extracts a combination of\n"
53
39
         "components from the `TimeStep'th time step\n"
54
40
         "in the view `iView'. If only `Expression0' is\n"
55
41
         "given (and `Expression1', ..., `Expression8' are\n"
69
55
         "in the view. If `iView' < 0, the plugin is run on\n"
70
56
         "the current view.\n"
71
57
         "\n"
72
 
         "Plugin(Extract) creates one new view.\n");
 
58
         "Plugin(Extract) creates one new view.\n";
73
59
}
74
60
 
75
61
int GMSH_ExtractPlugin::getNbOptions() const
92
78
  return &ExtractOptions_String[iopt];
93
79
}
94
80
 
95
 
void GMSH_ExtractPlugin::catchErrorMessage(char *errorMessage) const
 
81
static std::vector<double> *incrementList(PViewDataList *data, int numComp, 
 
82
                                          int type)
96
83
{
97
 
  strcpy(errorMessage, "Extract failed...");
 
84
  switch(type){
 
85
  case TYPE_PNT:
 
86
    if     (numComp == 1){ data->NbSP++; return &data->SP; }
 
87
    else if(numComp == 3){ data->NbVP++; return &data->VP; }
 
88
    else if(numComp == 9){ data->NbTP++; return &data->TP; }
 
89
    break;
 
90
  case TYPE_LIN:
 
91
    if     (numComp == 1){ data->NbSL++; return &data->SL; }
 
92
    else if(numComp == 3){ data->NbVL++; return &data->VL; }
 
93
    else if(numComp == 9){ data->NbTL++; return &data->TL; }
 
94
    break;
 
95
  case TYPE_TRI:
 
96
    if     (numComp == 1){ data->NbST++; return &data->ST; }
 
97
    else if(numComp == 3){ data->NbVT++; return &data->VT; }
 
98
    else if(numComp == 9){ data->NbTT++; return &data->TT; }
 
99
    break;
 
100
  case TYPE_QUA:
 
101
    if     (numComp == 1){ data->NbSQ++; return &data->SQ; }
 
102
    else if(numComp == 3){ data->NbVQ++; return &data->VQ; }
 
103
    else if(numComp == 9){ data->NbTQ++; return &data->TQ; }
 
104
    break;
 
105
  case TYPE_TET:
 
106
    if     (numComp == 1){ data->NbSS++; return &data->SS; }
 
107
    else if(numComp == 3){ data->NbVS++; return &data->VS; }
 
108
    else if(numComp == 9){ data->NbTS++; return &data->TS; }
 
109
    break;
 
110
  case TYPE_HEX:
 
111
    if     (numComp == 1){ data->NbSH++; return &data->SH; }
 
112
    else if(numComp == 3){ data->NbVH++; return &data->VH; }
 
113
    else if(numComp == 9){ data->NbTH++; return &data->TH; }
 
114
    break;
 
115
  case TYPE_PRI:
 
116
    if     (numComp == 1){ data->NbSI++; return &data->SI; }
 
117
    else if(numComp == 3){ data->NbVI++; return &data->VI; }
 
118
    else if(numComp == 9){ data->NbTI++; return &data->TI; }
 
119
    break;
 
120
  case TYPE_PYR:
 
121
    if     (numComp == 1){ data->NbSY++; return &data->SY; }
 
122
    else if(numComp == 3){ data->NbVY++; return &data->VY; }
 
123
    else if(numComp == 9){ data->NbTY++; return &data->TY; }
 
124
    break;
 
125
  }
 
126
  return 0;
98
127
}
99
128
 
100
 
static void extract(const char *expr[9], List_T *inList, int inNb, 
101
 
                    List_T *outListScalar, int *outNbScalar, 
102
 
                    List_T *outListVector, int *outNbVector, 
103
 
                    List_T *outListTensor, int *outNbTensor, 
104
 
                    int timeStep, int nbNod, int nbComp)
 
129
PView *GMSH_ExtractPlugin::execute(PView *view)
105
130
{
106
 
  if(!inNb)
107
 
    return;
108
 
 
109
 
  int outNbComp, *outNb;
110
 
  List_T *outList;
111
 
 
112
 
  if(strlen(expr[3]) || strlen(expr[4]) || strlen(expr[5]) || 
113
 
     strlen(expr[6]) || strlen(expr[7]) || strlen(expr[8])){
114
 
    outNbComp = 9;
115
 
    outNb = outNbTensor;
116
 
    outList = outListTensor;
 
131
  int timeStep = (int)ExtractOptions_Number[0].def;
 
132
  int iView = (int)ExtractOptions_Number[1].def;
 
133
  std::vector<std::string> expr(9);
 
134
  for(int i = 0; i < 9; i++) expr[i] = ExtractOptions_String[i].def;
 
135
  
 
136
  PView *v1 = getView(iView, view);
 
137
  if(!v1) return view;
 
138
  PViewData *data1 = v1->getData();
 
139
 
 
140
  if(data1->hasMultipleMeshes()){
 
141
    Msg::Error("Extract plugin cannot be applied to multi-mesh views");
 
142
    return view;
 
143
  }
 
144
 
 
145
  int numComp2;
 
146
  if(expr[3].size() || expr[4].size() || expr[5].size() || 
 
147
     expr[6].size() || expr[7].size() || expr[8].size()){
 
148
    numComp2 = 9;
117
149
    for(int i = 0; i < 9; i++)
118
 
      if(!strlen(expr[i])) expr[i] = "0";
 
150
      if(expr[i].empty()) expr[i] = "0";
119
151
  }
120
 
  else if(strlen(expr[1]) || strlen(expr[2])){
121
 
    outNbComp = 3;
122
 
    outNb = outNbVector;
123
 
    outList = outListVector;
 
152
  else if(expr[1].size() || expr[2].size()){
 
153
    numComp2 = 3;
124
154
    for(int i = 0; i < 3; i++)
125
 
      if(!strlen(expr[i])) expr[i] = "0";
 
155
      if(expr[i].empty()) expr[i] = "0";
126
156
  }
127
157
  else{
128
 
    outNbComp = 1;
129
 
    outNb = outNbScalar;
130
 
    outList = outListScalar;
131
 
  }
132
 
 
133
 
  // if we have MathEval, we can evaluate arbitrary expressions;
134
 
  // otherwise, we only allow to extract single components
135
 
 
136
 
#if defined(HAVE_MATH_EVAL)
137
 
  void *f[9] = { NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL };
138
 
  for(int i = 0; i < outNbComp; i++){
139
 
    f[i] = evaluator_create((char*)expr[i]);
140
 
    if(!f[i]){
141
 
      Msg::Error("Invalid expression '%s'", expr[i]);
142
 
      for(int j = 0; j < i; j++)
143
 
        if(f[j]) evaluator_destroy(f[j]);
144
 
      return;
145
 
    }
146
 
  }
147
 
#else
148
 
  int comp[9] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 };
149
 
  for(int i = 0; i < outNbComp; i++){
150
 
    if     (!strcmp(expr[i], "v0")) comp[i] = 0;
151
 
    else if(!strcmp(expr[i], "v1")) comp[i] = 1;
152
 
    else if(!strcmp(expr[i], "v2")) comp[i] = 2;
153
 
    else if(!strcmp(expr[i], "v3")) comp[i] = 3;
154
 
    else if(!strcmp(expr[i], "v4")) comp[i] = 4;
155
 
    else if(!strcmp(expr[i], "v5")) comp[i] = 5;
156
 
    else if(!strcmp(expr[i], "v6")) comp[i] = 6;
157
 
    else if(!strcmp(expr[i], "v7")) comp[i] = 7;
158
 
    else if(!strcmp(expr[i], "v8")) comp[i] = 8;
159
 
    else{
160
 
      Msg::Error("Invalid expression '%s'", expr[i]);
161
 
      return;
162
 
    }
163
 
  }
164
 
#endif
165
 
 
166
 
  int timeBeg = (timeStep < 0) ? 0 : timeStep;
167
 
  int timeEnd = (timeStep < 0) ? -timeStep : timeStep + 1;
168
 
 
169
 
  int nb = List_Nbr(inList) / inNb;
170
 
  for(int i = 0; i < List_Nbr(inList); i += nb) {
171
 
    double *x = (double *)List_Pointer_Fast(inList, i);
172
 
    double *y = (double *)List_Pointer_Fast(inList, i + nbNod);
173
 
    double *z = (double *)List_Pointer_Fast(inList, i + 2 * nbNod);
174
 
    for(int j = 0; j < 3 * nbNod; j++)
175
 
      List_Add(outList, List_Pointer_Fast(inList, i + j));
176
 
    for(int j = timeBeg; j < timeEnd; j++){
177
 
      for(int k = 0; k < nbNod; k++){
178
 
        double xx = x[k];
179
 
        double yy = y[k];
180
 
        double zz = z[k];
181
 
        double d[9] = {0., 0., 0., 0., 0., 0., 0., 0., 0.};
182
 
        for(int l = 0; l < nbComp; l++)
183
 
          List_Read(inList, i + 3 * nbNod + nbNod * nbComp * j + nbComp * k + l, &d[l]);
184
 
        for(int l = 0; l < outNbComp; l++){
185
 
#if defined(HAVE_MATH_EVAL)
186
 
          char *names[] = { "x", "y", "z", "v0", "v1", "v2", "v3", "v4", "v5", "v6", "v7", "v8" };
187
 
          double values[] = { xx, yy, zz, d[0], d[1], d[2], d[3], d[4], d[5], d[6], d[7], d[8] };
188
 
          double res = evaluator_evaluate(f[l], sizeof(names)/sizeof(names[0]), names, values);
189
 
#else
190
 
          double res = d[comp[l]];
191
 
#endif
192
 
          List_Add(outList, &res);
 
158
    numComp2 = 1;
 
159
  }
 
160
  expr.resize(numComp2);
 
161
 
 
162
  const char *names[] = 
 
163
    { "x", "y", "z", "v0", "v1", "v2", "v3", "v4", "v5", "v6", "v7", "v8" };
 
164
  unsigned int numVariables = sizeof(names) / sizeof(names[0]);
 
165
  std::vector<std::string> variables(numVariables);
 
166
  for(unsigned int i = 0; i < numVariables; i++) variables[i] = names[i];
 
167
  mathEvaluator f(expr, variables);
 
168
  if(expr.empty()) return view;
 
169
  std::vector<double> values(numVariables), res(numComp2);
 
170
          
 
171
  PView *v2 = new PView();
 
172
  PViewDataList *data2 = getDataList(v2);
 
173
 
 
174
  if(timeStep < 0){
 
175
    timeStep = - data1->getNumTimeSteps();
 
176
  }
 
177
  else if(timeStep > data1->getNumTimeSteps() - 1){
 
178
    Msg::Error("Invalid time step (%d) in View[%d]: using all steps instead",
 
179
               timeStep, v1->getIndex());
 
180
    timeStep = - data1->getNumTimeSteps();
 
181
  }
 
182
 
 
183
  for(int ent = 0; ent < data1->getNumEntities(0); ent++){
 
184
    for(int ele = 0; ele < data1->getNumElements(0, ent); ele++){
 
185
      if(data1->skipElement(0, ent, ele)) continue;
 
186
      int numNodes = data1->getNumNodes(0, ent, ele);
 
187
      int type = data1->getType(0, ent, ele);
 
188
      int numComp = data1->getNumComponents(0, ent, ele);
 
189
      std::vector<double> *out = incrementList(data2, numComp2, type);
 
190
      std::vector<double> x(numNodes), y(numNodes), z(numNodes);
 
191
      for(int nod = 0; nod < numNodes; nod++)
 
192
        data1->getNode(0, ent, ele, nod, x[nod], y[nod], z[nod]);
 
193
      for(int nod = 0; nod < numNodes; nod++) out->push_back(x[nod]); 
 
194
      for(int nod = 0; nod < numNodes; nod++) out->push_back(y[nod]); 
 
195
      for(int nod = 0; nod < numNodes; nod++) out->push_back(z[nod]); 
 
196
      int timeBeg = (timeStep < 0) ? 0 : timeStep;
 
197
      int timeEnd = (timeStep < 0) ? -timeStep : timeStep + 1;
 
198
      for(int step = timeBeg; step < timeEnd; step++){
 
199
        for(int nod = 0; nod < numNodes; nod++){
 
200
          std::vector<double> v(std::max(9, numComp), 0.);
 
201
          for(int comp = 0; comp < numComp; comp++)
 
202
            data1->getValue(step, ent, ele, nod, comp, v[comp]);
 
203
          values[0] = x[nod]; values[1] = y[nod]; values[2] = z[nod];
 
204
          for(int i = 0; i < 9; i++) values[3 + i] = v[i];
 
205
          if(f.eval(values, res))
 
206
            for(int i = 0; i < numComp2; i++)
 
207
              out->push_back(res[i]);
193
208
        }
194
209
      }
195
210
    }
196
 
    (*outNb)++;
197
 
  }
198
 
 
199
 
#if defined(HAVE_MATH_EVAL)
200
 
  for(int i = 0; i < outNbComp; i++)
201
 
    evaluator_destroy(f[i]);
202
 
#endif
203
 
}
204
 
 
205
 
PView *GMSH_ExtractPlugin::execute(PView *v)
206
 
{
207
 
  int step = (int)ExtractOptions_Number[0].def;
208
 
  int iView = (int)ExtractOptions_Number[1].def;
209
 
  const char *expr[9] = { ExtractOptions_String[0].def.c_str(), 
210
 
                          ExtractOptions_String[1].def.c_str(),
211
 
                          ExtractOptions_String[2].def.c_str(),
212
 
                          ExtractOptions_String[3].def.c_str(),
213
 
                          ExtractOptions_String[4].def.c_str(),
214
 
                          ExtractOptions_String[5].def.c_str(),
215
 
                          ExtractOptions_String[6].def.c_str(),
216
 
                          ExtractOptions_String[7].def.c_str(),
217
 
                          ExtractOptions_String[8].def.c_str() };
218
 
 
219
 
  PView *v1 = getView(iView, v);
220
 
  if(!v1) return v;
221
 
 
222
 
  PViewDataList *data1 = getDataList(v1);
223
 
  if(!data1) return v;
224
 
 
225
 
  PView *v2 = new PView(true, data1->getNumElements());
226
 
 
227
 
  PViewDataList *data2 = getDataList(v2);
228
 
  if(!data2) return v;
229
 
 
230
 
  if(step < 0){
231
 
    step = - data1->getNumTimeSteps();
232
 
  }
233
 
  else if(step > data1->getNumTimeSteps() - 1){
234
 
    Msg::Error("Invalid time step (%d) in View[%d]: using all steps instead",
235
 
        step, v1->getIndex());
236
 
    step = - data1->getNumTimeSteps();
237
 
  }
238
 
 
239
 
  // points
240
 
  extract(expr, data1->SP, data1->NbSP, data2->SP, &data2->NbSP, 
241
 
          data2->VP, &data2->NbVP, data2->TP, &data2->NbTP, step, 1, 1);
242
 
  extract(expr, data1->VP, data1->NbVP, data2->SP, &data2->NbSP,
243
 
          data2->VP, &data2->NbVP, data2->TP, &data2->NbTP, step, 1, 3);
244
 
  extract(expr, data1->TP, data1->NbTP, data2->SP, &data2->NbSP,
245
 
          data2->VP, &data2->NbVP, data2->TP, &data2->NbTP, step, 1, 9);
246
 
  // lines                                                                              
247
 
  extract(expr, data1->SL, data1->NbSL, data2->SL, &data2->NbSL,
248
 
          data2->VL, &data2->NbVL, data2->TL, &data2->NbTL, step, 2, 1);
249
 
  extract(expr, data1->VL, data1->NbVL, data2->SL, &data2->NbSL,
250
 
          data2->VL, &data2->NbVL, data2->TL, &data2->NbTL, step, 2, 3);
251
 
  extract(expr, data1->TL, data1->NbTL, data2->SL, &data2->NbSL,
252
 
          data2->VL, &data2->NbVL, data2->TL, &data2->NbTL, step, 2, 9);
253
 
  // triangles                                                                          
254
 
  extract(expr, data1->ST, data1->NbST, data2->ST, &data2->NbST,
255
 
          data2->VT, &data2->NbVT, data2->TT, &data2->NbTT, step, 3, 1);
256
 
  extract(expr, data1->VT, data1->NbVT, data2->ST, &data2->NbST,
257
 
          data2->VT, &data2->NbVT, data2->TT, &data2->NbTT, step, 3, 3);
258
 
  extract(expr, data1->TT, data1->NbTT, data2->ST, &data2->NbST,
259
 
          data2->VT, &data2->NbVT, data2->TT, &data2->NbTT, step, 3, 9);
260
 
  // quadrangles                                                                        
261
 
  extract(expr, data1->SQ, data1->NbSQ, data2->SQ, &data2->NbSQ,
262
 
          data2->VQ, &data2->NbVQ, data2->TQ, &data2->NbTQ, step, 4, 1);
263
 
  extract(expr, data1->VQ, data1->NbVQ, data2->SQ, &data2->NbSQ,
264
 
          data2->VQ, &data2->NbVQ, data2->TQ, &data2->NbTQ, step, 4, 3);
265
 
  extract(expr, data1->TQ, data1->NbTQ, data2->SQ, &data2->NbSQ,
266
 
          data2->VQ, &data2->NbVQ, data2->TQ, &data2->NbTQ, step, 4, 9);
267
 
  // tets                                                                               
268
 
  extract(expr, data1->SS, data1->NbSS, data2->SS, &data2->NbSS,
269
 
          data2->VS, &data2->NbVS, data2->TS, &data2->NbTS, step, 4, 1);
270
 
  extract(expr, data1->VS, data1->NbVS, data2->SS, &data2->NbSS,
271
 
          data2->VS, &data2->NbVS, data2->TS, &data2->NbTS, step, 4, 3);
272
 
  extract(expr, data1->TS, data1->NbTS, data2->SS, &data2->NbSS,
273
 
          data2->VS, &data2->NbVS, data2->TS, &data2->NbTS, step, 4, 9);
274
 
  // hexas                                                                              
275
 
  extract(expr, data1->SH, data1->NbSH, data2->SH, &data2->NbSH,
276
 
          data2->VH, &data2->NbVH, data2->TH, &data2->NbTH, step, 8, 1);
277
 
  extract(expr, data1->VH, data1->NbVH, data2->SH, &data2->NbSH,
278
 
          data2->VH, &data2->NbVH, data2->TH, &data2->NbTH, step, 8, 3);
279
 
  extract(expr, data1->TH, data1->NbTH, data2->SH, &data2->NbSH,
280
 
          data2->VH, &data2->NbVH, data2->TH, &data2->NbTH, step, 8, 9);
281
 
  // prisms                                                                             
282
 
  extract(expr, data1->SI, data1->NbSI, data2->SI, &data2->NbSI,
283
 
          data2->VI, &data2->NbVI, data2->TI, &data2->NbTI, step, 6, 1);
284
 
  extract(expr, data1->VI, data1->NbVI, data2->SI, &data2->NbSI,
285
 
          data2->VI, &data2->NbVI, data2->TI, &data2->NbTI, step, 6, 3);
286
 
  extract(expr, data1->TI, data1->NbTI, data2->SI, &data2->NbSI,
287
 
          data2->VI, &data2->NbVI, data2->TI, &data2->NbTI, step, 6, 9);
288
 
  // pyramids                                                                           
289
 
  extract(expr, data1->SY, data1->NbSY, data2->SY, &data2->NbSY,
290
 
          data2->VY, &data2->NbVY, data2->TY, &data2->NbTY, step, 5, 1);
291
 
  extract(expr, data1->VY, data1->NbVY, data2->SY, &data2->NbSY,
292
 
          data2->VY, &data2->NbVY, data2->TY, &data2->NbTY, step, 5, 3);
293
 
  extract(expr, data1->TY, data1->NbTY, data2->SY, &data2->NbSY,
294
 
          data2->VY, &data2->NbVY, data2->TY, &data2->NbTY, step, 5, 9);
295
 
 
296
 
  if(step < 0)
297
 
    for(int i = 0; i < List_Nbr(data1->Time); i++)
298
 
      List_Add(data2->Time, List_Pointer(data1->Time, i));
 
211
  }
 
212
  
 
213
  if(timeStep < 0){
 
214
    for(int i = 0; i < data1->getNumTimeSteps(); i++)
 
215
      data2->Time.push_back(data1->getTime(i));
 
216
  }
299
217
  else
300
 
    List_Add(data2->Time, List_Pointer(data1->Time, step));
301
 
 
 
218
    data2->Time.push_back(data1->getTime(timeStep));
 
219
        
302
220
  data2->setName(data1->getName() + "_Extract");
303
221
  data2->setFileName(data1->getName() + "_Extract.pos");
304
222
  data2->finalize();