[cig-commits] r13520 - cs/cigma/trunk/src
luis at geodynamics.org
luis at geodynamics.org
Tue Dec 9 18:13:11 PST 2008
Author: luis
Date: 2008-12-09 18:13:10 -0800 (Tue, 09 Dec 2008)
New Revision: 13520
Added:
cs/cigma/trunk/src/core_args.cpp
cs/cigma/trunk/src/core_args.h
Modified:
cs/cigma/trunk/src/core_base_op.cpp
cs/cigma/trunk/src/core_base_op.h
cs/cigma/trunk/src/core_compare_op.cpp
cs/cigma/trunk/src/core_compare_op.h
cs/cigma/trunk/src/core_eval_op.cpp
cs/cigma/trunk/src/core_eval_op.h
cs/cigma/trunk/src/core_extract_op.cpp
cs/cigma/trunk/src/core_extract_op.h
cs/cigma/trunk/src/core_list_op.cpp
cs/cigma/trunk/src/core_list_op.h
cs/cigma/trunk/src/core_readers.cpp
cs/cigma/trunk/src/core_readers.h
cs/cigma/trunk/src/core_writers.cpp
cs/cigma/trunk/src/core_writers.h
Log:
Updates to core objects & functions
Added: cs/cigma/trunk/src/core_args.cpp
===================================================================
--- cs/cigma/trunk/src/core_args.cpp (rev 0)
+++ cs/cigma/trunk/src/core_args.cpp 2008-12-10 02:13:10 UTC (rev 13520)
@@ -0,0 +1,144 @@
+#include "core_args.h"
+#include "MeshPart.h"
+
+using namespace cigma;
+
+// ----------------------------------------------------------------------------
+
+MeshInfo::MeshInfo() : use_locator(false) {}
+MeshInfo::~MeshInfo() {}
+
+MeshInfo::MeshInfo(const MeshInfo& other)
+{
+ p_mesh = other.p_mesh;
+ p_nc = other.p_nc;
+ p_eb = other.p_eb;
+ use_locator = other.use_locator;
+ cell_type_name = other.cell_type_name;
+}
+
+MeshInfo& MeshInfo::operator=(const MeshInfo& other)
+{
+ p_mesh = other.p_mesh;
+ p_nc = other.p_nc;
+ p_eb = other.p_eb;
+ use_locator = other.use_locator;
+ cell_type_name = other.cell_type_name;
+}
+
+DataPath MeshInfo::get_nc_path() const
+{
+ // allow override
+ if (!p_nc.empty())
+ {
+ return p_nc;
+ }
+
+ DataPath path;
+ if (!p_mesh.empty())
+ {
+ path.set_filename(p_mesh.filename());
+ path.set_location(p_mesh.location() + "/coordinates");
+ }
+ return path;
+}
+
+DataPath MeshInfo::get_eb_path() const
+{
+ // allow override
+ if (!p_eb.empty())
+ {
+ return p_eb;
+ }
+ DataPath path;
+ if (!p_mesh.empty())
+ {
+ path.set_filename(p_mesh.filename());
+ path.set_location(p_mesh.location() + "/connectivity");
+ }
+ return path;
+}
+
+// ----------------------------------------------------------------------------
+
+QuadratureInfo::QuadratureInfo() {}
+QuadratureInfo::~QuadratureInfo() {}
+
+QuadratureInfo::QuadratureInfo(const QuadratureInfo& other)
+{
+ p_quadrature = other.p_quadrature;
+ p_weights = other.p_weights;
+ p_points = other.p_points;
+ cell_type_name = other.cell_type_name;
+}
+
+QuadratureInfo& QuadratureInfo::operator=(const QuadratureInfo& other)
+{
+ p_quadrature = other.p_quadrature;
+ p_weights = other.p_weights;
+ p_points = other.p_points;
+ cell_type_name = other.cell_type_name;
+}
+
+bool QuadratureInfo::single_source() const
+{
+ const bool q = !p_quadrature.empty();
+ const bool no_w = p_weights.empty();
+ const bool no_x = p_points.empty();
+ return q && (no_w || no_x);
+}
+
+bool QuadratureInfo::double_source() const
+{
+ const bool w = !p_weights.empty();
+ const bool x = !p_points.empty();
+ const bool no_q = p_quadrature.empty();
+ return no_q && (w && x);
+}
+
+QuadratureInfo::operator bool() const
+{
+ const bool s = this->single_source();
+ const bool d = this->double_source();
+ return (s && !d) || (!s && d);
+}
+
+// ----------------------------------------------------------------------------
+
+FE_Info::FE_Info() {}
+FE_Info::~FE_Info() {}
+
+FE_Info::FE_Info(const FE_Info& other)
+{
+ q_info = other.q_info;
+ p_fe_basis = other.p_fe_basis;
+ p_fe_basis_jet = other.p_fe_basis_jet;
+}
+
+FE_Info& FE_Info::operator=(const FE_Info& other)
+{
+ q_info = other.q_info;
+ p_fe_basis = other.p_fe_basis;
+ p_fe_basis_jet = other.p_fe_basis_jet;
+}
+
+// ----------------------------------------------------------------------------
+
+FieldInfo::FieldInfo() {}
+FieldInfo::~FieldInfo() {}
+
+FieldInfo::FieldInfo(const FieldInfo& other)
+{
+ p_field = other.p_field;
+ fe_info = other.fe_info;
+ mesh_info = other.mesh_info;
+}
+
+FieldInfo& FieldInfo::operator=(const FieldInfo& other)
+{
+ p_field = other.p_field;
+ fe_info = other.fe_info;
+ mesh_info = other.mesh_info;
+}
+
+// ----------------------------------------------------------------------------
Added: cs/cigma/trunk/src/core_args.h
===================================================================
--- cs/cigma/trunk/src/core_args.h (rev 0)
+++ cs/cigma/trunk/src/core_args.h 2008-12-10 02:13:10 UTC (rev 13520)
@@ -0,0 +1,82 @@
+#ifndef __CIGMA_CORE_ARGS_H__
+#define __CIGMA_CORE_ARGS_H__
+
+#include "Exception.h"
+#include "DataPath.h"
+#include "Cell.h"
+//#include "io_file_reader.h"
+
+namespace cigma
+{
+
+ struct MeshInfo
+ {
+ DataPath p_mesh;
+ DataPath p_nc;
+ DataPath p_eb;
+ bool use_locator;
+ std::string cell_type_name;
+
+ MeshInfo();
+ ~MeshInfo();
+ MeshInfo(const MeshInfo& other);
+ MeshInfo& operator=(const MeshInfo& other);
+ DataPath get_nc_path() const;
+ DataPath get_eb_path() const;
+ //bool single_file() const;
+ //FileReader::ReaderType get_nc_reader_type() const;
+ //FileReader::ReaderType get_eb_reader_type() const;
+ };
+
+
+ struct QuadratureInfo
+ {
+ DataPath p_quadrature;
+ DataPath p_weights;
+ DataPath p_points;
+ std::string cell_type_name;
+
+ QuadratureInfo();
+ ~QuadratureInfo();
+ QuadratureInfo(const QuadratureInfo& other);
+ QuadratureInfo& operator=(const QuadratureInfo& other);
+ bool single_source() const;
+ bool double_source() const;
+ operator bool() const;
+ //DataPath get_weights_path() const;
+ //DataPath get_points_path() const;
+ //FileReader::ReaderType get_weights_reader_type() const;
+ //FileReader::ReaderType get_points_reader_type() const;
+ };
+
+
+ struct FE_Info
+ {
+ QuadratureInfo q_info;
+ DataPath p_fe_basis;
+ DataPath p_fe_basis_jet;
+
+ FE_Info();
+ ~FE_Info();
+ FE_Info(const FE_Info& other);
+ FE_Info& operator=(const FE_Info& other);
+ };
+
+
+ struct FieldInfo
+ {
+ DataPath p_field;
+ FE_Info fe_info;
+ MeshInfo mesh_info;
+
+ FieldInfo();
+ ~FieldInfo();
+ FieldInfo(const FieldInfo& other);
+ FieldInfo& operator=(const FieldInfo& other);
+ };
+
+
+}
+
+
+#endif
Modified: cs/cigma/trunk/src/core_base_op.cpp
===================================================================
--- cs/cigma/trunk/src/core_base_op.cpp 2008-12-10 02:13:09 UTC (rev 13519)
+++ cs/cigma/trunk/src/core_base_op.cpp 2008-12-10 02:13:10 UTC (rev 13520)
@@ -4,6 +4,7 @@
BaseOp::BaseOp()
{
verbose = false;
+ freq = 0;
}
BaseOp::~BaseOp()
Modified: cs/cigma/trunk/src/core_base_op.h
===================================================================
--- cs/cigma/trunk/src/core_base_op.h 2008-12-10 02:13:09 UTC (rev 13519)
+++ cs/cigma/trunk/src/core_base_op.h 2008-12-10 02:13:10 UTC (rev 13520)
@@ -13,12 +13,19 @@
class cigma::BaseOp : private boost::noncopyable
{
public:
+
BaseOp();
virtual ~BaseOp();
+ virtual void configure() = 0;
+ virtual int run() = 0;
+
public:
+
bool verbose;
ProgressTimer timer;
+ int freq;
+
};
#endif
Modified: cs/cigma/trunk/src/core_compare_op.cpp
===================================================================
--- cs/cigma/trunk/src/core_compare_op.cpp 2008-12-10 02:13:09 UTC (rev 13519)
+++ cs/cigma/trunk/src/core_compare_op.cpp 2008-12-10 02:13:10 UTC (rev 13520)
@@ -1,9 +1,19 @@
-#include <iostream>
#include "core_compare_op.h"
+#include "tri_logger.hpp"
+#include "Numeric.h"
+#include "fn_zero.h"
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+#include <cassert>
+#include <boost/shared_ptr.hpp>
+
using namespace std;
using namespace cigma;
+using boost::shared_ptr;
+
CompareOp::CompareOp()
{
}
@@ -12,8 +22,194 @@
{
}
-void CompareOp::run()
+void CompareOp::configure()
{
- cout << "CompareOp::run()" << endl;
+ TRI_LOG_STR("CompareOp::configure()");
+
+ if (!first)
+ {
+ if (first_name != "")
+ {
+ first = Function::New(first_name.c_str());
+ }
+ else
+ {
+ first = ReadField(first_info);
+ }
+ }
+ if (!first)
+ {
+ string msg("First field not specified");
+ throw cigma::Exception("CompareOp::configure", msg);
+ }
+
+
+ if (!second)
+ {
+ if (second_name != "")
+ {
+ second = Function::New(second_name.c_str());
+ }
+ else
+ {
+ second = ReadField(second_info);
+ }
+ }
+ if (!second)
+ {
+ string msg("Second field not specified");
+ throw cigma::Exception("CompareOp::configure", msg);
+ }
+
+ if (!domain)
+ {
+ domain = ReadField(domain_info);
+ }
+ assert(domain);
+ assert(domain->mesh);
+ assert(domain->mesh->coords);
+ assert(domain->mesh->connect);
+ assert(domain->fe);
+ assert(domain->fe->cell);
+ assert(domain->fe->quadrature);
+
+ if (first_name == "zero")
+ {
+ ZeroFunction *fn = static_cast<ZeroFunction*>(&(*first));
+ fn->setShape(second->n_dim(), second->n_rank());
+ }
+
+ if (second_name == "zero")
+ {
+ ZeroFunction *fn = static_cast<ZeroFunction*>(&(*second));
+ fn->setShape(first->n_dim(), first->n_rank());
+ }
+
+ if (first->n_dim() != second->n_dim())
+ {
+ string msg("Dimension mismatch for domains of first and second functions");
+ throw cigma::Exception("CompareOp::configure", msg);
+ }
+
+ if (first->n_rank() != second->n_rank())
+ {
+ string msg("Rank mismatch for first and second functions");
+ throw cigma::Exception("CompareOp::configure", msg);
+ }
+
}
+int CompareOp::run()
+{
+ TRI_LOG_STR("CompareOp::run()");
+
+ assert(first);
+ assert(second);
+ assert(domain);
+ assert(domain->mesh);
+ assert(domain->mesh->coords);
+ assert(domain->mesh->connect);
+ assert(domain->fe);
+ assert(domain->fe->cell);
+ assert(domain->fe->quadrature);
+
+ // convenient references
+ shared_ptr<MeshPart> mesh = domain->mesh;
+ shared_ptr<FE> fe = domain->fe;
+ shared_ptr<Cell> cell = fe->cell;
+ shared_ptr<Quadrature> quadrature = fe->quadrature;
+
+ 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 *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->setMesh(domain->mesh);
+ residuals->resetAccumulator();
+
+
+ if (verbose)
+ {
+ // XXX: timer_start(nel, "elts");
+ cout << setprecision(5);
+ timer.printHeader(cout, "elts");
+ timer.start(nel);
+ timer.update(0);
+ cout << timer;
+ }
+
+ 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
+
+ // map fe (local quadrature points) to global coords
+ for (q = 0; q < nq; q++)
+ {
+ double uvw[3] = {0,0,0};
+ for (j = 0; j < ndim; j++)
+ {
+ uvw[j] = quadrature->getPoint(q,j);
+ }
+ double *xyz = &qpts[q*3];
+ cell->uvw2xyz(uvw, xyz);
+ }
+
+ // evaluate both functions at q-th point (in global coordinates)
+ for (q = 0; q < nq; q++)
+ {
+ double *globalPoint = &qpts[q*3];
+ first->eval(globalPoint, &vals1[q*nval]);
+ second->eval(globalPoint, &vals2[q*nval]);
+ }
+
+ // calculate norm on local cell
+ double err2 = 0;
+ for (q = 0; q < nq; q++)
+ {
+ 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);
+ }
+ }
+
+ // calculate volume of cell
+ double vol = cell->volume();
+
+ // update residuals
+ residuals->update(e, sqrt(err2), vol);
+
+ // update timer
+ if (verbose && ((e+1) % freq == 0))
+ {
+ // XXX: timer_update(e)
+ timer.update(e+1);
+ cout << timer;
+ }
+ }
+
+ if (verbose)
+ {
+ // XXX: timer_end()
+ timer.update(nel);
+ cout << timer << endl;
+ }
+
+ delete [] vals1;
+ delete [] vals2;
+ delete [] qpts;
+
+ return 0;
+}
+
Modified: cs/cigma/trunk/src/core_compare_op.h
===================================================================
--- cs/cigma/trunk/src/core_compare_op.h 2008-12-10 02:13:09 UTC (rev 13519)
+++ cs/cigma/trunk/src/core_compare_op.h 2008-12-10 02:13:10 UTC (rev 13520)
@@ -1,10 +1,11 @@
#ifndef __CIGMA_COMPARE_OP_H__
#define __CIGMA_COMPARE_OP_H__
+#include <boost/shared_ptr.hpp>
#include "core_base_op.h"
+#include "core_readers.h"
#include "Function.h"
-#include "MeshPart.h"
-#include "Quadrature.h"
+#include "Field.h"
#include "Residuals.h"
namespace cigma
@@ -12,27 +13,32 @@
class CompareOp;
}
+
class cigma::CompareOp : public cigma::BaseOp
{
public:
+
CompareOp();
~CompareOp();
- /*
- void setMesh(..);
- void setFirstFunction(...);
- void setSecondFunction(...);
- void computeResiduals();
- */
+ void configure();
+ int run();
- void run();
+public:
-public:
+ std::string first_name;
+ std::string second_name;
+
+ FieldInfo first_info;
+ FieldInfo second_info;
+ FieldInfo domain_info;
+
boost::shared_ptr<Function> first;
boost::shared_ptr<Function> second;
- boost::shared_ptr<MeshPart> mesh;
- boost::shared_ptr<Quadrature> quadrature;
+ boost::shared_ptr<Field> domain; // integration region (container for mesh & fe objects, w/ respective coords as the dof_values)
+
boost::shared_ptr<Residuals> residuals;
+
};
#endif
Modified: cs/cigma/trunk/src/core_eval_op.cpp
===================================================================
--- cs/cigma/trunk/src/core_eval_op.cpp 2008-12-10 02:13:09 UTC (rev 13519)
+++ cs/cigma/trunk/src/core_eval_op.cpp 2008-12-10 02:13:10 UTC (rev 13520)
@@ -1,8 +1,12 @@
+#include "core_eval_op.h"
+#include "core_readers.h"
+#include "tri_logger.hpp"
+#include <boost/shared_ptr.hpp>
#include <iostream>
-#include "core_eval_op.h"
using namespace std;
using namespace cigma;
+using boost::shared_ptr;
EvalOp::EvalOp()
{
@@ -12,8 +16,81 @@
{
}
-void EvalOp::run()
+
+void EvalOp::configure()
{
- cout << "EvalOp::run()" << endl;
+ TRI_LOG_STR("EvalOp::configure()");
+
+ if (!function)
+ {
+ assert(fn_name != ""); // XXX: throw cigma::Exception instead
+ function = Function::New(fn_name.c_str());
+ }
+ if (!function)
+ {
+ string msg("No function specified");
+ throw cigma::Exception("EvalOp::configure", msg);
+ }
+
+ if (!coords)
+ {
+ coords = ReadNodeCoordinates(nc_path);
+ }
+ if (!coords)
+ {
+ string msg("No points specified");
+ throw cigma::Exception("EvalOp::configure", msg);
+
+ }
+
+ if (coords->n_dim() != function->n_dim())
+ {
+ string msg("Dimension mismatch for function domain and coordinates");
+ throw cigma::Exception("EvalOp::configure", msg);
+ }
+
}
+
+int EvalOp::run()
+{
+ TRI_LOG_STR("EvalOp::run()");
+
+ assert(function);
+ assert(coords);
+
+ const int nno = coords->n_points();
+ const int dim = coords->n_dim();
+ const int rank = function->n_rank();
+
+ // allocate enough space for the values
+ values.reinit(nno, rank);
+
+
+ // temporary array
+ bool *err = new bool[nno];
+
+ int i,j,k;
+ for (i = 0; i < nno; i++)
+ {
+ double pt[dim];
+ double val[rank];
+
+ for (j = 0; j < dim; j++)
+ {
+ pt[j] = coords->getPoint(i,j);
+ }
+
+ err[i] = function->eval(pt, val);
+
+ for (k = 0; k < rank; k++)
+ {
+ values(i,k) = val[k];
+ }
+ }
+
+ delete [] err;
+
+ return 0;
+}
+
Modified: cs/cigma/trunk/src/core_eval_op.h
===================================================================
--- cs/cigma/trunk/src/core_eval_op.h 2008-12-10 02:13:09 UTC (rev 13519)
+++ cs/cigma/trunk/src/core_eval_op.h 2008-12-10 02:13:10 UTC (rev 13520)
@@ -2,8 +2,10 @@
#define __CIGMA_EVAL_OP_H__
#include "core_base_op.h"
+#include "DataPath.h"
#include "Function.h"
#include "NodeCoordinates.h"
+#include "core_array.h"
namespace cigma
{
@@ -13,15 +15,23 @@
class cigma::EvalOp : public cigma::BaseOp
{
public:
+
EvalOp();
~EvalOp();
- void run();
+ void configure();
+ int run();
public:
+
+ std::string fn_name;
+ DataPath nc_path;
+
boost::shared_ptr<Function> function;
boost::shared_ptr<NodeCoordinates> coords;
- //boost::shared_ptr<cigma::array> points; // XXX: array of points (2d or 3d)
+
+ cigma::array<double> values; // array of points (2d or 3d)
+
};
#endif
Modified: cs/cigma/trunk/src/core_extract_op.cpp
===================================================================
--- cs/cigma/trunk/src/core_extract_op.cpp 2008-12-10 02:13:09 UTC (rev 13519)
+++ cs/cigma/trunk/src/core_extract_op.cpp 2008-12-10 02:13:10 UTC (rev 13520)
@@ -1,9 +1,14 @@
+#include "core_extract_op.h"
+#include "tri_logger.hpp"
+#include <boost/shared_ptr.hpp>
#include <iostream>
-#include "core_extract_op.h"
+#include <cassert>
using namespace std;
using namespace cigma;
+using boost::shared_ptr;
+
ExtractOp::ExtractOp()
{
}
@@ -12,8 +17,127 @@
{
}
-void ExtractOp::run()
+void ExtractOp::configure()
{
- cout << "ExtractOp::run()" << endl;
+ TRI_LOG_STR("ExtractOp::configure()");
+
+ if (!mesh)
+ {
+ mesh = ReadMeshPart(mesh_info);
+ }
+
+ assert(mesh);
+ assert(mesh->coords);
+ assert(mesh->connect);
+
+ if (!quadrature)
+ {
+ quadrature = ReadQuadrature(quadrature_info);
+ }
+
+ assert(quadrature);
+
}
+int ExtractOp::run()
+{
+ TRI_LOG_STR("ExtractOp::run()");
+
+ assert(mesh);
+ assert(mesh->coords);
+ assert(mesh->connect);
+ assert(quadrature);
+
+ const int nel = mesh->connect->n_cells();
+ const int ndofs = mesh->connect->n_dofs();
+
+ const int nno = mesh->coords->n_points();
+ const int ndim = mesh->coords->n_dim();
+
+ const int nq = quadrature->n_points();
+ const int q_ndim = quadrature->n_dim();
+
+ assert(ndim == q_ndim);
+
+ // allocate enough space for points & weights
+ weights.reinit(nel*nq, 1);
+ points.reinit(nel*nq, ndim);
+
+ //weights.init_value(1);
+ //points.init_value(2);
+
+ shared_ptr<FE> fe(new FE);
+ fe->init(quadrature);
+ assert(fe->cell);
+ assert(fe->quadrature);
+ assert(fe->jxw != 0);
+ assert(fe->basis_tab != 0);
+
+ int i,j;
+ double *global_qpts = new double[nq * 3];
+ for (i = 0; i < nq*3; i++) global_qpts[i] = 0;
+
+ int e,q;
+ for (e = 0; e < nel; e++)
+ {
+ // load e-th cell from mesh
+ mesh->getCell(e, *(fe->cell));
+ fe->update_jxw();
+
+ // map fe (local qpts) to global coords
+ for (q = 0; q < nq; q++)
+ {
+ double uvw[3] = {0,0,0};
+ for (j = 0; j < ndim; j++)
+ {
+ uvw[j] = quadrature->getPoint(q,j);
+ }
+ double *xyz = &global_qpts[q*3];
+ fe->cell->uvw2xyz(uvw, xyz);
+ }
+
+ // copy global coords to output arrays
+ for (q = 0; q < nq; q++)
+ {
+ weights(e*nq + q, 0) = fe->jxw[q];
+ for (j = 0; j < q_ndim; j++)
+ points(e*nq + q, j) = global_qpts[q*3 + j];
+ }
+ }
+
+ delete [] global_qpts;
+
+
+ /*
+ Quadrature global_qpts(nq, 3);
+ //cigma::array<double> qpts(nq,3);
+ shared_ptr<FE> fe(new FE(quadrature));
+
+ ApplyReferenceMap(fe, global_quadrature);
+
+ for (e = 0; e < nel; e++) // XXX: use cell_iterator
+ {
+ // load cell
+ mesh->getCell(e, *(fe->cell));
+ fe->update_jxw();
+
+ // apply refmap to local quadrature points
+ fe->apply_refmap(global_qpts);
+
+ // copy data to
+ for (q = 0; q < nq; q++)
+ {
+ const int n = e*nq + q;
+ weights->data(n) = quadrature->getWeight(q) * jxw[q];
+ for (j = 0; j < ndim; j++)
+ {
+ points->data(n,j) = qpts(q,j);
+ }
+ }
+ }
+
+ // */
+
+ return 0;
+}
+
Modified: cs/cigma/trunk/src/core_extract_op.h
===================================================================
--- cs/cigma/trunk/src/core_extract_op.h 2008-12-10 02:13:09 UTC (rev 13519)
+++ cs/cigma/trunk/src/core_extract_op.h 2008-12-10 02:13:10 UTC (rev 13520)
@@ -1,12 +1,13 @@
#ifndef __CIGMA_EXTRACT_OP_H__
#define __CIGMA_EXTRACT_OP_H__
+#include <boost/shared_ptr.hpp>
+
#include "core_base_op.h"
+#include "core_readers.h"
#include "MeshPart.h"
#include "Quadrature.h"
-//#include "Points.h"
-//#include "core_array.h"
-#include "nc_array.h"
+#include "core_array.h"
namespace cigma
{
@@ -16,16 +17,24 @@
class cigma::ExtractOp : public cigma::BaseOp
{
public:
+
ExtractOp();
~ExtractOp();
- void run();
+ void configure();
+ int run();
public:
+
+ MeshInfo mesh_info;
+ QuadratureInfo quadrature_info;
+
boost::shared_ptr<MeshPart> mesh;
boost::shared_ptr<Quadrature> quadrature;
- boost::shared_ptr<nc_array> points;
- //boost::shared_ptr<cigma::array> points;
+
+ cigma::array<double> weights;
+ cigma::array<double> points;
+
};
#endif
Modified: cs/cigma/trunk/src/core_list_op.cpp
===================================================================
--- cs/cigma/trunk/src/core_list_op.cpp 2008-12-10 02:13:09 UTC (rev 13519)
+++ cs/cigma/trunk/src/core_list_op.cpp 2008-12-10 02:13:10 UTC (rev 13520)
@@ -1,5 +1,7 @@
#include "core_list_op.h"
+#include "Common.h"
#include "Filesystem.h"
+//#include "io_hdf5_reader.h"
//#include "io_vtk_reader.h"
// std
@@ -127,6 +129,17 @@
}
+void ListOp::configure()
+{
+ fs::path p(filename);
+ if (!fs::exists(p))
+ {
+ cout << "File " << filename << " does not exist.";
+ throw cigma::Exception("ListOp::configure", "Could not list file");
+ }
+}
+
+
int ListOp::run()
{
cout << "ListOp::run()" << endl;
Modified: cs/cigma/trunk/src/core_list_op.h
===================================================================
--- cs/cigma/trunk/src/core_list_op.h 2008-12-10 02:13:09 UTC (rev 13519)
+++ cs/cigma/trunk/src/core_list_op.h 2008-12-10 02:13:10 UTC (rev 13520)
@@ -15,6 +15,7 @@
ListOp();
~ListOp();
+ void configure();
int run();
public:
Modified: cs/cigma/trunk/src/core_readers.cpp
===================================================================
--- cs/cigma/trunk/src/core_readers.cpp 2008-12-10 02:13:09 UTC (rev 13519)
+++ cs/cigma/trunk/src/core_readers.cpp 2008-12-10 02:13:10 UTC (rev 13520)
@@ -3,6 +3,8 @@
#include "nc_array.h"
#include "eb_array.h"
#include "AnnLocator.h"
+#include "Common.h"
+#include "io_file_reader.h"
#include <string>
using namespace std;
@@ -11,21 +13,6 @@
// ----------------------------------------------------------------------------
-DataPath MeshInfo::get_nc_path() const
-{
- // XXX: combine p_mesh w/ p_nc override
- return p_nc;
-}
-
-DataPath MeshInfo::get_eb_path() const
-{
- // XXX: combine p_mesh w/ p_eb override
- return p_eb;
-}
-
-// ----------------------------------------------------------------------------
-
-
static cigma::array<double>* ReadArray(const DataPath& array_path)
{
cigma::array<double>* a = 0;
@@ -46,7 +33,7 @@
else if (rt == FileReader::TEXT_FILE_READER)
{
a = new cigma::array<double>();
- int status = reader->getDataset(0, &(a->_data), &(a->npts), &(a->ndim));
+ int status = reader->getDataset("", &(a->_data), &(a->npts), &(a->ndim));
if (status < 0)
{
delete a;
@@ -64,9 +51,11 @@
// ----------------------------------------------------------------------------
-
shared_ptr<NodeCoordinates> cigma::ReadNodeCoordinates(const DataPath &nc_path)
{
+ TRI_LOG_STR("cigma::ReadNodeCoordinates()");
+ TRI_LOG(nc_path);
+
shared_ptr<NodeCoordinates> nc_ptr;
string filename = nc_path.filename();
string location = nc_path.location();
@@ -76,8 +65,13 @@
const FileReader::ReaderType rt = reader->getReaderType();
if (rt == FileReader::HDF5_FILE_READER)
{
- int status;
- throw cigma::Exception("ReadNodeCoordinates", "Error in HDF5 reader");
+ shared_ptr<nc_array> nc(new nc_array);
+ int status = reader->getCoordinates(location.c_str(), &(nc->coords), &(nc->nno), &(nc->nsd));
+ if (status < 0)
+ {
+ throw cigma::Exception("ReadNodeCoordinates", "Error in HDF5 reader");
+ }
+ nc_ptr = nc;
}
else if (rt == FileReader::VTK_FILE_READER)
{
@@ -87,16 +81,16 @@
else if (rt == FileReader::TEXT_FILE_READER)
{
shared_ptr<nc_array> nc(new nc_array);
- int status = reader->getCoordinates(0, &(nc->coords), &(nc->nno), &(nc->nsd));
+ int status = reader->getCoordinates("", &(nc->coords), &(nc->nno), &(nc->nsd));
if (status < 0)
{
- throw cigma::Exception("ReadNodeCoordinates", "error in text reader");
+ throw cigma::Exception("ReadNodeCoordinates", "Error in text reader");
}
nc_ptr = nc;
}
else
{
- throw cigma::Exception("ReadNodeCoordinates", "could not load data path");
+ throw cigma::Exception("ReadNodeCoordinates", "Could not load data path");
}
return nc_ptr;
@@ -105,9 +99,11 @@
// ----------------------------------------------------------------------------
-
shared_ptr<ElementBlock> cigma::ReadElementBlock(const DataPath& eb_path)
{
+ TRI_LOG_STR("cigma::ReadElementBlock()");
+ TRI_LOG(eb_path);
+
shared_ptr<ElementBlock> eb_ptr;
string filename = eb_path.filename();
string location = eb_path.location();
@@ -117,8 +113,13 @@
const FileReader::ReaderType rt = reader->getReaderType();
if (rt == FileReader::HDF5_FILE_READER)
{
- int status;
- throw cigma::Exception("ReadElementBlock", "Error in HDF5 reader");
+ shared_ptr<eb_array> eb(new eb_array);
+ int status = reader->getConnectivity(location.c_str(), &(eb->connect), &(eb->nel), &(eb->ndofs));
+ if (status < 0)
+ {
+ throw cigma::Exception("ReadElementBlock", "Error in HDF5 reader");
+ }
+ eb_ptr = eb;
}
else if (rt == FileReader::VTK_FILE_READER)
{
@@ -128,27 +129,27 @@
else if (rt == FileReader::TEXT_FILE_READER)
{
shared_ptr<eb_array> eb(new eb_array);
- int status = reader->getConnectivity(0, &(eb->connect), &(eb->nel), &(eb->ndofs));
+ int status = reader->getConnectivity("", &(eb->connect), &(eb->nel), &(eb->ndofs));
if (status < 0)
{
- throw cigma::Exception("ReadElementBlock", "error in text reader");
+ throw cigma::Exception("ReadElementBlock", "Error in text reader");
}
eb_ptr = eb;
}
else
{
- throw cigma::Exception("ReadElementBlock", "could not load data path");
+ throw cigma::Exception("ReadElementBlock", "Could not load data path");
}
return eb_ptr;
}
-
// ----------------------------------------------------------------------------
-
shared_ptr<MeshPart> cigma::ReadMeshPart(const MeshInfo& mesh_info)
{
+ TRI_LOG_STR("cigma::ReadMeshPart()");
+
shared_ptr<MeshPart> mesh_ptr;
shared_ptr<NodeCoordinates> nc_ptr;
shared_ptr<ElementBlock> eb_ptr;
@@ -157,11 +158,27 @@
// load coordinates
DataPath nc_path = mesh_info.get_nc_path();
- nc_ptr = ReadNodeCoordinates(nc_path);
+ try
+ {
+ nc_ptr = ReadNodeCoordinates(nc_path);
+ }
+ catch (cigma::Exception& e)
+ {
+ string msg("Could not read NodeCoordinates");
+ throw cigma::Exception("ReadMeshPart", msg);
+ }
// load elements
DataPath eb_path = mesh_info.get_eb_path();
- eb_ptr = ReadElementBlock(eb_path);
+ try
+ {
+ eb_ptr = ReadElementBlock(eb_path);
+ }
+ catch (cigma::Exception& e)
+ {
+ string msg("Could not read ElementBlock");
+ throw cigma::Exception("ReadMeshPart", msg);
+ }
// load cell type from name
cell_type = Cell::string2type(mesh_info.cell_type_name);
@@ -182,13 +199,12 @@
return mesh_ptr;
}
-
-
// ----------------------------------------------------------------------------
-
shared_ptr<Quadrature> cigma::ReadQuadrature(const QuadratureInfo& quadrature_info)
{
+ TRI_LOG_STR("cigma::ReadQuadrature()");
+
shared_ptr<Quadrature> q_ptr(new Quadrature);
Cell::type cell_type = Cell::NONE;
@@ -196,47 +212,56 @@
{
cell_type = Cell::string2type(quadrature_info.cell_type_name);
}
+ q_ptr->setCellType(cell_type);
if (quadrature_info.p_quadrature.empty() &&
quadrature_info.p_weights.empty() &&
quadrature_info.p_points.empty())
{
+ //
+ // No paths were specified, so we should create a
+ // new Quadrature object based on the given cell type
+ //
if (cell_type != Cell::NONE)
{
q_ptr = Quadrature::default_rule(cell_type);
}
else
{
- throw cigma::Exception("ReadQuadrature", "");
+ throw cigma::Exception("ReadQuadrature", "No cell type specified!");
}
}
else if (quadrature_info.p_quadrature.empty())
{
+ //
+ // Since the p_quadrature path is empty, we use the other
+ // two paths p_weights and p_points to create the new Quadrature
+ //
if (quadrature_info.p_weights.empty())
{
- throw cigma::Exception("ReadQuadrature", "");
+ throw cigma::Exception("ReadQuadrature", "No path to quadrature weights given.");
}
if (quadrature_info.p_points.empty())
{
- throw cigma::Exception("ReadQuadrature", "");
+ throw cigma::Exception("ReadQuadrature", "No path to quadrature points given.");
}
- // load weights
+ // load weights array from path
cigma::array<double> *wts = ReadArray(quadrature_info.p_weights);
if (wts == 0)
{
- throw cigma::Exception("ReadQuadrature", "");
+ throw cigma::Exception("ReadQuadrature", "failed to read quadrature weights!");
}
q_ptr->weights = wts->_data;
wts->_data = 0;
delete wts;
- // load points
+ // load points array from path
cigma::array<double> *pts = ReadArray(quadrature_info.p_points);
if (pts == 0)
{
- throw cigma::Exception("ReadQuadrature", "");
+ throw cigma::Exception("ReadQuadrature", "failed to read quadrature points!");
}
q_ptr->points = pts->_data;
pts->_data = 0;
@@ -244,14 +269,17 @@
}
else
{
+ //
+ // Now, we know p_quadrature is not empty so we use it exclusively
+ //
if (!quadrature_info.p_weights.empty())
{
- throw cigma::Exception("ReadQuadrature", "");
+ throw cigma::Exception("ReadQuadrature", "quadrature path already specified");
}
if (!quadrature_info.p_points.empty())
{
- throw cigma::Exception("ReadQuadrature", "");
+ throw cigma::Exception("ReadQuadrature", "quadrature path already specified");
}
// load weights & points simultaneously
@@ -273,12 +301,11 @@
return q_ptr;
}
-
// ----------------------------------------------------------------------------
-
shared_ptr<FE> cigma::ReadFE(const FE_Info& fe_info)
{
+ TRI_LOG_STR("cigma::ReadFE()");
assert(fe_info.q_info.cell_type_name != "");
shared_ptr<FE> fe_ptr;
@@ -303,12 +330,13 @@
return fe_ptr;
}
-
// ----------------------------------------------------------------------------
-
shared_ptr<DofHandler> cigma::ReadDofHandler(const DataPath& dofs_path)
{
+ TRI_LOG_STR("cigma::ReadDofHandler()");
+ TRI_LOG(dofs_path);
+
shared_ptr<DofHandler> dofs_ptr;
string filename = dofs_path.filename();
string location = dofs_path.location();
@@ -318,8 +346,12 @@
const FileReader::ReaderType rt = reader->getReaderType();
if (rt == FileReader::HDF5_FILE_READER)
{
- int status;
- throw cigma::Exception("ReadDofHandler", "Error in HDF5 reader");
+ dofs_ptr = shared_ptr<DofHandler>(new DofHandler);
+ int status = reader->getDataset(location.c_str(), &(dofs_ptr->dofs), &(dofs_ptr->nno), &(dofs_ptr->rank));
+ if (status < 0)
+ {
+ throw cigma::Exception("ReadDofHandler", "Error in HDF5 reader");
+ }
}
else if (rt == FileReader::VTK_FILE_READER)
{
@@ -329,7 +361,7 @@
else if (rt == FileReader::TEXT_FILE_READER)
{
dofs_ptr = shared_ptr<DofHandler>(new DofHandler);
- int status = reader->getDataset(0, &(dofs_ptr->dofs), &(dofs_ptr->nno), &(dofs_ptr->rank));
+ int status = reader->getDataset("", &(dofs_ptr->dofs), &(dofs_ptr->nno), &(dofs_ptr->rank));
if (status < 0)
{
throw cigma::Exception("ReadDofHandler", "Error in text reader");
@@ -343,24 +375,24 @@
return dofs_ptr;
}
-
-
// ----------------------------------------------------------------------------
shared_ptr<Field> cigma::ReadField(const FieldInfo& field_info)
{
+ TRI_LOG_STR("cigma::ReadField()");
+
shared_ptr<Field> field_ptr;
shared_ptr<MeshPart> mesh_ptr;
shared_ptr<FE> fe_ptr;
shared_ptr<DofHandler> dofs_ptr;
- // load mesh
+ // load MeshPart from mesh_info
mesh_ptr = ReadMeshPart(field_info.mesh_info);
- // load FE
+ // load FE from fe_info
fe_ptr = ReadFE(field_info.fe_info);
- // load DofHandler
+ // load DofHandler from path p_field
dofs_ptr = ReadDofHandler(field_info.p_field);
// initialize Field
Modified: cs/cigma/trunk/src/core_readers.h
===================================================================
--- cs/cigma/trunk/src/core_readers.h 2008-12-10 02:13:09 UTC (rev 13519)
+++ cs/cigma/trunk/src/core_readers.h 2008-12-10 02:13:10 UTC (rev 13520)
@@ -11,120 +11,11 @@
#include "Field.h"
#include "Cell.h"
#include "DataPath.h"
-#include "io_file_reader.h"
-#include "Exception.h"
+#include "core_args.h"
-
namespace cigma
{
- struct MeshInfo
- {
- DataPath p_mesh;
- DataPath p_nc;
- DataPath p_eb;
- bool use_locator;
- std::string cell_type_name;
-
- MeshInfo() : use_locator(false) {}
- ~MeshInfo() {}
- MeshInfo(const MeshInfo& other)
- {
- p_mesh = other.p_mesh;
- p_nc = other.p_nc;
- p_eb = other.p_eb;
- use_locator = other.use_locator;
- cell_type_name = other.cell_type_name;
- }
- MeshInfo& operator=(const MeshInfo& other)
- {
- p_mesh = other.p_mesh;
- p_nc = other.p_nc;
- p_eb = other.p_eb;
- use_locator = other.use_locator;
- cell_type_name = other.cell_type_name;
- }
- //bool single_file() const;
- DataPath get_nc_path() const;
- DataPath get_eb_path() const;
- //FileReader::ReaderType get_nc_reader_type() const;
- //FileReader::ReaderType get_eb_reader_type() const;
- };
-
- struct QuadratureInfo
- {
- DataPath p_quadrature;
- DataPath p_weights;
- DataPath p_points;
- std::string cell_type_name;
-
- QuadratureInfo() {}
- ~QuadratureInfo() {}
- QuadratureInfo(const QuadratureInfo& other)
- {
- p_quadrature = other.p_quadrature;
- p_weights = other.p_weights;
- p_points = other.p_points;
- cell_type_name = other.cell_type_name;
- }
- QuadratureInfo& operator=(const QuadratureInfo& other)
- {
- p_quadrature = other.p_quadrature;
- p_weights = other.p_weights;
- p_points = other.p_points;
- cell_type_name = other.cell_type_name;
- }
- //bool single_file() const;
- //DataPath get_weights_path() const;
- //DataPath get_points_path() const;
- //FileReader::ReaderType get_weights_reader_type() const;
- //FileReader::ReaderType get_points_reader_type() const;
- };
-
- struct FE_Info
- {
- QuadratureInfo q_info;
- DataPath p_fe_basis;
- DataPath p_fe_basis_jet;
-
- FE_Info() {}
- ~FE_Info() {}
- FE_Info(const FE_Info& other)
- {
- q_info = other.q_info;
- p_fe_basis = other.p_fe_basis;
- p_fe_basis_jet = other.p_fe_basis_jet;
- }
- FE_Info& operator=(const FE_Info& other)
- {
- q_info = other.q_info;
- p_fe_basis = other.p_fe_basis;
- p_fe_basis_jet = other.p_fe_basis_jet;
- }
- };
-
- struct FieldInfo
- {
- DataPath p_field;
- FE_Info fe_info;
- MeshInfo mesh_info;
-
- FieldInfo() {}
- ~FieldInfo() {}
- FieldInfo(const FieldInfo& other)
- {
- p_field = other.p_field;
- fe_info = other.fe_info;
- mesh_info = other.mesh_info;
- }
- FieldInfo& operator=(const FieldInfo& other)
- {
- p_field = other.p_field;
- fe_info = other.fe_info;
- mesh_info = other.mesh_info;
- }
- };
-
boost::shared_ptr<NodeCoordinates> ReadNodeCoordinates(const DataPath& nc_path);
boost::shared_ptr<ElementBlock> ReadElementBlock(const DataPath& eb_path);
@@ -138,6 +29,7 @@
boost::shared_ptr<DofHandler> ReadDofHandler(const DataPath& dofs_path);
boost::shared_ptr<Field> ReadField(const FieldInfo& field_info);
+
}
#endif
Modified: cs/cigma/trunk/src/core_writers.cpp
===================================================================
--- cs/cigma/trunk/src/core_writers.cpp 2008-12-10 02:13:09 UTC (rev 13519)
+++ cs/cigma/trunk/src/core_writers.cpp 2008-12-10 02:13:10 UTC (rev 13520)
@@ -1,2 +1,196 @@
#include "core_writers.h"
+#include "io_text_writer.h"
+#include "io_hdf5_writer.h"
+#include "tri_logger.hpp"
+using namespace cigma;
+#include <boost/shared_ptr.hpp>
+using boost::shared_ptr;
+
+using namespace std;
+
+
+void cigma::WriteArray(const DataPath& path, const cigma::array<double>& x)
+{
+ TRI_LOG_STR("cigma::WriteArray()");
+ TRI_LOG(path);
+
+ std::string filename = path.filename();
+ shared_ptr<FileWriter> writer = FileWriter::New(filename, "w");
+ const FileWriter::WriterType wt = writer->getWriterType();
+ if (wt == FileWriter::HDF5_FILE_WRITER)
+ {
+ throw cigma::Exception("WriteArray", "need HDF5 writer");
+ }
+ else if (wt == FileWriter::VTK_FILE_WRITER)
+ {
+ throw cigma::Exception("WriteArray", "need VTK writer");
+ }
+ else if (wt == FileWriter::TEXT_FILE_WRITER)
+ {
+ const double *data = x._data;
+ int status = writer->writeDataset("", data, x.n_points(), x.n_dim());
+ if (status < 0)
+ {
+ throw cigma::Exception("WriteArray", "could not write dataset to text file");
+ }
+ }
+ else
+ {
+ throw cigma::Exception("WriteArray", "could not write array to path");
+ }
+}
+
+void cigma::WriteNodeCoordinates(const DataPath& path, shared_ptr<NodeCoordinates> nc)
+{
+ TRI_LOG_STR("cigma::WriteNodeCoordinates()");
+ TRI_LOG(path);
+ throw cigma::Exception("WriteNodeCoordinates", "need implementation");
+}
+
+void cigma::WriteElementBlock(const DataPath& path, shared_ptr<ElementBlock> eb)
+{
+ TRI_LOG_STR("cigma::WriteElementBlock()");
+ TRI_LOG(path);
+ throw cigma::Exception("WriteElementBlock", "need implementation");
+}
+
+void cigma::WriteQuadrature(const DataPath& path, shared_ptr<Quadrature> Q)
+{
+ //throw cigma::Exception("WriteQuadrature", "need implementation");
+
+ TRI_LOG_STR("cigma::WriteQuadrature()");
+ TRI_LOG(path);
+
+ std::string filename = path.filename();
+ shared_ptr<FileWriter> writer = FileWriter::New(filename, "w");
+ const FileWriter::WriterType wt = writer->getWriterType();
+ if (wt == FileWriter::HDF5_FILE_WRITER)
+ {
+ int status;
+ string pts_loc = path.location() + "/points";
+ string wts_loc = path.location() + "/weights";
+ HDF5_Writer *hdf5_writer = static_cast<HDF5_Writer*>(&(*writer));
+
+ status = hdf5_writer->writeDataset(wts_loc.c_str(), Q->weights, Q->npts, 1);
+ if (status < 0)
+ {
+ throw cigma::Exception("WriteQuadrature", "could not write weights to HDF5 file");
+ }
+
+ status = hdf5_writer->writeDataset(pts_loc.c_str(), Q->points, Q->npts, Q->ndim);
+ if (status < 0)
+ {
+ throw cigma::Exception("WriteQuadrature", "could not write points to HDF5 file");
+ }
+ }
+ else if (wt == FileWriter::VTK_FILE_WRITER)
+ {
+ throw cigma::Exception("WriteQuadrature", "need VTK writer");
+ }
+ else if (wt == FileWriter::TEXT_FILE_WRITER)
+ {
+ TextWriter *text_writer = static_cast<TextWriter*>(&(*writer));
+ FILE *fp = text_writer->fp;
+ if (fp == NULL)
+ {
+ throw cigma::Exception("WriteQuadrature", "could not write dataset to text file");
+ }
+ const int nq = Q->n_points();
+ const int dim = Q->n_dim();
+ fprintf(fp, "%d %d\n", nq, 1+dim);
+ for (int q = 0; q < nq; q++)
+ {
+ fprintf(fp, " %g", Q->getWeight(q));
+ for (int j = 0; j < dim; j++)
+ fprintf(fp, " %g", Q->getPoint(q,j));
+ fprintf(fp, "\n");
+ }
+ }
+ else
+ {
+ throw cigma::Exception("WriteQuadrature", "could not write quadrature to path");
+ }
+}
+
+void cigma::WriteQPoints(const DataPath& path, const cigma::array<double>& w, const cigma::array<double>& x)
+{
+ shared_ptr<Quadrature> Q(new Quadrature);
+
+ assert(w.n_points() > 0);
+ assert(x.n_points() > 0);
+ assert(w.n_points() == x.n_points());
+ assert(w.n_dim() == 1);
+
+ Q->npts = x.n_points();
+ Q->ndim = x.n_dim();
+ Q->points = x._data;
+ Q->weights = w._data;
+
+ try
+ {
+ WriteQuadrature(path, Q);
+ }
+ catch (cigma::Exception& e)
+ {
+ Q->points = 0;
+ Q->weights = 0;
+ throw e;
+ }
+
+ Q->points = 0;
+ Q->weights = 0;
+}
+
+void cigma::WriteDofs(const DataPath& path, shared_ptr<DofHandler> dofs)
+{
+ TRI_LOG_STR("cigma::WriteDofs()");
+ TRI_LOG(path);
+ throw cigma::Exception("WriteDofs", "need implementation");
+}
+
+void cigma::WriteResiduals(const DataPath& path, shared_ptr<Residuals> residuals)
+{
+ TRI_LOG_STR("cigma::WriteResiduals()");
+ TRI_LOG(path);
+
+ std::string filename = path.filename();
+ 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
+
+ TRI_LOG(eps_loc);
+
+ const double* epsilon = residuals->epsilon;
+ status = hdf5_writer->writeDataset(eps_loc.c_str(), epsilon, residuals->nel, 1);
+ if (status < 0)
+ {
+ throw cigma::Exception("WriteResiduals", "could not write residuals to HDF5 file");
+ }
+ }
+ else if (wt == FileWriter::VTK_FILE_WRITER)
+ {
+ throw cigma::Exception("WriteResiduals", "need VTK writer");
+ }
+ else if (wt == FileWriter::TEXT_FILE_WRITER)
+ {
+ const double *epsilon = residuals->epsilon;
+ int status = writer->writeDataset("", epsilon, residuals->nel, 1);
+ if (status < 0)
+ {
+ throw cigma::Exception("WriteResiduals", "could not write residuals to text file");
+ }
+ }
+ else
+ {
+ throw cigma::Exception("WriteResiduals", "could not write residuals to path");
+ }
+}
+
Modified: cs/cigma/trunk/src/core_writers.h
===================================================================
--- cs/cigma/trunk/src/core_writers.h 2008-12-10 02:13:09 UTC (rev 13519)
+++ cs/cigma/trunk/src/core_writers.h 2008-12-10 02:13:10 UTC (rev 13520)
@@ -2,10 +2,31 @@
#define __CIGMA_WRITERS_H__
#include <boost/shared_ptr.hpp>
+#include "DataPath.h"
+#include "core_array.h"
+#include "NodeCoordinates.h"
+#include "ElementBlock.h"
+#include "Quadrature.h"
+#include "DofHandler.h"
+//#include "MeshPart.h"
+#include "Residuals.h"
+#include "io_file_writer.h"
namespace cigma
{
- void WriteResiduals(boost::shared<Residuals> residuals, const DataPath& path);
+ void WriteArray(const DataPath& path, const cigma::array<double>& x);
+
+ void WriteNodeCoordinates(const DataPath& path, boost::shared_ptr<NodeCoordinates> nc);
+
+ void WriteElementBlock(const DataPath& path, boost::shared_ptr<ElementBlock> eb);
+
+ void WriteQuadrature(const DataPath& path, boost::shared_ptr<Quadrature> Q);
+
+ void WriteQPoints(const DataPath& path, const cigma::array<double>& w, const cigma::array<double>& x);
+
+ void WriteDofs(const DataPath& path, boost::shared_ptr<DofHandler> dofs);
+
+ void WriteResiduals(const DataPath& path, boost::shared_ptr<Residuals> residuals);
}
#endif
More information about the CIG-COMMITS
mailing list