GridWriter extended to vector of Dof's and add control for precission for output

AMDiS/src/GridWriter.h

@@ -35,7 +35,7 @@ namespace AMDiS {
   class GridWriter {
   public:
     /** \brief
-     * Produces a grid-based output of a dof-vector
+     * Produces a grid-based output of a dof-vector or vector of dof-vectors 
      * \param p array of world coordinates. 
      * - p[0] defines origin of grid-system.
      * - p[1] defines first basis vector by p[1] - p[0].
@@ -51,12 +51,21 @@ namespace AMDiS {
      * - dist[2] distance in the third direction
      * \param vec DOFVector which is interpolated to the grid.
      * \param filename file to which the grid will be written
+     * \param outFilePrecision precision used to write file
      */
     static void writeGrid(const WorldVector<double> *p,
 			  int               *numPoints, 
 			  double            *dist,
 			  DOFVector<T>      *vec,
-			  char*              filename); 
+			  const char*              filename,
+	                  int outFilePrecision=6); 
+
+    static void writeGrid(const WorldVector<double> *p,
+			  int               *numPoints, 
+			  double            *dist,
+			  std::vector< DOFVector<T> * >     vec,
+			  const char*              filename,
+	                  int outFilePrecision=6); 
   };
 
 }

AMDiS/src/GridWriter.hh

@@ -12,7 +12,8 @@ namespace AMDiS {
 				int *numPoints, 
 				double *dist,
 				DOFVector<T> *vec,
-				char *filename)
+				const char *filename,
+                                int outFilePrecision)
   {
     FUNCNAME("GridWriter<T>::writeGrid()");
 
@@ -20,7 +21,9 @@ namespace AMDiS {
     TEST_EXIT(filename)("no filename\n");
     TEST_EXIT(vec->getFESpace()->getBasisFcts())("no basis fcts\n");
  
-    std::ofstream outFile = std::ofstream(filename);
+    std::ofstream outFile(filename);
+    outFile.precision(outFilePrecision);
+    outFile.setf(std::ios_base::scientific);    
     const Mesh *mesh = vec->getFESpace()->getMesh();
     int dim = mesh->getDim();
 
@@ -115,4 +118,119 @@ namespace AMDiS {
     delete [] step;
   }
 
+
+  template<typename T>
+  void GridWriter<T>::writeGrid(const WorldVector<double> *p,
+				int *numPoints, 
+				double *dist,
+				std::vector< DOFVector<T> * >     vec,
+				const char *filename,
+                                int outFilePrecision)
+  {
+    FUNCNAME("GridWriter<T>::writeGrid()");
+
+    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 < 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 = vec[0]->getFESpace()->getBasisFcts();
+    const double *uhLoc;
+
+    // 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
+	  }
+
+	  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++) {
+		  uhLoc = vec[rr]->getLocalVector(elp, NULL);
+		  value = basFcts->evalUh(bary, uhLoc);
+		  /* write value */
+		  outFile << value <<"\t";
+	      }
+	      outFile << std::endl;
+	  }
+	}
+	if (localNumPoints[2] > 1) 
+	  outFile << std::endl;
+      }
+      if (localNumPoints[1] > 1) 
+	outFile << std::endl;
+    }
+  
+    delete [] basis;
+    delete [] lengthBasis;
+    delete [] step;
+  }
+
 }