[cig-commits] r14083 - cs/cigma/trunk/src
luis at geodynamics.org
luis at geodynamics.org
Wed Feb 18 08:14:36 PST 2009
Author: luis
Date: 2009-02-18 08:14:36 -0800 (Wed, 18 Feb 2009)
New Revision: 14083
Modified:
cs/cigma/trunk/src/core_writers.cpp
cs/cigma/trunk/src/core_writers.h
Log:
Add WriteCellCentroids() and updated WriteResiduals()
Modified: cs/cigma/trunk/src/core_writers.cpp
===================================================================
--- cs/cigma/trunk/src/core_writers.cpp 2009-02-18 16:14:34 UTC (rev 14082)
+++ cs/cigma/trunk/src/core_writers.cpp 2009-02-18 16:14:36 UTC (rev 14083)
@@ -6,20 +6,20 @@
#include "nc_array.h"
#include "eb_array.h"
#include "Common.h"
-using namespace cigma;
#include <boost/shared_ptr.hpp>
-using boost::shared_ptr;
-
#include <string>
#include <sstream>
#include <valarray>
-using namespace std;
#ifdef HAVE_VTK
#include "vtkCellType.h"
#endif
+using namespace cigma;
+using namespace std;
+using boost::shared_ptr;
+
void cigma::WriteArray(const DataPath& path, const cigma::array<double>& x, bool overwrite)
{
TRI_LOG_STR("cigma::WriteArray()");
@@ -396,8 +396,74 @@
}
}
-void cigma::WriteResiduals(const DataPath& path, shared_ptr<Residuals> residuals, bool overwrite)
+void cigma::WriteCellCentroids(const DataPath& path, boost::shared_ptr<MeshPart> M, bool overwrite)
{
+ TRI_LOG_STR("cigma::WriteCellCentroids()");
+ TRI_LOG(path);
+
+ // calculate centroids (XXX: move this calculation to MeshPart class)
+ int ncells = M->n_cells();
+ int ndim = M->n_dim();
+ std::valarray<double> centroids;
+ shared_ptr<Cell> cell = M->getCell();
+ if (cell)
+ {
+ centroids.resize(ncells * ndim);
+ for (int e = 0; e < ncells; e++)
+ {
+ M->getCell(e, *cell);
+
+ double c[3] = {0,0,0};
+ cell->centroid(c);
+ for (int i = 0; i < ndim; i++)
+ {
+ centroids[ndim*e + i] = c[i];
+ }
+ }
+ }
+
+ // write it out
+ string filename = path.filename();
+ string location = path.location();
+
+ shared_ptr<FileWriter> writer = FileWriter::New(filename, "a");
+ writer->setOverwriteMode(overwrite);
+
+ const FileWriter::WriterType wt = writer->getWriterType();
+ if (wt == FileWriter::HDF5_FILE_WRITER)
+ {
+ int status = writer->writeDataset(location.c_str(), ¢roids[0], ncells, ndim);
+ if (status < 0)
+ {
+ std::ostringstream stream;
+ stream << "Could not write cell centroids to HDF5 file '" << filename << "'" << std::ends;
+ throw cigma::Exception("WriteCellCentroids", stream.str());
+ }
+ }
+ else if (wt == FileWriter::VTK_FILE_WRITER)
+ {
+ throw cigma::Exception("WriteCellCentroids", "Not implemented for VTK files");
+ }
+ else if (wt == FileWriter::TEXT_FILE_WRITER)
+ {
+ int status = writer->writeDataset(location.c_str(), ¢roids[0], ncells, ndim);
+ if (status < 0)
+ {
+ std::ostringstream stream;
+ stream << "Could not write cell centroids to the text file '" << filename << "'" << std::ends;
+ throw cigma::Exception("WriteCellCentroids", stream.str());
+ }
+ }
+ else
+ {
+ std::ostringstream stream;
+ stream << "Could not write cell centroids to path '" << path << "'" << std::ends;
+ throw cigma::Exception("WriteCellCentroids", stream.str());
+ }
+}
+
+void cigma::WriteResiduals(const DataPath& path, shared_ptr<Residuals> residuals, bool normalize, bool overwrite)
+{
TRI_LOG_STR("cigma::WriteResiduals()");
TRI_LOG(path);
@@ -409,13 +475,30 @@
std::string filename = path.filename();
std::string location = path.location();
+ valarray<double> eps;
+ eps.resize(residuals->epsilon.size());
+
+ if (normalize)
+ {
+ TRI_LOG_STR(">> Writing out the normalized residuals");
+ if (residuals->mesh->volume.size() == 0)
+ {
+ residuals->mesh->computeCellVolumes();
+ }
+ eps = (residuals->epsilon) / (residuals->mesh->volume);
+ }
+ else
+ {
+ eps = residuals->epsilon;
+ }
+
shared_ptr<FileWriter> writer = FileWriter::New(filename, "a");
const FileWriter::WriterType wt = writer->getWriterType();
if (wt == FileWriter::HDF5_FILE_WRITER)
{
HDF5_Writer *h5 = static_cast<HDF5_Writer*>(&(*writer));
- // prepare writer to overwrite exsting datasets
+ // prepare writer to overwrite existing datasets
h5->setOverwriteMode(overwrite);
string eps_loc = location;
@@ -425,7 +508,8 @@
}
TRI_LOG(eps_loc);
- int status = h5->writeDataset(eps_loc.c_str(), &(residuals->epsilon[0]), residuals->mesh->n_cells(), 1);
+ //int status = h5->writeDataset(eps_loc.c_str(), &(residuals->epsilon[0]), residuals->mesh->n_cells(), 1);
+ int status = h5->writeDataset(eps_loc.c_str(), &eps[0], eps.size(), 1);
TRI_LOG(status);
if (status < 0)
{
@@ -436,13 +520,15 @@
throw cigma::Exception("WriteResiduals", stream.str());
}
- double L2 = residuals->L2();
+ double volume = residuals->vol;
+ double Linf = residuals->infinity_norm();
+ double L2 = residuals->L2();
double L2_rel = residuals->relative_L2();
- double Linf = residuals->infinity_norm();
- double volume = residuals->vol;
+ double L2_vol = L2/volume;
status = write_scalar_attribute<double>(h5->file, eps_loc.c_str(), "INFINITY_NORM", residuals->infinity_norm());
status = write_scalar_attribute<double>(h5->file, eps_loc.c_str(), "L2_NORM", L2);
status = write_scalar_attribute<double>(h5->file, eps_loc.c_str(), "L2_NORM_RELATIVE", L2_rel);
+ status = write_scalar_attribute<double>(h5->file, eps_loc.c_str(), "L2_NORM_BY_VOLUME", L2_vol);
status = write_scalar_attribute<double>(h5->file, eps_loc.c_str(), "VOLUME", volume);
}
else if (wt == FileWriter::VTK_FILE_WRITER)
@@ -475,7 +561,7 @@
points.init_value(0);
for (i = 0; i < npts; i++)
{
- for (j = 0; j < 3; j++)
+ for (j = 0; j < ndim; j++)
{
points(i,j) = nc->getPoint(i,j);
}
@@ -513,11 +599,13 @@
}
TRI_LOG_STR(">>> Writing residuals as vtk cell data");
- vtk_writer->_writeCellData(location.c_str(), &(residuals->epsilon[0]), residuals->mesh->n_cells(), 1);
+ //vtk_writer->_writeCellData(location.c_str(), &(residuals->epsilon[0]), residuals->mesh->n_cells(), 1);
+ vtk_writer->_writeCellData(location.c_str(), &eps[0], eps.size(), 1);
}
else if (wt == FileWriter::TEXT_FILE_WRITER)
{
- int status = writer->writeDataset("", &(residuals->epsilon[0]), residuals->mesh->n_cells(), 1);
+ //int status = writer->writeDataset("", &(residuals->epsilon[0]), residuals->mesh->n_cells(), 1);
+ int status = writer->writeDataset("", &eps[0], eps.size(), 1);
TRI_LOG(status);
if (status < 0)
{
Modified: cs/cigma/trunk/src/core_writers.h
===================================================================
--- cs/cigma/trunk/src/core_writers.h 2009-02-18 16:14:34 UTC (rev 14082)
+++ cs/cigma/trunk/src/core_writers.h 2009-02-18 16:14:36 UTC (rev 14083)
@@ -22,8 +22,9 @@
void WriteQuadrature(const DataPath& path, boost::shared_ptr<Quadrature> Q, bool overwrite=false);
void WriteQPoints(const DataPath& path, const cigma::array<double>& w, const cigma::array<double>& x, bool overwrite=false);
void WriteDofs(const DataPath& path, boost::shared_ptr<DofHandler> dofs);
- void WriteResiduals(const DataPath& path, boost::shared_ptr<Residuals> residuals, bool overwrite=false);
+ void WriteResiduals(const DataPath& path, boost::shared_ptr<Residuals> residuals, bool normalize, bool overwrite=false);
void WriteCellVolumes(const DataPath& path, const std::valarray<double>& volume, bool overwrite=false);
+ void WriteCellCentroids(const DataPath& path, boost::shared_ptr<MeshPart> M, bool overwrite=false);
}
#endif
More information about the CIG-COMMITS
mailing list