GridWriter.hh 3.09 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
#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,
12
13
14
15
				int *numPoints, 
				double *dist,
				DOFVector<T> *vec,
				char *filename)
16
17
18
  {
    FUNCNAME("GridWriter<T>::writeGrid()");

19
20
21
22
23
    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 = std::ofstream(filename);
24
25
26
    const Mesh *mesh = vec->getFESpace()->getMesh();
    int dim = mesh->getDim();

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

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

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

39
40
41
42
43
    WorldVector<double> curCoord;
    DimVec<double> bary(dim, NO_INIT);
    Element *elp;
    const BasisFunction *basFcts = vec->getFESpace()->getBasisFcts();
    const double *uhLoc;
44

45
    // get basis of grid 
46
47
    for (int i = 0; i < dim; i++) {
      for (int j = 0; j < dim; j++) {
48
	basis[i][j] = p[i + 1][j] - p[0][j];
49
50
51
52
53
	lengthBasis[i] += basis[i][j] * basis[i][j];
      }
      lengthBasis[i] = sqrt(lengthBasis[i]);
    }

54
    // norm basis, get steps
55
56
57
58
59
60
61
62
    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 */
63
64
65
    int localNumPoints[3] = {1, 1, 1};
    for (int i = 0; i < dim; i++)
      localNumPoints[i] = numPoints[i];    
66
67
68
69
70
71
72

    // 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   
73
	  // clac current coords
74
75
76
77
78
79
80
	  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
	  }

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

86
87
88
	  // write coords 
	  for (int l = 0; l < dim; l++)
	    outFile << curCoord[l] << " ";	  
89
90
91
92
93
94

	  if (!inside) {
	    if (!warning) {
	      WARNING("Coords not in mesh domain\n");
	      warning = true;
	    }
95
96
	    // write value 
	    outFile << "0.0" << std::endl;
97
	  } else {
98
	    // get value at coords 
99
	    uhLoc = vec->getLocalVector(elp, NULL);
100
	    double value = basFcts->evalUh(bary, uhLoc);
101
	  
102
103
	    // write value
	    outFile << value << std::endl;
104
105
	  }
	}
106
107
	if (localNumPoints[2] > 1) 
	  outFile << std::endl;
108
      }
109
110
      if (localNumPoints[1] > 1) 
	outFile << std::endl;
111
112
    }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
113
    delete [] basis;
Thomas Witkowski's avatar
Thomas Witkowski committed
114
    delete [] lengthBasis;
Thomas Witkowski's avatar
Thomas Witkowski committed
115
    delete [] step;
116
117
118
  }

}