← Back to branch summary

~ubuntu-branches/ubuntu/wily/freefem++/wily

« back to all changes in this revision

Viewing changes to examples++-load/VTK_writer_3d.cpp

  • Committer: Package Import Robot
  • Author(s): Dimitrios Eftaxiopoulos, Dimitrios Eftaxiopoulos, Christophe Trophime
  • Date: 2013-09-12 00:02:58 UTC
  • mfrom: (1.2.1) (11.1.1 experimental)
  • Revision ID: package-import@ubuntu.com-20130912000258-aclq2zfa1svt0p3x
Tags: 3.25-1
http://bugs.debian.org/701161
[ Dimitrios Eftaxiopoulos ]
* Imported Upstream version 3.25 (Closes: #701161 #706714)
* Change installation directory of header-like *.idp files
  from /usr/lib/freefem++ to /usr/include/freefem++, in order
  to fix a lintian warning
* Update patch to examples++-load/Makefile.am in order to enable
  functioning of load *.so and include *.idp commands in *.edp
  scripts
* Delete patches to src/Graphics/sansgraph.cpp and
  src/Graphics/xglrgraph.cpp because they are not needed any more
* Fix lintian warning about missing LDFLAGS
* Override dh_auto_test in debian/rules, such that in case it is 
  used, it completes executing all *.edp example files, regardless
  of aborting on some of them
* Add libmetis-dev to build-deps in d/control
* Remove libparmetis-dev from build deps
* Add --parallel option to dh $@ in debian/rules
* Add hardening compilation flags to mpic++
* Allow testing of compiling and running the example files after build

[ Christophe Trophime ]
* update C. Trophime email
* add support for nlopt, ipopt - simplify debian/rules
* upload CT changes to 3.20
* add patch for configure
* add patch for examples++-mpi
* fix bamg install
* add corrected scripts to build plugins
* add patch for properly build examples++-load
* add lintian overrides for libfreefem++
* add some missing files
* update patches
* update rules
* reorder BuildDepends - comment out unsupported libs

Show diffs side-by-side

added added

removed removed

Lines of Context:
examples++-load/VTK_writer_3d.cpp
1
 
// -*- Mode : c++ -*-
2
 
//
3
 
// SUMMARY  : 
4
 
// USAGE    :        
5
 
// ORG      : 
6
 
// AUTHOR   : Cedric Ody (not an expert in c++)
7
 
// E-MAIL   : cedric.listes@gmail.com
8
 
// from the work of  Sala Lorenzo (Dxwriter)
9
 
 
10
 
#include "mode_open.hpp"
11
 
#include <iostream>
12
 
#include <cfloat>
13
 
#include <cmath>
14
 
#include <iterator>
15
 
using namespace std;
16
 
#include "ff++.hpp"
17
 
//using namespace Fem2D;
18
 
#include <set>
19
 
#include <vector>
20
 
//#include "msh3.hpp"
21
 
 
22
 
class VtkWriter 
23
 
{
24
 
 struct tsinfo
25
 
 {
26
 
  int imesh;//!<index of the mesh
27
 
  std::string name;
28
 
  std::vector<double> vecistant;
29
 
 };
30
 
 
31
 
private:
32
 
 std::vector<Mesh3*> _vecmesh;
33
 
 //std::vector<tsinfo> _vecofts;
34
 
 std::string _nameoffile;
35
 
 
36
 
 /*! This string contains the name of data file with \\ where there's a \ in the path*/
37
 
 std::string _nameofdatafile;
38
 
 
39
 
 //!files containing the data and the timeseries
40
 
 std::ofstream _ofdata;
41
 
 
42
 
public:
43
 
 VtkWriter() { std::cout << "Constructor of VtkWriter" << endl;  }
44
 
 void openfiles(const std::string& s)
45
 
  {
46
 
   _nameoffile=s;
47
 
   std::string tmp=s+".vtu";
48
 
   std::cout<<tmp<<" ";
49
 
   _ofdata.open(tmp.c_str(), std::ios_base::out);
50
 
   _nameofdatafile="";
51
 
   for(int i=0;i<tmp.length();++i)
52
 
    {
53
 
     if(tmp.at(i)=='\\')
54
 
      _nameofdatafile.append(1,'\\');
55
 
     _nameofdatafile.append(1,tmp.at(i));
56
 
    }
57
 
  }
58
 
 
59
 
 void addmesh(Mesh3* mesh)
60
 
  {
61
 
   Mesh3& Th(*mesh);
62
 
   _vecmesh.push_back(mesh);
63
 
   _ofdata.flags(std::ios_base::scientific);
64
 
   _ofdata.precision(15);
65
 
 
66
 
   _ofdata << "<?xml version=\"1.0\"?>" << std::endl;
67
 
   _ofdata << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" << std::endl;
68
 
   _ofdata << "<UnstructuredGrid>" << std::endl;
69
 
   _ofdata << "<Piece NumberOfPoints=\"" << Th.nv << "\" NumberOfCells=\"" << Th.nt << "\">" << std::endl;
70
 
   _ofdata << "<Points>" << std::endl;
71
 
   _ofdata << "<DataArray type=\"Float32\" Name=\"Position\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl;
72
 
   for(int k=0;k<Th.nv;++k) _ofdata << Th(k).x<<" "<<Th(k).y<< " "<<Th(k).z<<std::endl;
73
 
   _ofdata << "</DataArray>" << std::endl;
74
 
   _ofdata << "</Points>" << std::endl;
75
 
   _ofdata << "<Cells>" << std::endl;
76
 
   _ofdata << "<DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">" << std::endl;
77
 
   for(int i=0;i<Th.nt;++i)
78
 
    for (int j=0; j <4; j++) _ofdata << Th(i,j) << " " ;
79
 
   _ofdata << std::endl;
80
 
   _ofdata << "</DataArray>" << std::endl;
81
 
   _ofdata << "<DataArray type=\"Int32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"ascii\">" << std::endl;     
82
 
   for(int i=0;i<Th.nt;++i)  _ofdata << 4+4*(i) << " ";
83
 
   _ofdata << std::endl;
84
 
   _ofdata << "</DataArray>" << std::endl;
85
 
   _ofdata << "<DataArray type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"ascii\">" << std::endl;       
86
 
   for(int i=0;i<Th.nt;++i)  _ofdata << 10 << " ";
87
 
   _ofdata << std::endl;
88
 
   _ofdata << "</DataArray>" << std::endl;
89
 
   _ofdata << "</Cells>" << std::endl; 
90
 
   _ofdata << "<PointData >" << endl; 
91
 
  }
92
 
 
93
 
 /*!Add a field*/
94
 
 void addscalar(const string& nameoffield, Mesh3* mesh, const KN<double>&val)
95
 
  {
96
 
   _ofdata.flags(std::ios_base::scientific);
97
 
   _ofdata.precision(15);
98
 
 
99
 
   _ofdata << "<DataArray type=\"Float32\" Name=\""<<nameoffield<<"\" NumberOfComponents=\"1\" format=\"ascii\">" << std::endl;
100
 
   for(int i=0;i<val.size();++i) _ofdata<<checkprecision(val[i])<<std::endl;
101
 
   _ofdata << "</DataArray>" << std::endl;
102
 
   
103
 
   _ofdata.flush();
104
 
  }
105
 
 
106
 
 double checkprecision(double val)
107
 
  {
108
 
   double tmp;
109
 
   if ( val >= 0. ) tmp=max(0.,val);
110
 
   if ( val <  0. ) tmp=min(0.,val);
111
 
   return tmp;
112
 
  }
113
 
 
114
 
  /*!Add a field*/
115
 
 void addvector(const string& nameoffield, Mesh3* mesh, const KN<double>&val,
116
 
                const KN<double>&val2,const KN<double>&val3 )
117
 
  {
118
 
   _ofdata.flags(std::ios_base::scientific);
119
 
   _ofdata.precision(15);
120
 
 
121
 
   _ofdata << "<DataArray type=\"Float32\" Name=\""<<nameoffield<<"\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl;
122
 
   for(int i=0;i<val.size();++i) _ofdata<<checkprecision(val[i])<< " " << checkprecision(val2[i]) << " " << checkprecision(val3[i]) << std::endl;
123
 
   _ofdata << "</DataArray>" << std::endl;
124
 
   _ofdata.flush();
125
 
  }
126
 
 
127
 
 
128
 
 /*!Get the mesh associated with the series nameofts*/
129
 
 Mesh3* getmeshts(const string& nameofts)
130
 
  {
131
 
   return _vecmesh[0];
132
 
  }
133
 
 
134
 
 void init()
135
 
  {
136
 
   new(this)VtkWriter(); 
137
 
  }
138
 
 
139
 
 void destroy() 
140
 
  {
141
 
   if(_ofdata.is_open())
142
 
    {
143
 
     _ofdata << "</PointData>" << endl; 
144
 
     _ofdata << "<CellData>" << endl;   
145
 
     _ofdata << "</CellData>" << endl; 
146
 
     _ofdata << "</Piece>" << endl;
147
 
     _ofdata << "</UnstructuredGrid>" << endl;
148
 
     _ofdata << "</VTKFile>" << endl;
149
 
     _ofdata.close();           
150
 
    }
151
 
  } 
152
 
}; //End of class
153
 
 
154
 
class Vtkwritesol_Op : public E_F0mps 
155
 
{
156
 
public:
157
 
 typedef long  Result;
158
 
 Expression edx;
159
 
 Expression ename;//!<name of time series or field
160
 
 Expression et;//!<time
161
 
 long what; // 1 scalar, 2 vector, 3 symtensor
162
 
 long nbfloat; // 1 scalar, n vector (3D), n symtensor(3D)
163
 
 Expression evct,evct2,evct3;
164
 
 
165
 
public:
166
 
 Vtkwritesol_Op(const basicAC_F0 &  args) :  what(0), nbfloat(0)
167
 
  {
168
 
   evct=0;
169
 
   evct2=0;
170
 
   evct3=0;
171
 
   int nbofsol;
172
 
   int ddim=3;
173
 
   //There's no named parameter
174
 
   args.SetNameParam();
175
 
   if(args.size()!=3)
176
 
    {
177
 
     CompileError("Vtkwritesol accepts only 4 parameters");
178
 
    }
179
 
   if (BCastTo<VtkWriter *>(args[0])) edx = CastTo<VtkWriter *>(args[0]);
180
 
   if (BCastTo<string *>(args[1])) ename = CastTo<string *>(args[1]);
181
 
   
182
 
   if ( args[2].left()==atype<double>() )
183
 
    {
184
 
     what=1;
185
 
     nbfloat=1;
186
 
     evct=to<double>( args[2] );
187
 
    }
188
 
   else if ( args[2].left()==atype<double *>() )
189
 
    {
190
 
     what=1;
191
 
     nbfloat=1;
192
 
     evct=to<double>( args[2] );
193
 
    }
194
 
   else if ( BCastTo<pfer>(args[2]) )
195
 
    {
196
 
     what=1;
197
 
     nbfloat=1;
198
 
     evct=to<double>( args[2] );
199
 
    }
200
 
   else if ( args[2].left()==atype<E_Array>() )
201
 
    {
202
 
     std::cout << "Until now only scalar solution" << std::endl;
203
 
     
204
 
     int i=2;
205
 
     const E_Array * a0  = dynamic_cast<const E_Array *>( args[i].LeftValue() );
206
 
       
207
 
 
208
 
      if( a0->size() == 1){
209
 
      // scalar solution
210
 
      what=1;
211
 
      nbfloat=a0->size();
212
 
      evct = to<double>( (*a0)[0]);
213
 
      
214
 
     }
215
 
  
216
 
     if( a0->size() == ddim){
217
 
      // vector solution
218
 
      what=2;
219
 
      nbfloat=a0->size();
220
 
      evct = to<double>( (*a0)[0]);
221
 
      evct2 = to<double>( (*a0)[1]);
222
 
      evct3 = to<double>( (*a0)[2]);
223
 
      
224
 
     }
225
 
      cout << "Passed Until now only scalar solution" << std::endl;
226
 
    }
227
 
   else 
228
 
    {
229
 
     CompileError("savesol in 2D: Sorry no way to save this kind of data");
230
 
    }
231
 
   
232
 
  }
233
 
 static ArrayOfaType  typeargs() { return  ArrayOfaType(atype<VtkWriter *>(), atype<string *>(), true); }// all type
234
 
 static  E_F0 * f(const basicAC_F0 & args) { return new Vtkwritesol_Op(args);} 
235
 
 AnyType operator()(Stack stack)  const ;
236
 
}; // end of class
237
 
 
238
 
 
239
 
AnyType Vtkwritesol_Op::operator()(Stack stack)  const 
240
 
{ 
241
 
 MeshPoint *mp(MeshPointStack(stack)) , mps=*mp;
242
 
 VtkWriter &dx=*(GetAny<VtkWriter *>((*edx)(stack)));
243
 
 string &name=*(GetAny<string *>((*ename)(stack)));
244
 
 //double t=GetAny<double>((*et)(stack));
245
 
 Mesh3 &Th=*(dx.getmeshts(name));
246
 
 
247
 
 int nt = Th.nt;
248
 
 int nv = Th.nv;
249
 
 
250
 
 int nbsol=nv;
251
 
 long longdefault;
252
 
 
253
 
 KN<double> valsol(nbsol);
254
 
 valsol=0.;
255
 
 KN<int> takemesh(nbsol);
256
 
 takemesh=0;
257
 
 MeshPoint *mp3(MeshPointStack(stack));
258
 
 for (int it=0;it<nt;it++)
259
 
  {
260
 
  for(int iv=0;iv<4;iv++)
261
 
   {
262
 
    int i=Th(it,iv);
263
 
    mp3->setP(&Th,it,iv);                                       
264
 
    valsol[i] = valsol[i] + GetAny< double >((*evct)(stack));                   
265
 
    ++takemesh[i];
266
 
   }
267
 
  }
268
 
 for(int i=0; i<nbsol; i++)
269
 
  {
270
 
   valsol[i] /= takemesh[i]; 
271
 
  }
272
 
 
273
 
 //Writes valsol on the file file
274
 
 if (what==1) 
275
 
   dx.addscalar(name,&Th,valsol);
276
 
 
277
 
 if (what == 2)
278
 
  {
279
 
   KN<double> valsol2(nbsol);
280
 
   valsol2=0.;
281
 
   KN<int> takemesh(nbsol);
282
 
   takemesh=0;
283
 
   MeshPoint *mp3(MeshPointStack(stack));
284
 
   for (int it=0;it<nt;it++)
285
 
    {
286
 
     for(int iv=0;iv<4;iv++)
287
 
      {
288
 
       int i=Th(it,iv);
289
 
       mp3->setP(&Th,it,iv);                                    
290
 
       valsol2[i] = valsol2[i] + GetAny< double >((*evct2)(stack));                     
291
 
       ++takemesh[i];
292
 
      }
293
 
    }
294
 
   for(int i=0; i<nbsol; i++)
295
 
    {
296
 
     valsol2[i] /= takemesh[i]; 
297
 
    }
298
 
   
299
 
   {
300
 
    KN<double> valsol3(nbsol);
301
 
    valsol3=0.;
302
 
    KN<int> takemesh(nbsol);
303
 
    takemesh=0;
304
 
    MeshPoint *mp3(MeshPointStack(stack));
305
 
    for (int it=0;it<nt;it++)
306
 
     {
307
 
      for(int iv=0;iv<4;iv++)
308
 
       {
309
 
        int i=Th(it,iv);
310
 
        mp3->setP(&Th,it,iv);                                   
311
 
        valsol3[i] = valsol3[i] + GetAny< double >((*evct3)(stack));                    
312
 
        ++takemesh[i];
313
 
       }
314
 
     }
315
 
    for(int i=0; i<nbsol; i++)
316
 
     {
317
 
      valsol3[i] /= takemesh[i]; 
318
 
     }
319
 
 
320
 
   //Writes valsol on the file file
321
 
   dx.addvector(name,&Th,valsol,valsol2,valsol3);
322
 
 
323
 
      }
324
 
   
325
 
  }
326
 
 
327
 
 return longdefault;
328
 
 
329
 
}
330
 
 
331
 
 
332
 
 
333
 
// le vrai constructeur est la
334
 
VtkWriter* init_VtkWriter(VtkWriter * const &a, string * const & s)
335
 
{
336
 
 std::cout << "start init_VtkWriter" << std::endl;
337
 
 a->init();
338
 
 a->openfiles(*s);
339
 
 std::cout << "end init_VtkWriter" << std::endl;
340
 
 return a;
341
 
} 
342
 
 
343
 
void* call_addmesh( VtkWriter * const & mt, Mesh3* const & pTh) { mt->addmesh(pTh);}
344
 
 
345
 
//   Add the function name to the freefem++ table 
346
 
class Init 
347
 
{ 
348
 
public:
349
 
 Init();
350
 
};
351
 
 
352
 
LOADINIT(Init);
353
 
Init::Init()
354
 
{
355
 
 
356
 
 Dcl_Type<VtkWriter*>(InitP<VtkWriter>,Destroy<VtkWriter>); // declare deux nouveau type pour freefem++  un pointeur et 
357
 
 
358
 
 zzzfff->Add("VtkWriter",atype<VtkWriter*>()); // ajoute le type myType a freefem++ 
359
 
 // constructeur  d'un type myType  dans freefem 
360
 
 TheOperators->Add("<-", new OneOperator2_<VtkWriter*, VtkWriter* ,string*>(&init_VtkWriter));
361
 
 Global.Add("Vtkaddmesh","(",new OneOperator2_<void *, VtkWriter*, Mesh3*>(call_addmesh)); 
362
 
 Global.Add("Vtkaddscalar","(",new OneOperatorCode< Vtkwritesol_Op> );
363
 
 
364
 
}
 
1
// -*- Mode : c++ -*-
 
2
//
 
3
// SUMMARY  : 
 
4
// USAGE    :        
 
5
// ORG      : 
 
6
// AUTHOR   : Cedric Ody (not an expert in c++)
 
7
// E-MAIL   : cedric.listes@gmail.com
 
8
// from the work of  Sala Lorenzo (Dxwriter)
 
9
 
 
10
#include "mode_open.hpp"
 
11
#include <iostream>
 
12
#include <cfloat>
 
13
#include <cmath>
 
14
#include <iterator>
 
15
using namespace std;
 
16
#include "ff++.hpp"
 
17
//using namespace Fem2D;
 
18
#include <set>
 
19
#include <vector>
 
20
//#include "msh3.hpp"
 
21
 
 
22
class VtkWriter 
 
23
{
 
24
 struct tsinfo
 
25
 {
 
26
  int imesh;//!<index of the mesh
 
27
  std::string name;
 
28
  std::vector<double> vecistant;
 
29
 };
 
30
 
 
31
private:
 
32
 std::vector<Mesh3*> _vecmesh;
 
33
 //std::vector<tsinfo> _vecofts;
 
34
 std::string _nameoffile;
 
35
 
 
36
 /*! This string contains the name of data file with \\ where there's a \ in the path*/
 
37
 std::string _nameofdatafile;
 
38
 
 
39
 //!files containing the data and the timeseries
 
40
 std::ofstream _ofdata;
 
41
 
 
42
public:
 
43
 VtkWriter() { std::cout << "Constructor of VtkWriter" << endl;  }
 
44
 void openfiles(const std::string& s)
 
45
  {
 
46
   _nameoffile=s;
 
47
   std::string tmp=s+".vtu";
 
48
   std::cout<<tmp<<" ";
 
49
   _ofdata.open(tmp.c_str(), std::ios_base::out);
 
50
   _nameofdatafile="";
 
51
   for(int i=0;i<tmp.length();++i)
 
52
    {
 
53
     if(tmp.at(i)=='\\')
 
54
      _nameofdatafile.append(1,'\\');
 
55
     _nameofdatafile.append(1,tmp.at(i));
 
56
    }
 
57
  }
 
58
 
 
59
 void addmesh(Mesh3* mesh)
 
60
  {
 
61
   Mesh3& Th(*mesh);
 
62
   _vecmesh.push_back(mesh);
 
63
   _ofdata.flags(std::ios_base::scientific);
 
64
   _ofdata.precision(15);
 
65
 
 
66
   _ofdata << "<?xml version=\"1.0\"?>" << std::endl;
 
67
   _ofdata << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" << std::endl;
 
68
   _ofdata << "<UnstructuredGrid>" << std::endl;
 
69
   _ofdata << "<Piece NumberOfPoints=\"" << Th.nv << "\" NumberOfCells=\"" << Th.nt << "\">" << std::endl;
 
70
   _ofdata << "<Points>" << std::endl;
 
71
   _ofdata << "<DataArray type=\"Float32\" Name=\"Position\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl;
 
72
   for(int k=0;k<Th.nv;++k) _ofdata << Th(k).x<<" "<<Th(k).y<< " "<<Th(k).z<<std::endl;
 
73
   _ofdata << "</DataArray>" << std::endl;
 
74
   _ofdata << "</Points>" << std::endl;
 
75
   _ofdata << "<Cells>" << std::endl;
 
76
   _ofdata << "<DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">" << std::endl;
 
77
   for(int i=0;i<Th.nt;++i)
 
78
    for (int j=0; j <4; j++) _ofdata << Th(i,j) << " " ;
 
79
   _ofdata << std::endl;
 
80
   _ofdata << "</DataArray>" << std::endl;
 
81
   _ofdata << "<DataArray type=\"Int32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"ascii\">" << std::endl;     
 
82
   for(int i=0;i<Th.nt;++i)  _ofdata << 4+4*(i) << " ";
 
83
   _ofdata << std::endl;
 
84
   _ofdata << "</DataArray>" << std::endl;
 
85
   _ofdata << "<DataArray type=\"UInt8\" Name=\"types\" NumberOfComponents=\"1\" format=\"ascii\">" << std::endl;       
 
86
   for(int i=0;i<Th.nt;++i)  _ofdata << 10 << " ";
 
87
   _ofdata << std::endl;
 
88
   _ofdata << "</DataArray>" << std::endl;
 
89
   _ofdata << "</Cells>" << std::endl; 
 
90
   _ofdata << "<PointData >" << endl; 
 
91
  }
 
92
 
 
93
 /*!Add a field*/
 
94
 void addscalar(const string& nameoffield, Mesh3* mesh, const KN<double>&val)
 
95
  {
 
96
   _ofdata.flags(std::ios_base::scientific);
 
97
   _ofdata.precision(15);
 
98
 
 
99
   _ofdata << "<DataArray type=\"Float32\" Name=\""<<nameoffield<<"\" NumberOfComponents=\"1\" format=\"ascii\">" << std::endl;
 
100
   for(int i=0;i<val.size();++i) _ofdata<<checkprecision(val[i])<<std::endl;
 
101
   _ofdata << "</DataArray>" << std::endl;
 
102
   
 
103
   _ofdata.flush();
 
104
  }
 
105
 
 
106
 double checkprecision(double val)
 
107
  {
 
108
   double tmp;
 
109
   if ( val >= 0. ) tmp=max(0.,val);
 
110
   if ( val <  0. ) tmp=min(0.,val);
 
111
   return tmp;
 
112
  }
 
113
 
 
114
  /*!Add a field*/
 
115
 void addvector(const string& nameoffield, Mesh3* mesh, const KN<double>&val,
 
116
                const KN<double>&val2,const KN<double>&val3 )
 
117
  {
 
118
   _ofdata.flags(std::ios_base::scientific);
 
119
   _ofdata.precision(15);
 
120
 
 
121
   _ofdata << "<DataArray type=\"Float32\" Name=\""<<nameoffield<<"\" NumberOfComponents=\"3\" format=\"ascii\">" << std::endl;
 
122
   for(int i=0;i<val.size();++i) _ofdata<<checkprecision(val[i])<< " " << checkprecision(val2[i]) << " " << checkprecision(val3[i]) << std::endl;
 
123
   _ofdata << "</DataArray>" << std::endl;
 
124
   _ofdata.flush();
 
125
  }
 
126
 
 
127
 
 
128
 /*!Get the mesh associated with the series nameofts*/
 
129
 Mesh3* getmeshts(const string& nameofts)
 
130
  {
 
131
   return _vecmesh[0];
 
132
  }
 
133
 
 
134
 void init()
 
135
  {
 
136
   new(this)VtkWriter(); 
 
137
  }
 
138
 
 
139
 void destroy() 
 
140
  {
 
141
   if(_ofdata.is_open())
 
142
    {
 
143
     _ofdata << "</PointData>" << endl; 
 
144
     _ofdata << "<CellData>" << endl;   
 
145
     _ofdata << "</CellData>" << endl; 
 
146
     _ofdata << "</Piece>" << endl;
 
147
     _ofdata << "</UnstructuredGrid>" << endl;
 
148
     _ofdata << "</VTKFile>" << endl;
 
149
     _ofdata.close();           
 
150
    }
 
151
  } 
 
152
}; //End of class
 
153
 
 
154
class Vtkwritesol_Op : public E_F0mps 
 
155
{
 
156
public:
 
157
 typedef long  Result;
 
158
 Expression edx;
 
159
 Expression ename;//!<name of time series or field
 
160
 Expression et;//!<time
 
161
 long what; // 1 scalar, 2 vector, 3 symtensor
 
162
 long nbfloat; // 1 scalar, n vector (3D), n symtensor(3D)
 
163
 Expression evct,evct2,evct3;
 
164
 
 
165
public:
 
166
 Vtkwritesol_Op(const basicAC_F0 &  args) :  what(0), nbfloat(0)
 
167
  {
 
168
   evct=0;
 
169
   evct2=0;
 
170
   evct3=0;
 
171
   int nbofsol;
 
172
   int ddim=3;
 
173
   //There's no named parameter
 
174
   args.SetNameParam();
 
175
   if(args.size()!=3)
 
176
    {
 
177
     CompileError("Vtkwritesol accepts only 4 parameters");
 
178
    }
 
179
   if (BCastTo<VtkWriter *>(args[0])) edx = CastTo<VtkWriter *>(args[0]);
 
180
   if (BCastTo<string *>(args[1])) ename = CastTo<string *>(args[1]);
 
181
   
 
182
   if ( args[2].left()==atype<double>() )
 
183
    {
 
184
     what=1;
 
185
     nbfloat=1;
 
186
     evct=to<double>( args[2] );
 
187
    }
 
188
   else if ( args[2].left()==atype<double *>() )
 
189
    {
 
190
     what=1;
 
191
     nbfloat=1;
 
192
     evct=to<double>( args[2] );
 
193
    }
 
194
   else if ( BCastTo<pfer>(args[2]) )
 
195
    {
 
196
     what=1;
 
197
     nbfloat=1;
 
198
     evct=to<double>( args[2] );
 
199
    }
 
200
   else if ( args[2].left()==atype<E_Array>() )
 
201
    {
 
202
     std::cout << "Until now only scalar solution" << std::endl;
 
203
     
 
204
     int i=2;
 
205
     const E_Array * a0  = dynamic_cast<const E_Array *>( args[i].LeftValue() );
 
206
       
 
207
 
 
208
      if( a0->size() == 1){
 
209
      // scalar solution
 
210
      what=1;
 
211
      nbfloat=a0->size();
 
212
      evct = to<double>( (*a0)[0]);
 
213
      
 
214
     }
 
215
  
 
216
     if( a0->size() == ddim){
 
217
      // vector solution
 
218
      what=2;
 
219
      nbfloat=a0->size();
 
220
      evct = to<double>( (*a0)[0]);
 
221
      evct2 = to<double>( (*a0)[1]);
 
222
      evct3 = to<double>( (*a0)[2]);
 
223
      
 
224
     }
 
225
      cout << "Passed Until now only scalar solution" << std::endl;
 
226
    }
 
227
   else 
 
228
    {
 
229
     CompileError("savesol in 2D: Sorry no way to save this kind of data");
 
230
    }
 
231
   
 
232
  }
 
233
 static ArrayOfaType  typeargs() { return  ArrayOfaType(atype<VtkWriter *>(), atype<string *>(), true); }// all type
 
234
 static  E_F0 * f(const basicAC_F0 & args) { return new Vtkwritesol_Op(args);} 
 
235
 AnyType operator()(Stack stack)  const ;
 
236
}; // end of class
 
237
 
 
238
 
 
239
AnyType Vtkwritesol_Op::operator()(Stack stack)  const 
 
240
{ 
 
241
 MeshPoint *mp(MeshPointStack(stack)) , mps=*mp;
 
242
 VtkWriter &dx=*(GetAny<VtkWriter *>((*edx)(stack)));
 
243
 string &name=*(GetAny<string *>((*ename)(stack)));
 
244
 //double t=GetAny<double>((*et)(stack));
 
245
 Mesh3 &Th=*(dx.getmeshts(name));
 
246
 
 
247
 int nt = Th.nt;
 
248
 int nv = Th.nv;
 
249
 
 
250
 int nbsol=nv;
 
251
 long longdefault;
 
252
 
 
253
 KN<double> valsol(nbsol);
 
254
 valsol=0.;
 
255
 KN<int> takemesh(nbsol);
 
256
 takemesh=0;
 
257
 MeshPoint *mp3(MeshPointStack(stack));
 
258
 for (int it=0;it<nt;it++)
 
259
  {
 
260
  for(int iv=0;iv<4;iv++)
 
261
   {
 
262
    int i=Th(it,iv);
 
263
    mp3->setP(&Th,it,iv);                                       
 
264
    valsol[i] = valsol[i] + GetAny< double >((*evct)(stack));                   
 
265
    ++takemesh[i];
 
266
   }
 
267
  }
 
268
 for(int i=0; i<nbsol; i++)
 
269
  {
 
270
   valsol[i] /= takemesh[i]; 
 
271
  }
 
272
 
 
273
 //Writes valsol on the file file
 
274
 if (what==1) 
 
275
   dx.addscalar(name,&Th,valsol);
 
276
 
 
277
 if (what == 2)
 
278
  {
 
279
   KN<double> valsol2(nbsol);
 
280
   valsol2=0.;
 
281
   KN<int> takemesh(nbsol);
 
282
   takemesh=0;
 
283
   MeshPoint *mp3(MeshPointStack(stack));
 
284
   for (int it=0;it<nt;it++)
 
285
    {
 
286
     for(int iv=0;iv<4;iv++)
 
287
      {
 
288
       int i=Th(it,iv);
 
289
       mp3->setP(&Th,it,iv);                                    
 
290
       valsol2[i] = valsol2[i] + GetAny< double >((*evct2)(stack));                     
 
291
       ++takemesh[i];
 
292
      }
 
293
    }
 
294
   for(int i=0; i<nbsol; i++)
 
295
    {
 
296
     valsol2[i] /= takemesh[i]; 
 
297
    }
 
298
   
 
299
   {
 
300
    KN<double> valsol3(nbsol);
 
301
    valsol3=0.;
 
302
    KN<int> takemesh(nbsol);
 
303
    takemesh=0;
 
304
    MeshPoint *mp3(MeshPointStack(stack));
 
305
    for (int it=0;it<nt;it++)
 
306
     {
 
307
      for(int iv=0;iv<4;iv++)
 
308
       {
 
309
        int i=Th(it,iv);
 
310
        mp3->setP(&Th,it,iv);                                   
 
311
        valsol3[i] = valsol3[i] + GetAny< double >((*evct3)(stack));                    
 
312
        ++takemesh[i];
 
313
       }
 
314
     }
 
315
    for(int i=0; i<nbsol; i++)
 
316
     {
 
317
      valsol3[i] /= takemesh[i]; 
 
318
     }
 
319
 
 
320
   //Writes valsol on the file file
 
321
   dx.addvector(name,&Th,valsol,valsol2,valsol3);
 
322
 
 
323
      }
 
324
   
 
325
  }
 
326
 
 
327
 return longdefault;
 
328
 
 
329
}
 
330
 
 
331
 
 
332
 
 
333
// le vrai constructeur est la
 
334
VtkWriter* init_VtkWriter(VtkWriter * const &a, string * const & s)
 
335
{
 
336
 std::cout << "start init_VtkWriter" << std::endl;
 
337
 a->init();
 
338
 a->openfiles(*s);
 
339
 std::cout << "end init_VtkWriter" << std::endl;
 
340
 return a;
 
341
} 
 
342
 
 
343
void* call_addmesh( VtkWriter * const & mt, Mesh3* const & pTh) {
 
344
  mt->addmesh(pTh);
 
345
  return NULL;
 
346
}
 
347
 
 
348
//   Add the function name to the freefem++ table 
 
349
class Init 
 
350
{ 
 
351
public:
 
352
 Init();
 
353
};
 
354
 
 
355
LOADINIT(Init);
 
356
Init::Init()
 
357
{
 
358
 
 
359
 Dcl_Type<VtkWriter*>(InitP<VtkWriter>,Destroy<VtkWriter>); // declare deux nouveau type pour freefem++  un pointeur et 
 
360
 
 
361
 zzzfff->Add("VtkWriter",atype<VtkWriter*>()); // ajoute le type myType a freefem++ 
 
362
 // constructeur  d'un type myType  dans freefem 
 
363
 TheOperators->Add("<-", new OneOperator2_<VtkWriter*, VtkWriter* ,string*>(&init_VtkWriter));
 
364
 Global.Add("Vtkaddmesh","(",new OneOperator2_<void *, VtkWriter*, Mesh3*>(call_addmesh)); 
 
365
 Global.Add("Vtkaddscalar","(",new OneOperatorCode< Vtkwritesol_Op> );
 
366
 
 
367
}