[cig-commits] r13831 - cs/cigma/trunk/src
luis at geodynamics.org
luis at geodynamics.org
Mon Jan 12 14:09:56 PST 2009
Author: luis
Date: 2009-01-12 14:09:56 -0800 (Mon, 12 Jan 2009)
New Revision: 13831
Modified:
cs/cigma/trunk/src/MeshPart.cpp
cs/cigma/trunk/src/MeshPart.h
cs/cigma/trunk/src/ProgressTimer.cpp
cs/cigma/trunk/src/ProgressTimer.h
cs/cigma/trunk/src/Residuals.cpp
cs/cigma/trunk/src/Residuals.h
cs/cigma/trunk/src/cli_compare_cmd.cpp
cs/cigma/trunk/src/core_compare_op.cpp
cs/cigma/trunk/src/core_writers.cpp
cs/cigma/trunk/src/vtk-residuals.cpp
Log:
Move mesh-related methods out of Residuals and into MeshPart class
* Calculation of MaxCellDiameter & Volume goes in MeshPart
* Use std::valarray<double> to store volume, residuals
* Return reference in ProgressTimer methods
Modified: cs/cigma/trunk/src/MeshPart.cpp
===================================================================
--- cs/cigma/trunk/src/MeshPart.cpp 2009-01-12 22:09:53 UTC (rev 13830)
+++ cs/cigma/trunk/src/MeshPart.cpp 2009-01-12 22:09:56 UTC (rev 13831)
@@ -53,15 +53,15 @@
this->cell_type = cell_type;
}
-void MeshPart::getCell(int e, Cell& cell) const
+shared_ptr<Cell> MeshPart::getCell() const
{
- this->getCellCoords(e, cell.globverts, cell.n_celldim());
+ shared_ptr<Cell> cell = Cell::New(cell_type);
+ return cell;
}
-void MeshPart::getBoundingBox(double *minpt, double *maxpt) const
+void MeshPart::getCell(int e, Cell& cell) const
{
- assert(coords);
- coords->getBoundingBox(minpt, maxpt);
+ this->getCellCoords(e, cell.globverts, cell.n_celldim());
}
void MeshPart::getCellCoords(int cellIndex, double *globalCoords, int dim) const
@@ -164,3 +164,93 @@
}
+void MeshPart::getBoundingBox(double *minpt, double *maxpt) const
+{
+ assert(coords);
+ coords->getBoundingBox(minpt, maxpt);
+}
+
+static double dist(double dx, double dy, double dz)
+{
+ return sqrt(dx*dx + dy*dy + dz*dz);
+}
+
+static double BoxDiameter(double bounds[6])
+{
+ const double x0 = bounds[0];
+ const double x1 = bounds[1];
+ const double y0 = bounds[2];
+ const double y1 = bounds[3];
+ const double z0 = bounds[4];
+ const double z1 = bounds[5];
+
+ assert(x0 <= x1);
+ assert(y0 <= y1);
+ assert(z0 <= z1);
+
+ return dist(x1-x0, y1-y0, z1-z0);
+}
+
+double MeshPart::getMaxCellDiameter() const
+{
+ double hmax = 0.0;
+ shared_ptr<Cell> cell = this->getCell();
+ if (cell)
+ {
+ int emax = -1;
+ double bounds[6];
+ int ncells = connect->n_cells();
+ for (int e = 0; e < ncells; e++)
+ {
+ this->getCell(e, *cell);
+ cell->getBounds(bounds);
+ double h = BoxDiameter(bounds);
+ if (h > hmax)
+ {
+ hmax = h;
+ emax = e;
+ }
+ }
+ }
+ return hmax;
+}
+
+double MeshPart::getVolume() const
+{
+ double vol = 0;
+ if (volume.size() != 0)
+ {
+ vol = volume.sum();
+ }
+ else
+ {
+ int nel = this->n_cells();
+ shared_ptr<Cell> cell = this->getCell();
+ for (int e = 0; e < nel; e++)
+ {
+ this->getCell(e, *cell);
+ vol += cell->volume();
+ }
+ }
+ return vol;
+}
+
+void MeshPart::computeCellVolumes()
+{
+ TRI_LOG_STR("MeshPart::computeCellVolumes()");
+ shared_ptr<Cell> cell = this->getCell();
+ if (cell)
+ {
+ int nel = this->n_cells();
+ volume.resize(nel);
+ for (int e = 0; e < nel; e++)
+ {
+ this->getCell(e, *cell);
+ volume[e] = cell->volume();
+ }
+ TRI_LOG(nel);
+ TRI_LOG(volume.size());
+ TRI_LOG(volume.sum());
+ }
+}
+
Modified: cs/cigma/trunk/src/MeshPart.h
===================================================================
--- cs/cigma/trunk/src/MeshPart.h 2009-01-12 22:09:53 UTC (rev 13830)
+++ cs/cigma/trunk/src/MeshPart.h 2009-01-12 22:09:56 UTC (rev 13831)
@@ -10,6 +10,7 @@
#include "Cell.h"
#include "DataPath.h"
#include "core_args.h"
+#include <valarray>
namespace cigma
{
@@ -25,14 +26,22 @@
int n_cells() const;
int n_dim() const;
- void getCell(int e, Cell &cell) const; // XXX
- void getBoundingBox(double *minpt, double *maxpt) const;
- void getCellCoords(int cellIndex, double *globalCoords, int dim) const;
Cell::type getCellType() const;
+ boost::shared_ptr<Cell> getCell() const;
+ void getCell(int e, Cell &cell) const;
+ void getCellCoords(int cellIndex, double *globalCoords, int dim) const;
+
bool findCell(double globalPoint[3], int *cellIndex);
bool findCell2(double globalPoint[3], double uvw[3], int *cellIndex);
+ void getBoundingBox(double *minpt, double *maxpt) const;
+
+ double getMaxCellDiameter() const;
+
+ void computeCellVolumes();
+ double getVolume() const;
+
/* direct setters */
void setNodeCoordinates(const boost::shared_ptr<NodeCoordinates>& nc);
void setElementBlock(const boost::shared_ptr<ElementBlock>& eb);
@@ -48,10 +57,11 @@
static boost::shared_ptr<MeshPart> New(const MeshInfo& mesh_info);
public:
+ Cell::type cell_type;
boost::shared_ptr<ElementBlock> connect;
boost::shared_ptr<NodeCoordinates> coords;
boost::shared_ptr<Locator> locator;
- Cell::type cell_type;
+ std::valarray<double> volume;
};
inline int cigma::MeshPart::n_cells() const { return connect->n_cells(); }
Modified: cs/cigma/trunk/src/ProgressTimer.cpp
===================================================================
--- cs/cigma/trunk/src/ProgressTimer.cpp 2009-01-12 22:09:53 UTC (rev 13830)
+++ cs/cigma/trunk/src/ProgressTimer.cpp 2009-01-12 22:09:56 UTC (rev 13831)
@@ -1,8 +1,6 @@
#include "ProgressTimer.h"
+using namespace std;
-
-// ----------------------------------------------------------------------------
-
ProgressTimer::ProgressTimer()
{
t_0 = 0;
@@ -14,36 +12,34 @@
{
}
-
-// ----------------------------------------------------------------------------
-
-void ProgressTimer::printHeader(std::ostream &os, const char *firstcol)
+string ProgressTimer::header(const char *firstcol)
{
- assert(firstcol != 0);
- os << firstcol << " rate mins eta total progress\n";
+ string hdr(firstcol);
+ hdr += " rate mins eta total progress\n";
+ return hdr;
}
-void ProgressTimer::start(int total)
+ProgressTimer& ProgressTimer::start(int tot)
{
time(&t_0);
- this->total = total;
+ total = tot;
+ return update(0);
}
-void ProgressTimer::update(int current)
+ProgressTimer& ProgressTimer::update(int cur)
{
time(&t_n);
- this->current = current;
+ current = cur;
elapsed_mins = (t_n - t_0) / 60.0;
mins_per_num = elapsed_mins / current;
num_per_sec = (1.0/60.0) / mins_per_num;
remaining_mins = (total - current) * mins_per_num;
total_mins = total * mins_per_num;
progress = 100 * elapsed_mins / total_mins;
+ return *this;
}
-// ----------------------------------------------------------------------------
-
-std::ostream &operator<<(std::ostream &os, const ProgressTimer &T)
+std::ostream& operator<<(std::ostream& os, const ProgressTimer& T)
{
os << T.current << " "
<< T.num_per_sec << " "
@@ -57,4 +53,3 @@
<< std::flush;
}
-// ----------------------------------------------------------------------------
Modified: cs/cigma/trunk/src/ProgressTimer.h
===================================================================
--- cs/cigma/trunk/src/ProgressTimer.h 2009-01-12 22:09:53 UTC (rev 13830)
+++ cs/cigma/trunk/src/ProgressTimer.h 2009-01-12 22:09:56 UTC (rev 13831)
@@ -1,13 +1,24 @@
#ifndef __PROGRESS_TIMER_H__
#define __PROGRESS_TIMER_H__
+#include <ctime>
#include <iostream>
+#include <string>
#include <cassert>
-#include <ctime>
class ProgressTimer
{
public:
+ ProgressTimer();
+ ~ProgressTimer();
+
+ std::string header(const char *firstcol);
+
+ ProgressTimer& start(int total);
+ ProgressTimer& update(int current);
+
+
+public:
int total;
int current;
time_t t_0, t_n;
@@ -17,18 +28,10 @@
double remaining_mins;
double total_mins;
double progress;
-
-public:
- ProgressTimer();
- ~ProgressTimer();
-
- void printHeader(std::ostream &os, const char *firstcol);
- void start(int total);
- void update(int current);
};
-std::ostream &operator<<(std::ostream &os, const ProgressTimer &T);
+std::ostream& operator<<(std::ostream& os, const ProgressTimer& T);
#endif
Modified: cs/cigma/trunk/src/Residuals.cpp
===================================================================
--- cs/cigma/trunk/src/Residuals.cpp 2009-01-12 22:09:53 UTC (rev 13830)
+++ cs/cigma/trunk/src/Residuals.cpp 2009-01-12 22:09:56 UTC (rev 13831)
@@ -1,210 +1,69 @@
#include "Residuals.h"
-#include <iostream>
#include <cassert>
#include <cmath>
+#include <limits>
using namespace std;
using namespace cigma;
using boost::shared_ptr;
-
Residuals::Residuals()
{
- nel = 0;
- epsilon = 0;
- volume = 0;
-
- global_error = 0;
- total_volume = 0;
- max_residual = 0;
- max_diameter = 0;
+ max2 = 0;
+ res2 = 0;
+ vol = 0;
}
Residuals::~Residuals()
{
- nel = 0;
- if (epsilon != 0)
- {
- delete [] epsilon;
- epsilon = 0;
- }
- if (volume != 0)
- {
- delete [] volume;
- volume = 0;
- }
}
-void Residuals::setMesh(shared_ptr<MeshPart> meshPart)
+void Residuals::setMesh(shared_ptr<MeshPart> mesh)
{
- this->meshPart = meshPart;
+ this->mesh = mesh;
- if (meshPart)
+ if (mesh)
{
- nel = meshPart->n_cells();
-
- if (epsilon != 0) { delete [] epsilon; }
- if (volume != 0) { delete [] volume; }
-
- epsilon = new double[nel];
- volume = new double[nel];
+ epsilon.resize(mesh->n_cells());
+ mesh->computeCellVolumes();
+ vol = mesh->getVolume();
}
}
-void Residuals::resetAccumulator()
+void Residuals::reset()
{
- assert(nel > 0);
- assert(meshPart);
- assert(epsilon != 0);
- assert(volume != 0);
-
- global_error = 0.0;
- max_residual = 0.0;
- total_volume = 0.0;
-
- for (int e = 0; e < nel; e++)
- {
- epsilon[e] = 0.0;
- volume[e] = 0.0;
- }
+ epsilon = numeric_limits<double>::max();
+ max2 = 0;
+ res2 = 0;
}
-void Residuals::update(int e, double cell_residual, double cell_volume)
+void Residuals::update(int e, double cell_res2, double cell_max2)
{
- assert(cell_residual >= 0);
- assert(cell_volume > 0);
+ // remember cell residual
+ epsilon[e] = sqrt(cell_res2);
- // remember residual on this cell
- epsilon[e] = cell_residual;
+ // update global error
+ res2 += cell_res2;
- // keep track of domain volume
- volume[e] = cell_volume;
- total_volume += cell_volume;
-
- // keep track of max residual
- if (cell_residual > max_residual)
+ // update max error
+ if (cell_max2 > max2)
{
- max_residual = cell_residual;
+ max2 = cell_max2;
}
-
- // XXX: generalize to global_error = norm_sum(global_error, cell_residual)
- global_error += cell_residual * cell_residual;
}
-void Residuals::computeMaxResidual()
+double Residuals::infinity_norm() const
{
- assert(meshPart);
- assert(nel > 0);
-
- double maxres = 0.0;
- for (int e = 0; e < nel; e++)
- {
- double res = epsilon[e];
- if (res > maxres)
- {
- maxres = res;
- }
- }
- max_residual = maxres;
+ return sqrt(max2);
}
-static double dist(double dx, double dy, double dz)
+double Residuals::L2() const
{
- return sqrt(dx*dx + dy*dy + dz*dz);
+ return sqrt(res2);
}
-static double BoxDiameter(double bounds[6])
+double Residuals::relative_L2() const
{
- const double x0 = bounds[0];
- const double x1 = bounds[1];
- const double y0 = bounds[2];
- const double y1 = bounds[3];
- const double z0 = bounds[4];
- const double z1 = bounds[5];
-
- assert(x0 <= x1);
- assert(y0 <= y1);
- assert(z0 <= z1);
-
- return dist(x1-x0, y1-y0, z1-z0);
+ return sqrt(res2/vol);
}
-void Residuals::computeMaxCellDiameter()
-{
- assert(meshPart);
- assert(meshPart->coords);
- assert(meshPart->connect);
- assert(meshPart->getCellType() != Cell::NONE);
- assert(nel > 0);
-
- double bounds[6];
- int npts, ncells;
- int e, emax;
- double h, hmax;
-
- shared_ptr<Cell> cell = Cell::New(meshPart->getCellType());
- hmax = 0.0;
- emax = -1;
-
- for (e = 0; e < nel; e++)
- {
- meshPart->getCell(e, *cell);
- cell->getBounds(bounds);
- h = BoxDiameter(bounds);
- if (h > hmax)
- {
- hmax = h;
- emax = e;
- }
- }
- max_diameter = hmax;
-}
-
-void Residuals::computeVolume()
-{
- assert(meshPart);
- assert(meshPart->coords);
- assert(meshPart->connect);
- assert(meshPart->getCellType() != Cell::NONE);
-
- shared_ptr<Cell> cell = Cell::New(meshPart->getCellType());
-
- double total = 0.0;
- for (int e = 0; e < nel; e++)
- {
- double vol = cell->volume();
- volume[e] = vol;
- total += vol;
- }
- total_volume = total;
-}
-
-double Residuals::norm_Linf() const
-{
- return max_residual;
-}
-
-double Residuals::norm_L2() const
-{
- // XXX: generalizes to norm_pow(global_error)
- return sqrt(global_error);
-}
-
-double Residuals::relative_squared_norm_L2() const
-{
- return global_error/total_volume;
-}
-
-double Residuals::relative_norm_L2() const
-{
- return sqrt(global_error/total_volume);
-}
-
-double Residuals::max_cell_diameter()
-{
- if (max_diameter <= 0)
- {
- this->computeMaxCellDiameter();
- }
- return max_diameter;
-}
Modified: cs/cigma/trunk/src/Residuals.h
===================================================================
--- cs/cigma/trunk/src/Residuals.h 2009-01-12 22:09:53 UTC (rev 13830)
+++ cs/cigma/trunk/src/Residuals.h 2009-01-12 22:09:56 UTC (rev 13831)
@@ -4,7 +4,6 @@
#include "MeshPart.h"
#include <boost/noncopyable.hpp>
#include <boost/shared_ptr.hpp>
-//#include <valarray>
namespace cigma
{
@@ -17,34 +16,23 @@
Residuals();
~Residuals();
- void setMesh(boost::shared_ptr<MeshPart> meshPart);
+ void setMesh(boost::shared_ptr<MeshPart> mesh);
- void resetAccumulator();
- void update(int e, double cell_residual, double cell_volume);
+ void reset();
+ void update(int e, double cell_res2, double cell_max2);
- void computeMaxResidual();
- void computeMaxCellDiameter();
- void computeVolume();
+ double infinity_norm() const;
+ double L2() const;
+ double relative_L2() const;
- double norm_Linf() const;
- double norm_L2() const;
- double relative_squared_norm_L2() const;
- double relative_norm_L2() const;
- double max_cell_diameter();
-
public:
- boost::shared_ptr<MeshPart> meshPart;
+ boost::shared_ptr<MeshPart> mesh;
+ std::valarray<double> epsilon;
- int nel;
- double *epsilon;
- double *volume;
-
- double global_error;
- double total_volume;
- double max_residual;
- double max_diameter;
-
+ double max2;
+ double res2;
+ double vol;
};
Modified: cs/cigma/trunk/src/cli_compare_cmd.cpp
===================================================================
--- cs/cigma/trunk/src/cli_compare_cmd.cpp 2009-01-12 22:09:53 UTC (rev 13830)
+++ cs/cigma/trunk/src/cli_compare_cmd.cpp 2009-01-12 22:09:56 UTC (rev 13831)
@@ -158,11 +158,13 @@
int status = op.run();
- double L2 = op.residuals->norm_L2();
- double Linf = op.residuals->norm_Linf();
- double relative_L2 = op.residuals->relative_norm_L2();
- double volume = op.residuals->total_volume;
- double h = op.residuals->max_cell_diameter();
+ double Linf = op.residuals->infinity_norm();
+ double L2 = op.residuals->L2();
+ double relative_L2 = op.residuals->relative_L2();
+ double volume = op.residuals->mesh->getVolume();
+ double h = op.residuals->mesh->getMaxCellDiameter();
+ //double h1 = op.first->getMaxCellDiameter();
+ //double h2 = op.second->getMaxCellDiameter();
if (L2 > threshold)
{
@@ -180,9 +182,8 @@
cout << "Summary of comparison: " << endl;
cout << " L2 = " << L2 << endl;
cout << " Linf = " << Linf << endl;
- cout << " h = " << h << endl;
+ cout << " sqrt((L2)^2/volume) = " << relative_L2 << endl;
cout << " volume of domain = " << volume << endl;
- cout << " sqrt((L2)^2/volume) = " << relative_L2 << endl;
cout << endl;
// write out data
Modified: cs/cigma/trunk/src/core_compare_op.cpp
===================================================================
--- cs/cigma/trunk/src/core_compare_op.cpp 2009-01-12 22:09:53 UTC (rev 13830)
+++ cs/cigma/trunk/src/core_compare_op.cpp 2009-01-12 22:09:56 UTC (rev 13831)
@@ -182,37 +182,35 @@
const int nel = mesh->n_cells();
const int nq = quadrature->n_points();
+
const int ndim = first->n_dim();
const int nval = first->n_rank();
+ double *qpts = new double[nq * 3];
double *vals1 = new double[nq * nval];
double *vals2 = new double[nq * nval];
- double *qpts = new double[nq * 3]; // local quadrature points in global coordinates
-
- residuals = shared_ptr<Residuals>(new Residuals());
+ residuals = shared_ptr<Residuals>(new Residuals);
residuals->setMesh(domain->mesh);
- residuals->resetAccumulator();
+ residuals->reset();
-
if (verbose)
{
- // XXX: timer_start(nel, "elts");
+ cout << timer.header("elts");
cout << setprecision(5);
- timer.printHeader(cout, "elts");
- timer.start(nel);
- timer.update(0);
- cout << timer;
+ cout << timer.start(nel);
}
int e,q,j;
for (e = 0; e < nel; e++)
{
- // load e-th cell from mesh (XXX: use iterator)
- mesh->getCell(e, *cell); // load cell->globverts from mesh
- fe->update_jxw(); // update jxw factor
+ // load e-th cell from mesh
+ mesh->getCell(e, *cell);
- // map fe (local quadrature points) to global coords
+ // update jxw factors
+ fe->update_jxw();
+
+ // map quadrature points to global coords
for (q = 0; q < nq; q++)
{
double uvw[3] = {0,0,0};
@@ -224,7 +222,7 @@
cell->uvw2xyz(uvw, xyz);
}
- // evaluate both functions at q-th point (in global coordinates)
+ // evaluate both functions at q-th point
for (q = 0; q < nq; q++)
{
double *globalPoint = &qpts[q*3];
@@ -232,39 +230,45 @@
second->eval(globalPoint, &vals2[q*nval]);
}
- // calculate norm on local cell
- double err2 = 0;
+ // calculate norms on local cell
+ double res2 = 0;
+ double max2 = 0;
for (q = 0; q < nq; q++)
{
+ double err2 = 0;
+
for (j = 0; j < nval; j++)
{
int k = q*nval + j;
double f_qj = vals1[k];
double g_qj = vals2[k];
- err2 += fe->jxw[q] * SQR(f_qj - g_qj);
+ double d = f_qj - g_qj;
+ double d2 = d*d;
+ err2 += d2;
}
+
+ res2 += fe->jxw[q] * err2;
+
+ if (err2 > max2)
+ {
+ max2 = err2;
+ }
}
- // calculate volume of cell
- double vol = cell->volume();
-
// update residuals
- residuals->update(e, sqrt(err2), vol);
+ residuals->update(e, res2, max2);
// update timer
if (verbose && ((e+1) % freq == 0))
{
- // XXX: timer_update(e)
- timer.update(e+1);
- cout << timer;
+ cout << timer.update(e+1);
}
}
+ // final update to timer
if (verbose)
{
- // XXX: timer_end()
- timer.update(nel);
- cout << timer << endl;
+ cout << timer.update(nel) << endl;
}
delete [] vals1;
Modified: cs/cigma/trunk/src/core_writers.cpp
===================================================================
--- cs/cigma/trunk/src/core_writers.cpp 2009-01-12 22:09:53 UTC (rev 13830)
+++ cs/cigma/trunk/src/core_writers.cpp 2009-01-12 22:09:56 UTC (rev 13831)
@@ -362,9 +362,6 @@
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();
@@ -384,7 +381,7 @@
}
TRI_LOG(eps_loc);
- int status = h5->writeDataset(eps_loc.c_str(), residuals->epsilon, residuals->nel, 1);
+ int status = h5->writeDataset(eps_loc.c_str(), &(residuals->epsilon[0]), residuals->mesh->n_cells(), 1);
TRI_LOG(status);
if (status < 0)
{
@@ -395,11 +392,14 @@
throw cigma::Exception("WriteResiduals", stream.str());
}
- status = write_scalar_attribute<double>(h5->file, eps_loc.c_str(), "L2_NORM", residuals->norm_L2());
- status = write_scalar_attribute<double>(h5->file, eps_loc.c_str(), "RELATIVE_L2", residuals->relative_norm_L2());
- status = write_scalar_attribute<double>(h5->file, eps_loc.c_str(), "VOLUME", residuals->total_volume);
- status = write_scalar_attribute<double>(h5->file, eps_loc.c_str(), "LINF_NORM", residuals->norm_Linf());
- status = write_scalar_attribute<double>(h5->file, eps_loc.c_str(), "MAX_CELL_DIAMETER", residuals->max_cell_diameter());
+ double L2 = residuals->L2();
+ double L2_rel = residuals->relative_L2();
+ double Linf = residuals->infinity_norm();
+ double volume = residuals->vol;
+ 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(), "VOLUME", volume);
}
else if (wt == FileWriter::VTK_FILE_WRITER)
{
@@ -412,7 +412,7 @@
int status;
- if (residuals->meshPart)
+ if (residuals->mesh)
{
int i,j;
if (location == "")
@@ -422,11 +422,11 @@
// coordinates first
TRI_LOG_STR(">>> Writing vtk points");
- const int npts = residuals->meshPart->coords->n_points();
- const int ndim = residuals->meshPart->coords->n_dim();
+ const int npts = residuals->mesh->coords->n_points();
+ const int ndim = residuals->mesh->coords->n_dim();
TRI_LOG(npts);
TRI_LOG(ndim);
- NodeCoordinates *nc = static_cast<NodeCoordinates*>(&(*(residuals->meshPart->coords)));
+ NodeCoordinates *nc = static_cast<NodeCoordinates*>(&(*(residuals->mesh->coords)));
cigma::array<double> points(npts, 3);
points.init_value(0);
for (i = 0; i < npts; i++)
@@ -441,11 +441,11 @@
// cells next..
TRI_LOG_STR(">>> Writing vtk cells");
- const int nel = residuals->meshPart->connect->n_cells();
- const int nno = residuals->meshPart->connect->n_dofs();
+ const int nel = residuals->mesh->connect->n_cells();
+ const int nno = residuals->mesh->connect->n_dofs();
TRI_LOG(nel);
TRI_LOG(nno);
- ElementBlock *eb = static_cast<ElementBlock*>(&(*(residuals->meshPart->connect)));
+ ElementBlock *eb = static_cast<ElementBlock*>(&(*(residuals->mesh->connect)));
cigma::array<int> cells(nel, nno);
for (i = 0; i < nel; i++)
{
@@ -459,7 +459,7 @@
// finally, cell types
TRI_LOG_STR(">>> Writing vtk cell types");
- int vtkCellType = vtkcelltype(residuals->meshPart->getCellType());
+ int vtkCellType = vtkcelltype(residuals->mesh->getCellType());
TRI_LOG(vtkCellType);
vtk_writer->_writeCellTypes(nel, vtkCellType);
}
@@ -469,11 +469,11 @@
}
TRI_LOG_STR(">>> Writing residuals as vtk cell data");
- vtk_writer->_writeCellData(location.c_str(), residuals->epsilon, residuals->nel, 1);
+ vtk_writer->_writeCellData(location.c_str(), &(residuals->epsilon[0]), residuals->mesh->n_cells(), 1);
}
else if (wt == FileWriter::TEXT_FILE_WRITER)
{
- int status = writer->writeDataset("", residuals->epsilon, residuals->nel, 1);
+ int status = writer->writeDataset("", &(residuals->epsilon[0]), residuals->mesh->n_cells(), 1);
TRI_LOG(status);
if (status < 0)
{
Modified: cs/cigma/trunk/src/vtk-residuals.cpp
===================================================================
--- cs/cigma/trunk/src/vtk-residuals.cpp 2009-01-12 22:09:53 UTC (rev 13830)
+++ cs/cigma/trunk/src/vtk-residuals.cpp 2009-01-12 22:09:56 UTC (rev 13831)
@@ -1,6 +1,7 @@
#include <cstdlib>
#include <iostream>
#include <string>
+#include <valarray>
#include <boost/program_options.hpp>
#include <boost/filesystem.hpp>
@@ -132,13 +133,8 @@
TRI_LOG_STR(">> Configuring residuals object");
shared_ptr<Residuals> residuals(new Residuals);
- residuals->meshPart = mesh;
- residuals->nel = mesh->n_cells();
+ residuals->setMesh(mesh);
- residuals->volume = new double[residuals->nel];
- residuals->computeMaxCellDiameter();
- residuals->computeVolume();
-
TRI_LOG_STR(">> Reading residuals from HDF5 file");
TRI_LOG(in);
cigma::array<double>* epsilon = 0;
@@ -174,11 +170,9 @@
delete epsilon;
return 4;
}
- residuals->epsilon = epsilon->_data;
- epsilon->_data = 0;
+ valarray<double> eps(epsilon->_data, epsilon->size());
+ residuals->epsilon = eps;
delete epsilon;
- epsilon = 0;
- residuals->computeMaxResidual();
TRI_LOG_STR(">> Writing residuals to VTK file");
cout << "Writing residuals to VTK file '" << out.filename() << "'" << endl;
More information about the CIG-COMMITS
mailing list