[cig-commits] r13743 - cs/cigma/trunk/src

[luis at geodynamics.org]

Author: luis
Date: 2008-12-17 02:32:51 -0800 (Wed, 17 Dec 2008)
New Revision: 13743

Modified:
   cs/cigma/trunk/src/core_writers.cpp
Log:
Updated WriteResiduals to include the mesh when writing to VTK file

Modified: cs/cigma/trunk/src/core_writers.cpp
===================================================================
--- cs/cigma/trunk/src/core_writers.cpp	2008-12-17 10:32:49 UTC (rev 13742)
+++ cs/cigma/trunk/src/core_writers.cpp	2008-12-17 10:32:51 UTC (rev 13743)
@@ -2,7 +2,7 @@
 #include "io_text_writer.h"
 #include "io_hdf5_writer.h"
 #include "io_vtk_writer.h"
-#include "tri_logger.hpp"
+#include "Common.h"
 using namespace cigma;
 
 #include <boost/shared_ptr.hpp>
@@ -12,6 +12,9 @@
 #include <sstream>
 using namespace std;
 
+#ifdef HAVE_VTK
+#include "vtkCellType.h"
+#endif
 
 void cigma::WriteArray(const DataPath& path, const cigma::array<double>& x)
 {
@@ -186,43 +189,134 @@
     throw cigma::Exception("WriteDofs", "Need implementation");
 }
 
+
+// XXX: move this elsewhere
+static int vtkcelltype(Cell::type cellType)
+{
+#ifdef HAVE_VTK
+    switch (cellType)
+    {
+    case Cell::HEX8:  return VTK_HEXAHEDRON;
+    case Cell::TET4:  return VTK_TETRA;
+    case Cell::QUAD4: return VTK_QUAD;
+    case Cell::TRI3:  return VTK_TRIANGLE;
+    }
+    return VTK_EMPTY_CELL;
+#else
+    return 0;
+#endif
+}
+
 void cigma::WriteResiduals(const DataPath& path, shared_ptr<Residuals> residuals)
 {
     TRI_LOG_STR("cigma::WriteResiduals()");
     TRI_LOG(path);
 
+    if (!residuals)
+    {
+        throw cigma::Exception("WriteResiduals", "No residuals to write");
+    }
+
+    assert(residuals->epsilon != 0);
+    assert(residuals->nel > 0);
+
     std::string filename = path.filename();
+    std::string location = path.location();
+
     shared_ptr<FileWriter> writer = FileWriter::New(filename, "w");
     const FileWriter::WriterType wt = writer->getWriterType();
     if (wt == FileWriter::HDF5_FILE_WRITER)
     {
-        int status;
         HDF5_Writer *hdf5_writer = static_cast<HDF5_Writer*>(&(*writer));
 
-        string eps_loc = path.location();
-        if (eps_loc == "") { eps_loc += "/residuals"; }
-        // XXX: check whether eps_loc is a group or a dataset
-
+        string eps_loc = location;
+        if (hdf5_writer->locationIsGroup(eps_loc.c_str()))
+        {
+            eps_loc += "/residuals";
+        }
         TRI_LOG(eps_loc);
 
-        const double* epsilon = residuals->epsilon;
-        status = hdf5_writer->writeDataset(eps_loc.c_str(), epsilon, residuals->nel, 1);
+        int status = hdf5_writer->writeDataset(eps_loc.c_str(), residuals->epsilon, residuals->nel, 1);
+        TRI_LOG(status);
         if (status < 0)
         {
             ostringstream stream;
-            stream << "Could not write residuals to HDF5 file '" << filename << "'" << std::ends;
+            stream << "Could not write residuals to location '"
+                   << eps_loc << "' in the HDF5 file '"
+                   << filename << "'" << std::ends;
             throw cigma::Exception("WriteResiduals", stream.str());
         }
     }
     else if (wt == FileWriter::VTK_FILE_WRITER)
     {
+        // 
         // XXX: Build up a cell-based scalar array
-        throw cigma::Exception("WriteResiduals", "Need VTK writer");
+        // Move this code elsewhere, later.
+        //
+
+        VtkWriter *vtk_writer = static_cast<VtkWriter*>(&(*writer));
+
+        int status;
+
+        if (residuals->meshPart)
+        {
+            int i,j;
+
+            // coordinates first
+            TRI_LOG_STR("...writing vtk points");
+            const int npts = residuals->meshPart->coords->n_points();
+            const int ndim = residuals->meshPart->coords->n_dim();
+            TRI_LOG(npts);
+            TRI_LOG(ndim);
+            NodeCoordinates *nc = static_cast<NodeCoordinates*>(&(*(residuals->meshPart->coords)));
+            cigma::array<double> points(npts, 3);
+            points.init_value(0);
+            for (i = 0; i < npts; i++)
+            {
+                for (j = 0; j < 3; j++)
+                {
+                    points(i,j) = nc->getPoint(i,j);
+                }
+            }
+            vtk_writer->_writePoints(points._data, npts, 3);
+            points.reinit(0,0);
+
+            // cells next..
+            TRI_LOG_STR("...writing vtk cells");
+            const int nel = residuals->meshPart->connect->n_cells();
+            const int nno = residuals->meshPart->connect->n_dofs();
+            TRI_LOG(nel);
+            TRI_LOG(nno);
+            ElementBlock *eb = static_cast<ElementBlock*>(&(*(residuals->meshPart->connect)));
+            cigma::array<int> cells(nel, nno);
+            for (i = 0; i < nel; i++)
+            {
+                for (j = 0; j < nno; j++)
+                {
+                    cells(i,j) = eb->getId(i,j);
+                }
+            }
+            vtk_writer->_writeCells(cells._data, nel, nno);
+            cells.reinit(0,0);
+
+            // finally, cell types
+            TRI_LOG("...writing vtk cell types");
+            int vtkCellType = vtkcelltype(residuals->meshPart->getCellType());
+            TRI_LOG(vtkCellType);
+            vtk_writer->_writeCellTypes(nel, vtkCellType);
+        }
+        else
+        {
+            TRI_LOG_STR("Warning! VTK file will not contain mesh information");
+        }
+
+        TRI_LOG_STR("...writing residuals dataset as vtk cell data");
+        vtk_writer->_writeCellData(location.c_str(), residuals->epsilon, residuals->nel, 1);
     }
     else if (wt == FileWriter::TEXT_FILE_WRITER)
     {
-        const double *epsilon = residuals->epsilon;
-        int status = writer->writeDataset("", epsilon, residuals->nel, 1);
+        int status = writer->writeDataset("", residuals->epsilon, residuals->nel, 1);
+        TRI_LOG(status);
         if (status < 0)
         {
             ostringstream stream;