GridWriter.hh 3.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#include <fstream>

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

namespace AMDiS {

  template<typename T>
  void GridWriter<T>::writeGrid(const WorldVector<double> *p,
				int               *numPoints, 
				double            *dist,
				DOFVector<T>      *vec,
				char*              filename)
  {
    FUNCNAME("GridWriter<T>::writeGrid()");
 
19
    std::ofstream *outFile = new std::ofstream(filename);
20
21
22
23
24
25
26
27
28
    TEST_EXIT(outFile)("can't open file %s\n", filename);

    const Mesh *mesh = vec->getFESpace()->getMesh();
    int dim = mesh->getDim();

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


Thomas Witkowski's avatar
Thomas Witkowski committed
29
30
31
    WorldVector<double> *basis = new WorldVector<double>[dim];
    double *lengthBasis = new double[dim];
    WorldVector<double> *step = new WorldVector<double>[3];
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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107

    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;
    }

    WorldVector<double>          curCoord;
    DimVec<double>       bary(dim, NO_INIT);
    Element             *elp;
    const BasisFunction *basFcts = NULL;
    int                  inside;
    double               value;
    const double        *uhLoc;

    TEST_EXIT(vec)("no dof vector\n");
    TEST_EXIT(filename)("no filename\n");
    TEST_EXIT(vec->getFESpace()->getBasisFcts())("no basis fcts\n");

    mesh = vec->getFESpace()->getMesh();
    basFcts = vec->getFESpace()->getBasisFcts();

    /* 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] - p[0][j];
	lengthBasis[i] += basis[i][j] * basis[i][j];
      }
      lengthBasis[i] = sqrt(lengthBasis[i]);
    }

    /* 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];
      }
    }

    /* 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
	  }

	  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 */
108
	    (*outFile) << "0.0" << std::endl;
109
110
111
112
113
114
	  } else {
	    /* get value at coords */
	    uhLoc = vec->getLocalVector(elp, NULL);
	    value = basFcts->evalUh(bary, uhLoc);
	  
	    /* write value */
115
	    (*outFile) << value << std::endl;
116
117
	  }
	}
118
	if (localNumPoints[2] > 1) (*outFile) << std::endl;
119
      }
120
      if (localNumPoints[1] > 1) (*outFile) << std::endl;
121
122
    }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
123
    delete [] basis;
Thomas Witkowski's avatar
Thomas Witkowski committed
124
    delete [] lengthBasis;
Thomas Witkowski's avatar
Thomas Witkowski committed
125
    delete [] step;
126
127
128
129
    delete outFile;
  }

}