GridWriter.hh 7.37 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
 * 
 ******************************************************************************/
20
21


22
23
24
25
26
27
28
#include <fstream>

#include "FixVec.h"
#include "DOFVector.h"
#include "BasisFunction.h"
#include "Mesh.h"

29
30
31
namespace AMDiS  { namespace io {
  
  namespace GridWriter
32
33
  {

34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
    template<typename T>
    void writeGrid(const WorldVector<double> *p,
				  int *numPoints, 
				  double *dist,
				  DOFVector<T> *vec,
				  const char *filename,
				  int outFilePrecision)
    {
      FUNCNAME("writeGrid()");

      TEST_EXIT(vec)("no dof vector\n");
      TEST_EXIT(filename)("no filename\n");
      TEST_EXIT(vec->getFeSpace()->getBasisFcts())("no basis fcts\n");
  
      std::ofstream outFile(filename);
      outFile.precision(outFilePrecision);
      outFile.setf(std::ios_base::scientific);    
      const Mesh *mesh = vec->getFeSpace()->getMesh();
      int dim = mesh->getDim();

      TEST_EXIT(dim == Global::getGeo(WORLD))("not for DIM != DIM_OF_WORLD\n");

      WorldVector<double> *basis = new WorldVector<double>[dim];
      double *lengthBasis = new double[dim];
      WorldVector<double> *step = new WorldVector<double>[3];
      for (int i = 0; i < 3; i++)
	step[i].set(0.0);

      for (int i = 0; i < dim; i++) {
	TEST_EXIT(numPoints[i] > 0)("numPoints < 1\n");
	TEST_EXIT(dist[i] >= 0)("dist < 0\n");
	lengthBasis[i] = 0;
66
67
      }

68
69
70
71
72
73
74
75
76
77
78
79
      WorldVector<double> curCoord;
      DimVec<double> bary(dim, NO_INIT);
      Element *elp;
      const BasisFunction *basFcts = vec->getFeSpace()->getBasisFcts();
      ElementVector uhLoc(basFcts->getNumber());

      // get basis of grid 
      for (int i = 0; i < dim; i++) {
	for (int j = 0; j < dim; j++) {
	  basis[i][j] = p[i + 1][j] ; 
	  lengthBasis[i] += basis[i][j] * basis[i][j];
	}
Praetorius, Simon's avatar
Praetorius, Simon committed
80
	lengthBasis[i] = std::sqrt(lengthBasis[i]);
81
82
      }

83
84
85
86
87
88
89
      // norm basis, get steps
      for (int i = 0; i < dim; i++) {
	for (int j = 0; j < dim; j++) {
	  basis[i][j] /= lengthBasis[i];
	  step[i][j] = basis[i][j] * dist[j];
	}
      }
90

91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
      /* write grid points */
      int localNumPoints[3] = {1, 1, 1};
      for (int i = 0; i < dim; i++)
	localNumPoints[i] = numPoints[i];    

      // Warning "Coords not in mesh domain" should be printed at most one time.
      bool warning = false;

      for (int i = 0; i < localNumPoints[0]; i++) {
	for (int j = 0; j < localNumPoints[1]; j++) {   // pseudo-loop for dim < 2
	  for (int k = 0; k < localNumPoints[2]; k++) { // pseudo-loop for dim < 3   
	    // clac current coords
	    for (int l = 0; l < dim; l++) {
	      curCoord[l] = p[0][l] 
		+ (i * step[0][l])
		+ (j * step[1][l])     // j = 0 for dim > 1
		+ (k * step[2][l]);    // k = 0 for dim > 2
	    }
109

110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
	    int inside = (const_cast<Mesh*>(mesh))->findElementAtPoint(curCoord, 
								      &elp, 
								      bary, 
								      NULL, NULL, NULL);

	    // write coords 
	    for (int l = 0; l < dim; l++)
	      outFile << curCoord[l] << " ";	  

	    if (!inside) {
	      if (!warning) {
		WARNING("Coords not in mesh domain\n");
		warning = true;
	      }
	      // write value 
	      outFile << "0.0" << std::endl;
	    } else {
	      // get value at coords 
	      vec->getLocalVector(elp, uhLoc);
	      double value = basFcts->evalUh(bary, uhLoc);
	    
	      // write value
	      outFile << value << std::endl;
133
134
	    }
	  }
135
136
	  if (localNumPoints[2] > 1) 
	    outFile << std::endl;
137
	}
138
	if (localNumPoints[1] > 1) 
139
	  outFile << std::endl;
140
      }
141
142
143
144
    
      delete [] basis;
      delete [] lengthBasis;
      delete [] step;
145
146
147
    }


148
149
150
151
152
153
154
155
156
    template<typename T>
    void writeGrid(const WorldVector<double> *p,
				  int *numPoints, 
				  double *dist,
				  std::vector<DOFVector<T>*> vec,
				  const char *filename,
				  int outFilePrecision)
    {
      FUNCNAME("writeGrid()");
157

158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
      TEST_EXIT(vec.size() > 0)("no dof vector\n");
      TEST_EXIT(filename)("no filename\n");
      TEST_EXIT(vec[0]->getFeSpace()->getBasisFcts())("no basis fcts\n");
  
      std::ofstream outFile(filename);
      outFile.precision(outFilePrecision);
      outFile.setf(std::ios_base::scientific);    
      const Mesh *mesh = vec[0]->getFeSpace()->getMesh();
      int dim = mesh->getDim();

      TEST_EXIT(dim == Global::getGeo(WORLD))("not for DIM != DIM_OF_WORLD\n");

      WorldVector<double> *basis = new WorldVector<double>[dim];
      double *lengthBasis = new double[dim];
      WorldVector<double> *step = new WorldVector<double>[3];
      for (int i = 0; i < 3; i++)
	step[i].set(0.0);

      for (int i = 0; i < dim; i++) {
	TEST_EXIT(numPoints[i] > 0)("numPoints < 1\n");
	TEST_EXIT(dist[i] >= 0)("dist < 0\n");
	lengthBasis[i] = 0;
      }
181

182
183
184
185
186
187
188
189
190
191
192
193
      WorldVector<double> curCoord;
      DimVec<double> bary(dim, NO_INIT);
      Element *elp;
      const BasisFunction *basFcts = vec[0]->getFeSpace()->getBasisFcts();
      ElementVector uhLoc(basFcts->getNumber());

      // get basis of grid 
      for (int i = 0; i < dim; i++) {
	for (int j = 0; j < dim; j++) {
	  basis[i][j] = p[i + 1][j] ;
	  lengthBasis[i] += basis[i][j] * basis[i][j];
	}
Praetorius, Simon's avatar
Praetorius, Simon committed
194
	lengthBasis[i] = std::sqrt(lengthBasis[i]);
195
      }
196

197
198
199
200
201
202
203
      // norm basis, get steps
      for (int i = 0; i < dim; i++) {
	for (int j = 0; j < dim; j++) {
	  basis[i][j] /= lengthBasis[i];
	  step[i][j] = basis[i][j] * dist[j];
	}
      }
204

205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
      /* write grid points */
      int localNumPoints[3] = {1, 1, 1};
      for (int i = 0; i < dim; i++)
	localNumPoints[i] = numPoints[i];    

      // Warning "Coords not in mesh domain" should be printed at most one time.
      bool warning = false;

      for (int i = 0; i < localNumPoints[0]; i++) {
	for (int j = 0; j < localNumPoints[1]; j++) {   // pseudo-loop for dim < 2
	  for (int k = 0; k < localNumPoints[2]; k++) { // pseudo-loop for dim < 3   
	    // clac current coords
	    for (int l = 0; l < dim; l++) {
	      curCoord[l] = p[0][l] 
		+ (i * step[0][l])
		+ (j * step[1][l])     // j = 0 for dim > 1
		+ (k * step[2][l]);    // k = 0 for dim > 2
222
	    }
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249

	    int inside = (const_cast<Mesh*>(mesh))->findElementAtPoint(curCoord, 
								      &elp, 
								      bary, 
								      NULL, NULL, NULL);

	    // write coords 
	    for (int l = 0; l < dim; l++)
	      outFile << curCoord[l] << " ";	  

	    if (!inside) {
	      if (!warning) {
		WARNING("Coords not in mesh domain\n");
		warning = true;
	      }
	      // write value 
	      outFile << "0.0" << std::endl;
	    } else {
	      // get value at coords
	      double value; 
	      for (int rr = 0; rr < vec.size(); rr++) {
		vec[rr]->getLocalVector(elp, uhLoc);
		value = basFcts->evalUh(bary, uhLoc);
		// write value
		outFile << value <<"\t";
	      }
	      outFile << std::endl;
250
	    }
251
	  }
252
253
	  if (localNumPoints[2] > 1) 
	    outFile << std::endl;
254
	}
255
	if (localNumPoints[1] > 1) 
256
257
	  outFile << std::endl;
      }
258
259
260
261
    
      delete [] basis;
      delete [] lengthBasis;
      delete [] step;
262
    }
263
264
265
    
  } // end namespace GridWriter
} } // end namespace io, AMDiS