← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/kernel/io/MFile.cpp

  • Committer: Johannes Ring
  • Date: 2008-03-05 22:43:06 UTC
  • Revision ID: johannr@simula.no-20080305224306-2npsdyhfdpl2esji
The BIG commit!

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/kernel/io/MFile.cpp
1
 
// Copyright (C) 2003-2007 Johan Hoffman and Anders Logg.
2
 
// Licensed under the GNU LGPL Version 2.1.
3
 
//
4
 
// Modified by Garth N. Wells 2005
5
 
//
6
 
// First added:  2003-05-06
7
 
// Last changed: 2007-05-02
8
 
 
9
 
#include <dolfin/dolfin_log.h>
10
 
#include <dolfin/Vector.h>
11
 
#include <dolfin/Matrix.h>
12
 
#include <dolfin/Mesh.h>
13
 
#include <dolfin/Function.h>
14
 
#include <dolfin/Sample.h>
15
 
#include <dolfin/Vertex.h>
16
 
#include <dolfin/Cell.h>
17
 
#include <dolfin/MFile.h>
18
 
 
19
 
using namespace dolfin;
20
 
 
21
 
//-----------------------------------------------------------------------------
22
 
MFile::MFile(const std::string filename) : 
23
 
  GenericFile(filename)
24
 
{
25
 
  type = "Octave/MATLAB";
26
 
}
27
 
//-----------------------------------------------------------------------------
28
 
MFile::~MFile()
29
 
{
30
 
  // Do nothing
31
 
}
32
 
//-----------------------------------------------------------------------------
33
 
void MFile::operator<<(Vector& x)
34
 
{
35
 
  error("Function output in Matlab/Octave format is broken. Need to access vector elements correctly.");
36
 
/*
37
 
  // Open file
38
 
  FILE *fp = fopen(filename.c_str(), "a");
39
 
  
40
 
  // Write vector
41
 
  fprintf(fp, "%s = [", x.name().c_str());
42
 
  for (unsigned int i = 0; i < x.size(); i++)
43
 
  {
44
 
    // FIXME: This is a slow way to access PETSc vectors. Need a fast way 
45
 
    //        which is consistent for different vector types.
46
 
    real temp = x(i);
47
 
    fprintf(fp, " %.15g", temp);
48
 
  }
49
 
  fprintf(fp, " ];\n");
50
 
  
51
 
  // Close file
52
 
  fclose(fp);
53
 
 
54
 
  message("Saved vector %s (%s) to file %s in Octave/MATLAB format.",
55
 
              x.name().c_str(), x.label().c_str(), filename.c_str());
56
 
*/
57
 
}
58
 
//-----------------------------------------------------------------------------
59
 
void MFile::operator<<(Mesh& mesh)
60
 
{
61
 
  Point p;
62
 
  
63
 
  // Open file
64
 
  FILE *fp = fopen(filename.c_str(), "a");
65
 
 
66
 
  // Create a list if we save the mesh a second time
67
 
  if ( counter == 1 )
68
 
  {
69
 
    fprintf(fp, "tmp = points;\n");
70
 
    fprintf(fp, "clear points\n");
71
 
    fprintf(fp, "points{1} = tmp;\n");
72
 
    fprintf(fp, "clear tmp\n\n");
73
 
 
74
 
    fprintf(fp, "tmp = cells;\n");
75
 
    fprintf(fp, "clear cells\n");
76
 
    fprintf(fp, "cells{1} = tmp;\n");
77
 
    fprintf(fp, "clear tmp\n\n");
78
 
 
79
 
    fprintf(fp, "tmp = edges;\n");
80
 
    fprintf(fp, "clear edges\n");
81
 
    fprintf(fp, "edges{1} = tmp;\n");
82
 
    fprintf(fp, "clear tmp\n\n");
83
 
  }
84
 
  
85
 
  // Write vertices
86
 
  if ( counter == 0 )
87
 
    fprintf(fp,"points = [");
88
 
  else
89
 
    fprintf(fp,"points{%u} = [", counter + 1);
90
 
  for (VertexIterator v(mesh); !v.end();)
91
 
  {
92
 
    p = v->point();
93
 
    
94
 
    ++v;
95
 
    if ( mesh.type().cellType() == CellType::triangle )
96
 
    {
97
 
      if ( v.end() )
98
 
        fprintf(fp,"%.15f %.15f]';\n", p.x(), p.y());
99
 
      else
100
 
        fprintf(fp,"%.15f %.15f\n", p.x(), p.y() );
101
 
    }
102
 
    else {
103
 
      if ( v.end() )
104
 
        fprintf(fp,"%.15f %.15f %.15f]';\n", p.x(), p.y(), p.z());
105
 
      else
106
 
        fprintf(fp,"%.15f %.15f %.15f\n", p.x(), p.y(), p.z());
107
 
    }
108
 
  }
109
 
  fprintf(fp,"\n");
110
 
  
111
 
  // Write cells
112
 
  if ( counter == 0 )
113
 
    fprintf(fp,"cells = [");
114
 
  else
115
 
    fprintf(fp,"cells{%u} = [", counter + 1);
116
 
  for (CellIterator c(mesh); !c.end();)
117
 
  {
118
 
    for (VertexIterator v(*c); !v.end(); ++v)
119
 
      fprintf(fp, "%u ", (v->index()) + 1 );    
120
 
    ++c;
121
 
    if ( c.end() )
122
 
      fprintf(fp, "]';\n");
123
 
    else
124
 
      fprintf(fp, "\n");
125
 
  }
126
 
  fprintf(fp,"\n");
127
 
 
128
 
  // FIXME: Save all edges correctly
129
 
  // Write edges (to make the pdeplot routines happy)
130
 
  if ( counter == 0 )
131
 
    fprintf(fp,"edges = [1;2;0;0;0;0;0];\n\n");
132
 
  else
133
 
    fprintf(fp,"edges{%u} = [1;2;0;0;0;0;0];\n\n", counter + 1);
134
 
  
135
 
  // Close file
136
 
  fclose(fp);
137
 
 
138
 
//  // Increase the number of times we have saved the mesh
139
 
//  // FIXME: Count number of meshes saved to this file, rather
140
 
//  // than the number of times this specific mesh has been saved.
141
 
//  ++mesh;
142
 
  
143
 
  // Increase the number of meshes saved to this file
144
 
  counter++;
145
 
 
146
 
  message(1, "Saved mesh %s (%s) to file %s in Octave/MATLAB format.",
147
 
          mesh.name().c_str(), mesh.label().c_str(), filename.c_str());
148
 
}
149
 
//-----------------------------------------------------------------------------
150
 
void MFile::operator<<(Function& u)
151
 
{
152
 
  error("Function output in Matlab/Octave format not implemented for new Function.");
153
 
/*
154
 
  // Write mesh the first time
155
 
  if ( counter1 == 0 )
156
 
    *this << u.mesh();
157
 
  
158
 
  // Open file
159
 
  FILE *fp = fopen(filename.c_str(), "a");
160
 
  
161
 
  // Move old vector into list if we are saving a new value
162
 
  if ( counter1 == 1 )
163
 
  {
164
 
    fprintf(fp, "tmp = %s;\n", u.name().c_str());
165
 
    fprintf(fp, "clear %s\n", u.name().c_str());
166
 
    fprintf(fp, "%s{1} = tmp;\n", u.name().c_str());
167
 
    fprintf(fp, "clear tmp\n\n");
168
 
  }
169
 
 
170
 
  // Write vector
171
 
  if ( counter1 == 0 )
172
 
  {
173
 
    fprintf(fp, "%s = [", u.name().c_str());
174
 
    for (unsigned int i = 0; i < u.vectordim(); i++)
175
 
    { 
176
 
      for (VertexIterator v(u.mesh()); !v.end(); ++v)
177
 
        fprintf(fp, " %.15f", u(*v, i));
178
 
        fprintf(fp, ";");
179
 
    }
180
 
    fprintf(fp, " ]';\n\n");
181
 
  }
182
 
  else
183
 
  {
184
 
    fprintf(fp, "%s{%u} = [", u.name().c_str(), counter1 + 1);
185
 
    for (unsigned int i = 0; i < u.vectordim(); i++)
186
 
    { 
187
 
      for (VertexIterator v(u.mesh()); !v.end(); ++v)
188
 
        fprintf(fp, " %.15f", u(*v, i));
189
 
        fprintf(fp, ";");
190
 
    }
191
 
    fprintf(fp, " ]';\n\n");
192
 
  }
193
 
  
194
 
  // Close file
195
 
  fclose(fp);
196
 
  
197
 
  // Increase the number of times we have saved the function
198
 
  counter1++;
199
 
 
200
 
  cout << "Saved function " << u.name() << " (" << u.label()
201
 
       << ") to file " << filename << " in Octave/Matlab format." << endl;
202
 
*/
203
 
}
204
 
//-----------------------------------------------------------------------------
205
 
void MFile::operator<<(Sample& sample)
206
 
{
207
 
  // Open file
208
 
  FILE *fp = fopen(filename.c_str(), "a");
209
 
 
210
 
  // Initialize data structures first time
211
 
  if ( counter2 == 0 )
212
 
  {
213
 
    fprintf(fp, "t = [];\n");
214
 
    fprintf(fp, "%s = [];\n", sample.name().c_str());
215
 
    fprintf(fp, "k = [];\n");
216
 
    fprintf(fp, "r = [];\n");
217
 
    fprintf(fp, "\n");
218
 
  }
219
 
  
220
 
  // Save time
221
 
  fprintf(fp, "t = [t %.15e];\n", sample.t());
222
 
 
223
 
  // Save solution
224
 
  fprintf(fp, "tmp = [ ");
225
 
  for (unsigned int i = 0; i < sample.size(); i++)
226
 
    fprintf(fp, "%.15e ", sample.u(i));  
227
 
  fprintf(fp, "];\n");
228
 
  fprintf(fp, "%s = [%s tmp'];\n", 
229
 
          sample.name().c_str(), sample.name().c_str());
230
 
  //fprintf(fp, "clear tmp;\n");
231
 
 
232
 
  // Save time steps
233
 
  fprintf(fp, "tmp = [ ");
234
 
  for (unsigned int i = 0; i < sample.size(); i++)
235
 
    fprintf(fp, "%.15e ", sample.k(i));
236
 
 
237
 
  fprintf(fp, "];\n");
238
 
  fprintf(fp, "k = [k tmp'];\n");
239
 
  //fprintf(fp, "clear tmp;\n");
240
 
 
241
 
  // Save residuals
242
 
  fprintf(fp, "tmp = [ ");
243
 
  for (unsigned int i = 0; i < sample.size(); i++)
244
 
    fprintf(fp, "%.15e ", sample.r(i));
245
 
 
246
 
  fprintf(fp, "];\n");
247
 
  fprintf(fp, "r = [r tmp'];\n");
248
 
  fprintf(fp, "clear tmp;\n");
249
 
 
250
 
  // Increase frame counter
251
 
  counter2++;
252
 
  
253
 
  // Close file
254
 
  fclose(fp);
255
 
}
256
 
//-----------------------------------------------------------------------------