[cig-commits] r11480 - in cs/benchmark/cigma/trunk/src: . hold
luis at geodynamics.org
luis at geodynamics.org
Wed Mar 19 12:00:16 PDT 2008
Author: luis
Date: 2008-03-19 12:00:15 -0700 (Wed, 19 Mar 2008)
New Revision: 11480
Added:
cs/benchmark/cigma/trunk/src/hold/
cs/benchmark/cigma/trunk/src/hold/CompareCmd.cpp
cs/benchmark/cigma/trunk/src/hold/CompareCmd.h
cs/benchmark/cigma/trunk/src/hold/EvalCmd.cpp
cs/benchmark/cigma/trunk/src/hold/EvalCmd.h
cs/benchmark/cigma/trunk/src/hold/ExtractCmd.cpp
cs/benchmark/cigma/trunk/src/hold/ExtractCmd.h
cs/benchmark/cigma/trunk/src/hold/FieldIO.cpp
cs/benchmark/cigma/trunk/src/hold/FieldIO.h
cs/benchmark/cigma/trunk/src/hold/HdfArray.cpp
cs/benchmark/cigma/trunk/src/hold/HdfArray.h
cs/benchmark/cigma/trunk/src/hold/HdfGroup.cpp
cs/benchmark/cigma/trunk/src/hold/HdfGroup.h
cs/benchmark/cigma/trunk/src/hold/HdfReader.cpp
cs/benchmark/cigma/trunk/src/hold/HdfReader.h
cs/benchmark/cigma/trunk/src/hold/HdfWriter.cpp
cs/benchmark/cigma/trunk/src/hold/HdfWriter.h
cs/benchmark/cigma/trunk/src/hold/Makefile.am
cs/benchmark/cigma/trunk/src/hold/MeshIO.cpp
cs/benchmark/cigma/trunk/src/hold/MeshIO.h
cs/benchmark/cigma/trunk/src/hold/QuadratureIO.cpp
cs/benchmark/cigma/trunk/src/hold/QuadratureIO.h
cs/benchmark/cigma/trunk/src/hold/Reader.cpp
cs/benchmark/cigma/trunk/src/hold/Reader.h
cs/benchmark/cigma/trunk/src/hold/TextReader.cpp
cs/benchmark/cigma/trunk/src/hold/TextReader.h
cs/benchmark/cigma/trunk/src/hold/TextWriter.cpp
cs/benchmark/cigma/trunk/src/hold/TextWriter.h
cs/benchmark/cigma/trunk/src/hold/VtkList.cpp
cs/benchmark/cigma/trunk/src/hold/VtkList.h
cs/benchmark/cigma/trunk/src/hold/VtkReader.cpp
cs/benchmark/cigma/trunk/src/hold/VtkReader.h
cs/benchmark/cigma/trunk/src/hold/VtkWriter.cpp
cs/benchmark/cigma/trunk/src/hold/VtkWriter.h
cs/benchmark/cigma/trunk/src/hold/VtkXmlReader.cpp
cs/benchmark/cigma/trunk/src/hold/VtkXmlReader.h
cs/benchmark/cigma/trunk/src/hold/Writer.cpp
cs/benchmark/cigma/trunk/src/hold/Writer.h
Log:
Temporarily place all the I/O classes on hold
Added: cs/benchmark/cigma/trunk/src/hold/CompareCmd.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/CompareCmd.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/CompareCmd.cpp 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,462 @@
+#include <iostream>
+#include <cassert>
+#include "CompareCmd.h"
+#include "StringUtils.h"
+
+#include "MeshPart.h"
+#include "Tet.h"
+#include "Hex.h"
+#include "Numeric.h"
+#include "AnnLocator.h"
+
+#include "HdfReader.h"
+#include "VtkReader.h"
+#include "VtkWriter.h"
+
+#include "Timer.h"
+#include "Misc.h"
+
+
+using namespace std;
+using namespace cigma;
+
+
+
+// ---------------------------------------------------------------------------
+
+cigma::CompareCmd::CompareCmd()
+{
+ name = "compare";
+
+ // integrating mesh
+ mesh = 0;
+ quadrature = 0;
+ qr = 0;
+
+ // fields
+ field_a = 0;
+ field_b = 0;
+
+ // residuals
+ residuals = new ResidualField();
+
+ // parameters
+ verbose = false;
+ output_frequency = 0;
+}
+
+cigma::CompareCmd::~CompareCmd()
+{
+ /* XXX: for now, don't delete
+ if (qr != 0)
+ {
+ delete qr; // XXX
+ } // */
+
+ if (residuals != 0) delete residuals;
+}
+
+
+// ---------------------------------------------------------------------------
+
+void cigma::CompareCmd::setupOptions(AnyOption *opt)
+{
+ //std::cout << "Calling cigma::CompareCmd::setupOptions()" << std::endl;
+
+ assert(opt != 0);
+
+ /* setup usage */
+ opt->addUsage("Usage:");
+ opt->addUsage(" cigma compare [options]");
+ opt->addUsage(" -a --first First field location");
+ opt->addUsage(" -b --second Second field location");
+ opt->addUsage(" -r --rule Quadrature rule location");
+ opt->addUsage(" -o --output Output file");
+
+ /* setup flags and options */
+ opt->setFlag("help", 'h');
+ opt->setFlag("debug"); // XXX: uses specific defaults
+
+ // options for mesh
+ opt->setOption("mesh");
+ opt->setOption("mesh-coordinates");
+ opt->setOption("mesh-connectivity");
+
+ // options for quadrature
+ opt->setOption("rule",'r');
+ opt->setOption("rule-order");
+ opt->setOption("rule-points");
+ opt->setOption("rule-weights");
+
+ // options for first field, in the form /path/to/file:/path/to/dset
+ opt->setOption("first",'a');
+ opt->setOption("first-mesh");
+ opt->setOption("first-mesh-coordinates");
+ opt->setOption("first-mesh-connectivity");
+
+ // option for second field, in the form /path/to/file:/path/to/dset
+ opt->setOption("second",'b');
+ opt->setOption("second-mesh");
+ opt->setOption("second-mesh-coordinates");
+ opt->setOption("second-mesh-connectivity");
+
+ // options for output
+ opt->setOption("output",'o');
+
+ // other options
+ opt->setFlag("verbose",'v');
+ opt->setOption("output-frequency",'f');
+}
+
+
+void cigma::CompareCmd::configure(AnyOption *opt)
+{
+ //std::cout << "Calling cigma::CompareCmd::configure()" << std::endl;
+
+ string field_prefix;
+ std::string inputstr;
+ char *in;
+
+
+ /* Check if --verbose was enabled */
+
+ verbose = opt->getFlag("verbose");
+ if (verbose)
+ {
+ output_frequency = 1000;
+ }
+
+ /* Determine the --output-frequency option */
+
+ in = opt->getValue("output-frequency");
+ if (in != 0)
+ {
+ inputstr = in;
+ string_to_int(inputstr, output_frequency);
+ }
+
+ if (output_frequency == 0)
+ {
+ // XXX: emit warning, or quit?
+ if (opt->getValue("output-frequency") != 0)
+ {
+ cerr << "compare: Warning: ignoring option --output-frequency" << endl;
+ }
+ verbose = false;
+ }
+
+ /* Determine the --output option */
+
+ in = opt->getValue("output");
+ if (in == 0)
+ {
+ if (opt->getFlag("debug"))
+ {
+ // XXX: provide default name when in debug mode
+ in = (char *)"foo.vtk";
+ }
+ else
+ {
+ cerr << "compare: Please specify the option --output" << endl;
+ exit(1);
+ }
+ }
+ output_path = in;
+
+
+ /*
+ * Initialization order:
+ * Load First field
+ * If FE_Field
+ * Load DofsA
+ * Load MeshA if req'd
+ * Load RuleA if req'd
+ * Else
+ * Load Analytic Field
+ * Load Second field
+ * If FE_Field
+ * Load DofsB
+ * Load MeshB if req'd
+ * Load RuleB if req'd
+ * Else
+ * Load Analytic Field
+ * Load Integration mesh
+ * Load Quadrature rule
+ */
+
+
+ /* Gather up the expected command line arguments */
+
+ load_args(opt, &meshIO, "mesh");
+ load_args(opt, &quadratureIO, "rule");
+ load_args(opt, &firstIO, "first");
+ load_args(opt, &secondIO, "second");
+
+
+ /* Validate these arguments and complain about missing ones */
+
+ if (opt->getFlag("debug"))
+ {
+ //
+ // assign defaults if we're in debug mode. this overrides
+ // the command line settings from load_args()
+ //
+ if (firstIO.field_path == "")
+ firstIO.field_path = "./tests/strikeslip_tet4_1000m_t0.vtk:displacements_t0";
+
+ if (secondIO.field_path == "")
+ secondIO.field_path = "./tests/strikeslip_hex8_1000m_t0.vtk:displacements_t0";
+ }
+
+ validate_args(&firstIO, "compare");
+ validate_args(&secondIO, "compare");
+ validate_args(&meshIO, "compare");
+ validate_args(&quadratureIO, "compare");
+
+
+ /* Load the datasets into memory */
+
+ firstIO.load();
+ field_a = firstIO.field;
+ assert(field_a != 0);
+ //field_a->fe = new FE();
+ //field_a->fe->set_cell(field_a->meshPart->cell);
+ cout << "first field path = " << firstIO.field_path << endl;
+ cout << "first field dimensions = "
+ << field_a->meshPart->nel << " cells, "
+ << field_a->meshPart->nno << " nodes, "
+ << field_a->meshPart->cell->n_nodes() << " dofs/cell, "
+ << "rank " << field_a->n_rank() << endl;
+
+ secondIO.load();
+ field_b = secondIO.field;
+ assert(field_b != 0);
+ //field_b->fe = new FE();
+ //field_b->fe->set_cell(field_b->meshPart->cell);
+ cout << "second field path = " << secondIO.field_path << endl;
+ cout << "second field dimensions = "
+ << field_b->meshPart->nel << " cells, "
+ << field_b->meshPart->nno << " nodes, "
+ << field_b->meshPart->cell->n_nodes() << " dofs/cell, "
+ << "rank " << field_b->n_rank() << endl;
+
+
+ /* XXX: if no --mesh option was specified, get mesh from the
+ * first field. if the first field has no mesh, (e.g. specified
+ * by an analytic soln), then try using the second field's mesh.
+ * if we still don't have a mesh that we can use for the integration,
+ * then exit with an error suggesting the user to specify the --mesh
+ * option.
+ */
+ meshIO.load();
+ mesh = meshIO.meshPart;
+ if (mesh == 0)
+ {
+ mesh = firstIO.field->meshPart;
+ }
+ assert(mesh != 0);
+ residuals->set_mesh(mesh);
+
+
+ /* Now load the quadrature rule. If no rule is specified on the
+ * command line, a default rule is assigned based on the type of
+ * the cell. Also, an exception is thrown if the specified a rule
+ * does not conform geometrically to the interior of the cell.
+ */
+ quadratureIO.load(mesh->cell);
+ quadrature = quadratureIO.quadrature;
+ assert(quadrature != 0);
+ //field_a->fe->set_quadrature(quadrature); // XXX: eliminate need for this statement
+
+
+ field_a->fe = new FE();
+ field_a->fe->set_mesh(field_a->meshPart);
+ //field_a->fe->set_quadrature_points(field_a->meshPart->cell->default_quadrature_points);
+ field_a->fe->set_quadrature_points(quadrature); // XXX
+
+
+ field_b->fe = new FE();
+ field_b->fe->set_mesh(field_b->meshPart);
+ //field_b->fe->set_quadrature_points(field_b->meshPart->cell->default_quadrature_points);
+
+
+ // XXX: if field_a has a quadrature rule, reuse it
+ // XXX: however, this complicates our destructor...
+ qr = field_a->fe;
+ assert(field_a->meshPart == mesh); // XXX
+
+
+ /* XXX: when neither field_a or field_b have a mesh, need own QuadratureRule instance
+ assert(mesh != 0);
+ assert(field_a->meshPart == 0);
+ assert(field_b->meshPart == 0);
+ qr = new QuadratureRule();
+ qr->set_mesh(mesh);
+ qr->set_quadrature_points(quadrature);
+ // */
+
+
+
+ return;
+}
+
+
+// ---------------------------------------------------------------------------
+
+void compare(CompareCmd *env, MeshPart *mesh, FE_Field *field_a, FE_Field *field_b)
+{
+ assert(env != 0);
+ const int output_frequency = env->output_frequency;
+ const bool verbose = env->verbose;
+
+ // XXX: remove need for this assert stmt ...
+ assert(mesh == field_a->meshPart);
+
+ // XXX: remove need for env...consider moving this function back to run()
+ ResidualField *residuals = env->residuals;
+
+ // XXX: we need need to initialize qr much earlier, so we avoid
+ // assuming that we can load the quadrature points from field_a->fe object
+ //QuadraturePoints *quadrature = field_a->fe->quadrature;
+ //QuadratureRule *qr = new QuadratureRule();
+ //qr->meshPart = mesh;
+ //qr->set_quadrature_points(quadrature);
+ QuadratureRule *qr = env->qr; // XXX: replace MeshPart fn arg by QuadratureRule
+ assert(qr != 0);
+ assert(qr->points != 0);
+
+ // dimensions
+ int nel = qr->meshPart->nel;
+ int nq = qr->points->n_points();
+ int valdim = field_a->n_rank();
+
+ // field values at quadrature points
+ Points phi_a, phi_b;
+ double *values_a = new double[nq * valdim];
+ double *values_b = new double[nq * valdim];
+ phi_a.set_data(values_a, nq, valdim);
+ phi_b.set_data(values_b, nq, valdim);
+
+ // reset error accumulator
+ residuals->zero_out();
+
+ // start timer
+ Timer timer;
+ if (verbose)
+ {
+ cout << std::setprecision(5);
+ timer.print_header(cout, "elts");
+ timer.start(nel);
+ timer.update(0);
+ cout << timer;
+ }
+
+ for (int e = 0; e < nel; e++)
+ {
+ // update cell data
+ qr->select_cell(e);
+
+ // ... calculate phi_a
+ //field_a->Field::eval(*(qr->points), phi_a);
+ field_a->tabulate_element(e, phi_a.data); // XXX: call when field_b is an instance of FE_Field
+
+ // ... calculate phi_b
+ field_b->Field::eval(*(qr->points), phi_b);
+ //field_b->tabulate_element(e, phi_b.data); // XXX: call when field_b->meshPart, field_a->meshPart,
+ // and qr->meshPart are all identical...
+ // this will result in redundant call to get_cell_dofs().
+ // is it worth it to optimize this? it would avoid
+ // jumping around the connectivity array a second time..
+
+ // apply quadrature rule
+ double err = qr->L2(phi_a, phi_b);
+
+ // update error accumulator
+ residuals->update(e, err);
+
+ // debug info
+ if (verbose && ((e+1) % output_frequency == 0))
+ {
+ timer.update(e+1);
+ cout << timer;
+ }
+ }
+
+ if (verbose)
+ {
+ timer.update(nel);
+ cout << timer << endl;
+ }
+
+ // report global error
+ double L2 = residuals->L2();
+ cout << setprecision(12);
+ cout << "L2 = " << L2 << endl;
+
+ // write out data
+ residuals->write_vtk(env->output_path.c_str());
+
+ // clean up
+ delete [] values_a;
+ delete [] values_b;
+}
+
+/* XXX: first, time the effects of introducing branches in the main loop
+void compare(MeshPart *mesh, FE_Field *field_a, Field *field_b,
+ ResidualField *residuals, int output_frequency)
+{
+}
+
+void compare(FE_Field *field_a, Field *field_b,
+ ResidualField *residuals, int output_frequency)
+{
+}
+
+void compare(FE_Field *field_a, PointField *field_b,
+ ResidualField *residuals, int output_frequency)
+{
+}
+
+void compare(MeshPart *mesh, Field *field_a, PointField *field_b,
+ ResidualField *residuals, int output_frequency)
+{
+}
+// */
+
+
+//* XXX: new run() method
+int cigma::CompareCmd::run()
+{
+ //std::cout << "Calling cigma::CompareCmd::run()" << std::endl;
+
+ //
+ // XXX: need to fail gracefully at this point, instead of throwing
+ // a nasty assert statement whenever something unexpected happens.
+ // clean up memory management so that we can use exceptions, and
+ // verify that all resources are properly finalized.
+ //
+
+
+ // start with the obvious checks
+ assert(quadrature != 0);
+ assert(mesh != 0);
+ assert(qr != 0);
+ assert(field_a != 0);
+ assert(field_b != 0);
+
+
+ // make sure the field dimensions match
+ assert(field_a->n_dim() == field_a->n_dim());
+ assert(field_a->n_rank() == field_b->n_rank());
+
+
+ // XXX: call the appropriate compare method based on the field types
+ // XXX: enclose statements in a try-catch-finally clause
+ compare(this, mesh, field_a, field_b);
+
+
+ return 0;
+} // */
+
+
+// ---------------------------------------------------------------------------
Added: cs/benchmark/cigma/trunk/src/hold/CompareCmd.h
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/CompareCmd.h (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/CompareCmd.h 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,61 @@
+#ifndef __COMPARE_CMD_H__
+#define __COMPARE_CMD_H__
+
+#include <string>
+#include "Command.h"
+#include "MeshPart.h"
+#include "Field.h"
+#include "FE_Field.h"
+#include "QuadratureRule.h"
+#include "ResidualField.h"
+#include "Writer.h"
+#include "Reader.h"
+#include "MeshIO.h"
+#include "QuadratureIO.h"
+#include "FieldIO.h"
+
+namespace cigma
+{
+ class CompareCmd;
+}
+
+
+/*
+ * Callback object for `cigma compare [args ...]'
+ */
+class cigma::CompareCmd : public Command
+{
+public:
+ CompareCmd();
+ ~CompareCmd();
+
+public:
+ void setupOptions(AnyOption *opt);
+ void configure(AnyOption *opt);
+ int run();
+
+public:
+
+ QuadraturePoints *quadrature;
+ MeshPart *mesh;
+ QuadratureRule *qr;
+
+ FE_Field *field_a;
+ FE_Field *field_b;
+
+ ResidualField *residuals;
+
+public:
+ MeshIO meshIO;
+ QuadratureIO quadratureIO;
+ FieldIO firstIO;
+ FieldIO secondIO;
+ std::string output_path;
+
+public:
+ bool verbose;
+ int output_frequency;
+};
+
+
+#endif
Added: cs/benchmark/cigma/trunk/src/hold/EvalCmd.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/EvalCmd.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/EvalCmd.cpp 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,244 @@
+#include <iostream>
+#include <cassert>
+#include "EvalCmd.h"
+#include "StringUtils.h"
+
+#include "HdfReader.h"
+#include "VtkReader.h"
+#include "TextReader.h"
+
+#include "HdfWriter.h"
+#include "VtkWriter.h"
+#include "TextWriter.h"
+
+#include "MeshPart.h"
+
+using namespace std;
+using namespace cigma;
+
+// ---------------------------------------------------------------------------
+
+cigma::EvalCmd::EvalCmd()
+{
+ name = "eval";
+ field = 0;
+ points = 0;
+ values = 0;
+}
+
+cigma::EvalCmd::~EvalCmd()
+{
+}
+
+// ---------------------------------------------------------------------------
+
+void cigma::EvalCmd::setupOptions(AnyOption *opt)
+{
+ //cout << "Calling cigma::EvalCmd::setupOptions()" << endl;
+ assert(opt != 0);
+
+ /* setup usage */
+ opt->addUsage("Usage:");
+ opt->addUsage("");
+ opt->addUsage(" cigma eval [args ...]");
+ opt->addUsage(" --points Input points to evaluate on");
+ opt->addUsage(" --field Input field to evaluate");
+ opt->addUsage(" --values Target path for result of evaluations");
+
+ /* setup flags and options */
+ opt->setFlag("help", 'h');
+ opt->setFlag("verbose", 'v');
+
+ opt->setOption("field");
+ opt->setOption("field-mesh");
+ opt->setOption("field-mesh-coordinates");
+ opt->setOption("field-mesh-connectivity");
+
+ opt->setOption("points");
+ opt->setOption("values");
+ return;
+}
+
+
+void cigma::EvalCmd::configure(AnyOption *opt)
+{
+ //cout << "Calling cigma::EvalCmd::configure()" << endl;
+ assert(opt != 0);
+
+ string inputstr;
+ char *in;
+
+
+ // read verbose flag
+ verbose = opt->getFlag("verbose");
+
+
+ /* determine the output values path */
+ in = opt->getValue("values");
+ if (in == 0)
+ {
+ cerr << "eval: Please specify the output option --values" << endl;
+ exit(1);
+ }
+ inputstr = in;
+
+ // determine the extension and instantiate the appropriate writer object
+ string values_ext;
+ parse_dataset_path(inputstr, values_loc, values_file, values_ext);
+ new_writer(&values_writer, values_ext);
+ if (values_writer == 0)
+ {
+ cerr << "eval: Specified a bad extension (" << values_ext << ") for "
+ << "output file" << values_file << endl;
+ exit(1);
+ }
+ if ((values_loc == "") && (values_writer->getType() == Writer::HDF_WRITER))
+ {
+ values_loc = "/values";
+ }
+
+
+ /* determine the input points path */
+ in = opt->getValue("points");
+ if (in == 0)
+ {
+ cerr << "eval: Please specify the input option --points" << endl;
+ exit(1);
+ }
+ inputstr = in;
+
+ // determine the extension and instantiate the appropriate reader object
+ string points_ext;
+ parse_dataset_path(inputstr, points_loc, points_file, points_ext);
+ new_reader(&points_reader, points_ext);
+ if (points_reader == 0)
+ {
+ cerr << "eval: Specified a bad extension (" << points_ext << ") for "
+ << "input file" << points_file << endl;
+ exit(1);
+ }
+ if ((points_loc == "") && (points_reader->getType() == Reader::HDF_READER))
+ {
+ points_loc = "/points";
+ }
+
+
+ /* determine the input field */
+ load_args(opt, &fieldIO, "field");
+ validate_args(&fieldIO, "eval");
+ fieldIO.load();
+ field = fieldIO.field;
+ assert(field != 0);
+ cout << "field path = " << fieldIO.field_path << endl;
+ cout << "field rank = " << field->n_rank() << endl;
+ //cout << (*field) << endl;
+
+ return;
+}
+
+int cigma::EvalCmd::run()
+{
+ //cout << "Calling cigma::EvalCmd::run()" << endl;
+ assert(field != 0);
+ assert(points != 0);
+ assert(values != 0);
+ assert(field->n_dim() == points->n_dim());
+
+ // indices
+ int e,i;
+
+ // dimensions
+ int npts = points->n_points();
+ int valdim = field->n_rank();
+
+ // data
+ double *phi = new double[npts * valdim];
+
+
+ for (i = 0; i < npts; i++)
+ {
+ double *globalPoint = (*points)[i];
+ field->eval(globalPoint, &phi[valdim*i]);
+ }
+
+
+
+ /*
+ for (i = 0; i < npts; i++)
+ {
+ double *globalPoint = (*points)[i];
+ field->meshPart->find_cell(globalPoint, &e);
+ field->meshPart->get_cell_coords(e, field->fe->cell->globverts);
+ field->eval(globalPoint, &phi[valdim*i]);
+ }
+
+ std::cout << "Creating file " << output_filename << std::endl;
+ VtkWriter *writer = new VtkWriter();
+ writer->open(output_filename);
+ writer->write_header();
+ writer->write_point_data("values", phi, npts, valdim);
+ writer->close();
+ //delete writer;
+ */
+
+
+
+ int ierr;
+
+ cout << "Creating file " << values_file << endl;
+
+ if (values_writer->getType() == Writer::HDF_WRITER)
+ {
+ HdfWriter *writer = static_cast<HdfWriter*>(values_writer);
+ ierr = writer->open(values_file);
+ if (ierr < 0)
+ {
+ cerr << "Error: Could not open (or create) the HDF5 file " << values_file << endl;
+ exit(1);
+ }
+
+ ierr = writer->write_coordinates(values_loc.c_str(), phi, npts, valdim);
+ if (ierr < 0)
+ {
+ cerr << "Error: Could not write values dataset " << values_loc << endl;
+ exit(1);
+ }
+ writer->close();
+ }
+ else if (values_writer->getType() == Writer::TXT_WRITER)
+ {
+ TextWriter *writer = static_cast<TextWriter*>(values_writer);
+ ierr = writer->open(values_file);
+ if (ierr < 0)
+ {
+ cerr << "Error: Could not create output text file " << values_file << endl;
+ exit(1);
+ }
+ writer->write_coordinates(phi, npts, valdim);
+ }
+ else if (values_writer->getType() == Writer::VTK_WRITER)
+ {
+ VtkWriter *writer = static_cast<VtkWriter*>(writer);
+ ierr = vtkWriter->open(values_file);
+ if (ierr < 0)
+ {
+ cerr << "Error: Could not create output VTK file " << values_file << endl;
+ exit(1);
+ }
+ writer->write_header();
+ //writer->write_points(...);
+ writer->write_point_data("values", phi, npts, valdim);
+ writer->close();
+ }
+ else
+ {
+ /* this should be unreachable */
+ cerr << "Fatal Error: Unsupported extension in output filename?" << endl;
+ return 1;
+ }
+
+ // XXX: wrap this guy inside an auto_ptr
+ delete [] phi;
+
+ return 0;
+}
Added: cs/benchmark/cigma/trunk/src/hold/EvalCmd.h
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/EvalCmd.h (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/EvalCmd.h 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,53 @@
+#ifndef __EVAL_CMD_H__
+#define __EVAL_CMD_H__
+
+#include "Command.h"
+
+#include "Field.h"
+#include "Points.h"
+
+//#include "FieldIO.h"
+//#include "PointsIO.h"
+
+
+namespace cigma
+{
+ class EvalCmd;
+}
+
+
+/**
+ * @brief Callback object for `cigma eval [args ...]'
+ *
+ */
+class cigma::EvalCmd : public Command
+{
+public:
+ EvalCmd();
+ ~EvalCmd();
+
+public:
+ void setupOptions(AnyOption *opt);
+ void configure(AnyOption *opt);
+ int run();
+
+public:
+ FieldIO fieldIO;
+
+ Reader *points_reader;
+ std::string points_loc;
+ std::string points_file;
+
+ Writer *values_writer;
+ std::string values_loc;
+ std::string values_file;
+
+ bool verbose;
+
+public:
+ Field *field;
+ Points *points;
+ Points *values;
+};
+
+#endif
Added: cs/benchmark/cigma/trunk/src/hold/ExtractCmd.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/ExtractCmd.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/ExtractCmd.cpp 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,284 @@
+#include <iostream>
+#include <cassert>
+#include "ExtractCmd.h"
+#include "StringUtils.h"
+#include "HdfWriter.h"
+#include "TextWriter.h"
+#include "VtkWriter.h"
+#include "MeshPart.h"
+
+
+using namespace std;
+using namespace cigma;
+
+
+// ---------------------------------------------------------------------------
+
+cigma::ExtractCmd::ExtractCmd()
+{
+ name = "extract";
+ coordsField = 0;
+ writer = 0;
+}
+
+cigma::ExtractCmd::~ExtractCmd()
+{
+}
+
+
+// ---------------------------------------------------------------------------
+
+void cigma::ExtractCmd::setupOptions(AnyOption *opt)
+{
+ //std::cout << "Calling cigma::ExtractCmd::setupOptions()" << std::endl;
+ assert(opt != 0);
+
+ /* setup usage */
+ opt->addUsage("Usage:");
+ opt->addUsage(" cigma extract [options]");
+ opt->addUsage(" --help Display this screen");
+ opt->addUsage(" --mesh Input mesh file");
+ opt->addUsage(" --mesh-coordinates Path to mesh coordinates");
+ opt->addUsage(" --mesh-connectivity Path to mesh connectivity");
+ opt->addUsage(" --rule Quadrature rule");
+ opt->addUsage(" --output Output file");
+
+ /* setup flags and options */
+
+ opt->setFlag("help", 'h');
+ opt->setFlag("verbose", 'v');
+
+ opt->setOption("mesh");
+ opt->setOption("mesh-coordinates");
+ opt->setOption("mesh-connectivity");
+
+ opt->setOption("rule", 'r');
+ opt->setOption("rule-order");
+ opt->setOption("rule-points");
+ opt->setOption("rule-weights");
+
+ opt->setOption("output", 'o');
+}
+
+
+void cigma::ExtractCmd::configure(AnyOption *opt)
+{
+ //std::cout << "Calling cigma::ExtractCmd::configure()" << std::endl;
+ assert(opt != 0);
+
+
+ string inputstr;
+ char *in;
+
+
+ // read verbose flag
+ verbose = opt->getFlag("verbose");
+
+
+ /* determine the output option */
+ in = opt->getValue("output");
+ if (in == 0)
+ {
+ cerr << "extract: Please specify the option --output" << endl;
+ exit(1);
+ }
+ inputstr = in;
+
+
+ // determine the extension and instantiate the appropriate writer object
+ string ext;
+ parse_dataset_path(inputstr, points_path, output_filename, ext);
+ new_writer(&writer, ext);
+ if (writer == 0)
+ {
+ cerr << "extract: Specified bad extension (" << ext << ") for "
+ << "output file" << endl;
+ exit(1);
+ }
+ if ((points_path == "") && (writer->getType() == Writer::HDF_WRITER))
+ {
+ points_path = "/points";
+ }
+
+
+ /* gather up expected command line arguments */
+ load_args(opt, &meshIO, "mesh");
+ load_args(opt, &quadratureIO, "rule");
+
+ /* validate these arguments and complain about missing ones */
+ validate_args(&meshIO, "extract");
+ validate_args(&quadratureIO, "extract");
+
+ /* create empty field object */
+ coordsField = new FE_Field();
+
+ /* load mesh into memory */
+ meshIO.load();
+ if (meshIO.meshPart == 0)
+ {
+ if (!meshIO.has_paths())
+ {
+ cerr << "extract: Missing input option --mesh" << endl;
+ exit(1);
+ }
+ else
+ {
+ cerr << "extract: Invalid mesh options! Could not load mesh." << endl;
+ exit(1);
+ }
+ }
+ coordsField->meshPart = meshIO.meshPart;
+ coordsField->meshPart->set_cell();
+
+
+ /* now, load quadrature rule */
+ quadratureIO.load(coordsField->meshPart->cell);
+ if (quadratureIO.quadrature == 0)
+ {
+ cerr << "extract: Invalid quadrature rule options!" << endl;
+ exit(1);
+ }
+
+
+ /* finite element for interpolation on mesh coordinates */
+ coordsField->fe = new FE();
+ coordsField->fe->set_mesh(meshIO.meshPart);
+ coordsField->fe->set_quadrature_points(quadratureIO.quadrature);
+
+
+ return;
+}
+
+
+int cigma::ExtractCmd::run()
+{
+ //std::cout << "Calling cigma::ExtractCmd::run()" << std::endl;
+ assert(coordsField != 0);
+
+
+ // local data -- destructure field object
+
+ FE *fe = coordsField->fe;
+ assert(fe != 0);
+
+ MeshPart *meshPart = coordsField->meshPart;
+ assert(meshPart != 0);
+
+ QuadraturePoints *quadrature = fe->points;
+ assert(quadrature != 0);
+
+ Cell *cell = fe->meshPart->cell;
+ assert(cell != 0);
+
+
+ // indices
+ int e,q;
+
+
+ // dimensions
+ int nel = meshPart->nel;
+ int nsd = meshPart->nsd;
+ int nq = quadrature->n_points();
+
+ const int cell_nno = cell->n_nodes();
+
+ double *global_points = new double[(nq*nel) * nsd];
+
+
+ if (verbose)
+ {
+ cout << "Extracting "
+ << (nq * nel)
+ << " points from input mesh..."
+ << endl;
+ }
+
+ for (e = 0; e < nel; e++)
+ {
+ double globalCellCoords[cell_nno*3];
+ meshPart->get_cell_coords(e, globalCellCoords);
+ cell->set_global_vertices(globalCellCoords, cell_nno, 3);
+ quadrature->apply_refmap(cell);
+
+ for (q = 0; q < nq; q++)
+ {
+ for (int i = 0; i < nsd; i++)
+ {
+ int n = (nq*nsd)*e + nsd*q + i;
+ global_points[n] = (*quadrature)(q,i);
+ }
+ }
+ /*
+ for (q = 0; q < nq; q++)
+ {
+ std::cout << (*quadrature)(q,0) << " "
+ << (*quadrature)(q,1) << " "
+ << (*quadrature)(q,2) << std::endl;
+ } // */
+ }
+
+
+ int ierr;
+
+ cout << "Creating file " << output_filename << endl;
+
+ if (writer->getType() == Writer::HDF_WRITER)
+ {
+ HdfWriter *hdfWriter = static_cast<HdfWriter*>(writer);
+ ierr = hdfWriter->open(output_filename);
+ if (ierr < 0)
+ {
+ cerr << "Error: Could not open (or create) the HDF5 file " << output_filename << endl;
+ exit(1);
+ }
+
+ ierr = hdfWriter->write_coordinates(points_path.c_str(), global_points, nq*nel, nsd);
+ if (ierr < 0)
+ {
+ cerr << "Error: Could not write points dataset " << points_path << endl;
+ exit(1);
+ }
+
+ hdfWriter->close();
+ }
+ else if (writer->getType() == Writer::TXT_WRITER)
+ {
+ TextWriter *txtWriter = static_cast<TextWriter*>(writer);
+ ierr = txtWriter->open(output_filename);
+ if (ierr < 0)
+ {
+ cerr << "Error: Could not create output text file " << output_filename << endl;
+ exit(1);
+ }
+
+ txtWriter->write_coordinates(global_points, nq * nel, nsd);
+ txtWriter->close();
+ }
+ else if (writer->getType() == Writer::VTK_WRITER)
+ {
+ VtkWriter *vtkWriter = static_cast<VtkWriter*>(writer);
+ ierr = vtkWriter->open(output_filename);
+ if (ierr < 0)
+ {
+ cerr << "Error: Could not create output VTK file " << output_filename << endl;
+ exit(1);
+ }
+
+ vtkWriter->write_header();
+ vtkWriter->write_points(global_points, nq * nel, nsd);
+ vtkWriter->close();
+ }
+ else
+ {
+ /* this should be unreachable */
+ cerr << "Fatal Error: Unsupported extension in output filename?" << endl;
+ return 1;
+ }
+
+
+ // XXX: wrap this guy inside an auto_ptr
+ delete [] global_points;
+
+
+ return 0;
+}
Added: cs/benchmark/cigma/trunk/src/hold/ExtractCmd.h
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/ExtractCmd.h (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/ExtractCmd.h 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,44 @@
+#ifndef __EXTRACT_CMD_H__
+#define __EXTRACT_CMD_H__
+
+#include "Command.h"
+#include "FE_Field.h"
+#include "MeshIO.h"
+#include "QuadratureIO.h"
+#include "Writer.h"
+
+namespace cigma
+{
+ class ExtractCmd;
+}
+
+
+/**
+ * @brief Callback object for `cigma extract [args ...]'
+ *
+ */
+class cigma::ExtractCmd : public Command
+{
+public:
+ ExtractCmd();
+ ~ExtractCmd();
+
+public:
+ void setupOptions(AnyOption *opt);
+ void configure(AnyOption *opt);
+ int run();
+
+public:
+ MeshIO meshIO;
+ QuadratureIO quadratureIO;
+ std::string output_filename;
+ std::string points_path;
+ bool verbose;
+
+public:
+ FE_Field *coordsField;
+ Writer *writer;
+
+};
+
+#endif
Added: cs/benchmark/cigma/trunk/src/hold/FieldIO.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/FieldIO.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/FieldIO.cpp 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,171 @@
+#include <iostream>
+#include "FieldIO.h"
+#include "AnnLocator.h"
+#include "StringUtils.h"
+#include "Misc.h"
+
+using namespace std;
+using namespace cigma;
+
+
+// ---------------------------------------------------------------------------
+
+
+void load_args(AnyOption *opt, FieldIO *fieldIO, const char *opt_prefix)
+{
+ assert(opt != 0);
+ assert(fieldIO != 0);
+
+ // remember the original option name
+ fieldIO->field_option = "--";
+ fieldIO->field_option += opt_prefix;
+
+ // deduce all implied command line options
+ char *in;
+ string field_name = opt_prefix;
+ string mesh_name = field_name + "-mesh";
+
+
+ in = opt->getValue(field_name.c_str());
+ if (in != 0)
+ {
+ fieldIO->field_path = in;
+ }
+
+ load_args(opt, &(fieldIO->meshIO), mesh_name.c_str());
+
+}
+
+void validate_args(FieldIO *fieldIO, const char *cmd_name)
+{
+ assert(fieldIO != 0);
+ assert(fieldIO->field_option != "");
+
+ //
+ // Check for incompatible/inconsistent options
+ //
+
+ if (fieldIO->field_path == "")
+ {
+ cerr << cmd_name
+ << ": Please specify the option "
+ << fieldIO->field_option
+ << endl;
+ exit(1);
+ }
+
+ validate_args(&(fieldIO->meshIO), cmd_name);
+
+ // XXX: thinking about it...maybe Quadrature should be associated
+ // with a mesh, not with FE_Field. Of course, the field->fe object should
+ // only be initialized after we know which mesh to associate with
+ //validate_args(&(fieldIO->quadratureIO), cmd_name);
+
+
+}
+
+// ---------------------------------------------------------------------------
+
+
+FieldIO::FieldIO()
+{
+ field = 0;
+}
+
+
+FieldIO::~FieldIO()
+{
+ if (field != 0)
+ {
+ // XXX: traverse field structure and delete everything
+ }
+}
+
+
+// ---------------------------------------------------------------------------
+
+
+void FieldIO::load()
+{
+ //cout << "Calling FieldIO::load()" << endl;
+
+ int dofs_nno, dofs_valdim;
+ double *dofs;
+
+ string dofs_loc, dofs_file, dofs_ext;
+
+ if (field_path != "")
+ {
+ //
+ // XXX: How should one decide whether we need to instantiate
+ // FE_Field, PointField, or an ExternalField?
+ //
+
+
+ // XXX: for now, assume we are instantiating an FE_Field.
+ // note that we are also assuming that dofs_loc points to a dataset.
+ // if dofs_loc points to an HDF5 group, we'll need to read its metadata
+ // to determine how the data should be interpreted. even if dofs_loc
+ // is a dataset, we should check whether it has its mesh metadata
+ // properly set. lastly, the FunctionSpace attribute should always
+ // when reading from HDF5 files.
+ parse_dataset_path(field_path, dofs_loc, dofs_file, dofs_ext);
+ new_reader(&reader, dofs_ext);
+ reader->open(dofs_file);
+ reader->get_dataset(dofs_loc.c_str(), &dofs, &dofs_nno, &dofs_valdim);
+ //reader->close();
+
+
+ if (meshIO.mesh_path == "")
+ {
+ // XXX: for now, assume that only --field-mesh is set,
+ // and that it points to a file...i think we also need to
+ // set coords_path and connect_path. lastly, don't forget
+ // about --field-mesh-coordinates and --field-mesh-connectivity.
+ meshIO.mesh_path = dofs_file;
+ }
+
+ meshIO.load();
+ assert(meshIO.meshPart != 0);
+
+
+ field = new FE_Field();
+
+ field->dim = meshIO.meshPart->nsd;
+ field->rank = dofs_valdim;
+
+ field->meshPart = meshIO.meshPart;
+ field->meshPart->set_cell();
+ assert(field->meshPart->cell != 0);
+
+ // XXX: Instantiate Locator only when necessary!
+ // XXX: Decide threshold based on number of elements?
+ // XXX: For now, use ann locator for 3D meshes only (use quadtree locator for 2D meshes?)
+ if (field->meshPart->nel > 1000)
+ {
+ if (field->meshPart->nsd == 3)
+ {
+ AnnLocator *locator = new AnnLocator();
+ field->meshPart->set_locator(locator);
+ }
+ }
+
+ field->dofHandler = new DofHandler();
+ field->dofHandler->set_data(dofs, dofs_nno, dofs_valdim);
+ }
+
+ return;
+}
+
+
+// ---------------------------------------------------------------------------
+
+
+void FieldIO::save()
+{
+ cout << "Calling FieldIO::save()" << endl;
+ assert(field != 0);
+ assert(false);
+}
+
+// ---------------------------------------------------------------------------
Added: cs/benchmark/cigma/trunk/src/hold/FieldIO.h
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/FieldIO.h (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/FieldIO.h 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,34 @@
+#ifndef __FIELD_IO_H__
+#define __FIELD_IO_H__
+
+#include <string>
+#include "MeshIO.h"
+#include "QuadratureIO.h"
+#include "FE_Field.h"
+
+class FieldIO
+{
+public:
+ cigma::Reader *reader;
+ cigma::Writer *writer;
+
+public:
+ std::string field_option;
+ std::string field_path;
+ cigma::FE_Field *field;
+ MeshIO meshIO;
+ QuadratureIO quadratureIO;
+
+public:
+ FieldIO();
+ ~FieldIO();
+ void load();
+ void save();
+};
+
+
+void load_args(AnyOption *opt, FieldIO *fieldIO, const char *opt_prefix);
+void validate_args(FieldIO *fieldIO, const char *cmd_name);
+
+
+#endif
Added: cs/benchmark/cigma/trunk/src/hold/HdfArray.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/HdfArray.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/HdfArray.cpp 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,49 @@
+#include "HdfArray.h"
+#include <cassert>
+#include <cstdlib>
+
+
+// ---------------------------------------------------------------------------
+
+HdfArray::HdfArray()
+{
+ rank = 0;
+ shape = 0;
+ data = 0;
+}
+
+HdfArray::~HdfArray()
+{
+ if (shape != 0) free(shape);
+ //if (data != 0) free(data);
+}
+
+// ---------------------------------------------------------------------------
+
+void HdfArray::initialize(hid_t type_id, int rank, int *shape)
+{
+ int i;
+
+ /* initialize dataspace */
+ this->rank = rank;
+ this->shape = (hsize_t *)malloc(rank * sizeof(hsize_t));
+ if (shape != NULL)
+ {
+ for (i = 0; i < rank; i++)
+ {
+ this->shape[i] = shape[i];
+ }
+ }
+ else
+ {
+ for (i = 0; i < rank; i++)
+ {
+ this->shape[i] = 0;
+ }
+ }
+
+ /* initialize data type */
+ this->type_id = type_id;
+
+}
+
Added: cs/benchmark/cigma/trunk/src/hold/HdfArray.h
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/HdfArray.h (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/HdfArray.h 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,47 @@
+#ifndef __HDF_ARRAY_H__
+#define __HDF_ARRAY_H__
+
+#include "hdf5.h"
+
+
+class HdfArray
+{
+public:
+ HdfArray();
+ ~HdfArray();
+
+
+public:
+ void set_shape(hid_t type_id, int rank, int *shape);
+
+public:
+ int create(hid_t loc_id, const char *name, const char *title);
+ int open(hid_t loc_id, const char *name);
+
+public:
+ int read(hid_t loc_id, const char *name);
+ int write(hid_t loc_id, const char *name);
+
+public:
+ int slice_read(int start, int end);
+ int slice_write(int start, int end);
+
+public:
+ int rank;
+ hsize_t *shape;
+ hid_t type_id;
+ void *data;
+
+public:
+ hid_t dataset_id;
+ hsize_t *offset;
+ hsize_t *stride;
+ hsize_t *count;
+ hsize_t *block;
+
+private:
+ HdfArray(const HdfArray &);
+
+};
+
+#endif
Added: cs/benchmark/cigma/trunk/src/hold/HdfGroup.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/HdfGroup.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/HdfGroup.cpp 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,34 @@
+#include "HdfGroup.h"
+#include <cassert>
+
+// ---------------------------------------------------------------------------
+
+HdfGroup::HdfGroup(hid_t loc_id, const char *name)
+{
+ herr_t status;
+
+ group_id = H5Gcreate(loc_id, name, 0);
+ if (group_id < 0)
+ {
+
+ }
+}
+
+HdfGroup::~HdfGroup()
+{
+ herr_t status;
+
+ if (group_id < 0)
+ {
+ status = H5Gclose(group_id);
+ }
+
+ group_id = -1;
+}
+
+// ---------------------------------------------------------------------------
+
+
+
+// ---------------------------------------------------------------------------
+
Added: cs/benchmark/cigma/trunk/src/hold/HdfGroup.h
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/HdfGroup.h (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/HdfGroup.h 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,33 @@
+#ifndef __HDF_GROUP_H__
+#define __HDF_GROUP_H__
+
+#include "hdf5.h"
+
+class HdfGroup
+{
+public:
+ HdfGroup(hid_t loc_id, const char *name);
+ ~HdfGroup();
+
+private:
+ HdfGroup(const HdfGroup &);
+
+public:
+ bool exists();
+
+
+public:
+ hid_t group_id;
+
+};
+
+// ---------------------------------------------------------------------------
+
+inline bool HdfGroup::exists()
+{
+ return (group_id >= 0);
+}
+
+// ---------------------------------------------------------------------------
+
+#endif
Added: cs/benchmark/cigma/trunk/src/hold/HdfReader.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/HdfReader.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/HdfReader.cpp 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,127 @@
+#include "HdfReader.h"
+#include "h5io.h"
+#include <cstdlib>
+#include <cassert>
+
+// ---------------------------------------------------------------------------
+
+cigma::HdfReader::HdfReader()
+{
+ file_id = -1;
+}
+
+cigma::HdfReader::~HdfReader()
+{
+}
+
+// ---------------------------------------------------------------------------
+
+
+int cigma::HdfReader::
+open(std::string filename)
+{
+ file_id = h5io_file_open(filename.c_str(), "r");
+ assert(file_id >= 0);
+ return 0; // XXX: change return value instead of using assert
+}
+
+
+void cigma::HdfReader::
+close()
+{
+ herr_t status = H5Fclose(file_id);
+ if (status < 0)
+ {
+ // XXX: emit warning?
+ }
+}
+
+// ---------------------------------------------------------------------------
+
+void cigma::HdfReader::
+get_dataset(const char *loc, double **data, int *num, int *dim)
+{
+ int rank;
+ hid_t type_id;
+ hid_t dataset_id;
+ herr_t status;
+ int ierr = -1;
+
+ dataset_id = h5io_dset_open(file_id, loc, &type_id, &rank, NULL, NULL);
+ assert(H5Tget_class(type_id) == H5T_FLOAT);
+ assert(dataset_id >= 0);
+ assert((rank == 1) || (rank == 2));
+ status = H5Dclose(dataset_id);
+ assert(status >= 0); // XXX: emit warning?
+
+ if (rank == 2)
+ {
+ ierr = h5io_dset_read2(file_id, loc, H5T_NATIVE_DOUBLE,
+ (void **)data, num, dim);
+ }
+ else if (rank == 1)
+ {
+ ierr = h5io_dset_read1(file_id, loc, H5T_NATIVE_DOUBLE,
+ (void **)data, num);
+ *dim = 1;
+ }
+ assert(ierr >= 0);
+}
+
+void cigma::HdfReader::
+get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd)
+{
+ get_dataset(loc, coordinates, nno, nsd);
+}
+
+void cigma::HdfReader::
+get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs)
+{
+ int rank;
+ hid_t type_id;
+ hid_t connect_id;
+ herr_t status;
+
+ connect_id = h5io_dset_open(file_id, loc, &type_id, &rank, NULL, NULL);
+ assert(H5Tget_class(type_id) == H5T_INTEGER);
+ assert(connect_id >= 0);
+ assert(rank == 2);
+ status = H5Dclose(connect_id);
+ if (status < 0)
+ {
+ // XXX: emit warning
+ }
+
+ int ierr;
+ ierr = h5io_dset_read2(file_id, loc, H5T_NATIVE_INT,
+ (void **)connectivity, nel, ndofs);
+}
+
+
+// ---------------------------------------------------------------------------
+
+void cigma::HdfReader::
+get_mesh(MeshPart *meshPart, const char *loc)
+{
+ // XXX: read coords_loc & connect_loc from metadata
+}
+
+void cigma::HdfReader::
+get_mesh(MeshPart *meshPart, const char *coords_loc, const char *connect_loc)
+{
+ assert(meshPart != 0);
+ get_coordinates(coords_loc, &(meshPart->coords), &(meshPart->nno), &(meshPart->nsd));
+ get_connectivity(connect_loc, &(meshPart->connect), &(meshPart->nel), &(meshPart->nel));
+}
+
+void cigma::HdfReader::
+get_dofs(DofHandler *dofHandler, const char *loc)
+{
+ assert(dofHandler != 0);
+ assert(dofHandler->meshPart != 0);
+ get_dataset(loc, &(dofHandler->dofs), &(dofHandler->nno), &(dofHandler->ndim));
+}
+
+
+
+// ---------------------------------------------------------------------------
Added: cs/benchmark/cigma/trunk/src/hold/HdfReader.h
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/HdfReader.h (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/HdfReader.h 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,53 @@
+#ifndef __HDF_READER_H__
+#define __HDF_READER_H__
+
+#include <string>
+#include "hdf5.h"
+#include "Reader.h"
+
+#include "QuadraturePoints.h"
+#include "MeshPart.h"
+#include "DofHandler.h"
+#include "FE_Field.h"
+
+
+namespace cigma
+{
+ class HdfReader;
+}
+
+
+class cigma::HdfReader : public Reader
+{
+public:
+ HdfReader();
+ ~HdfReader();
+
+public:
+ ReaderType getType() { return HDF_READER; }
+
+public:
+ int open(std::string filename);
+ void close();
+
+public:
+ void get_dataset(const char *loc, double **data, int *num, int *dim);
+ void get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd);
+ void get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs);
+
+public:
+ void get_quadrature(QuadraturePoints *quadrature, const char *loc);
+ void get_quadrature(QuadraturePoints *quadrature, const char *points_loc, const char *weights_loc);
+
+public:
+ void get_field(FE_Field *field, const char *loc);
+ void get_mesh(MeshPart *meshPart, const char *coords_loc, const char *connect_loc);
+ void get_mesh(MeshPart *meshPart, const char *loc);
+ void get_dofs(DofHandler *dofHandler, const char *loc);
+
+public:
+ hid_t file_id;
+};
+
+
+#endif
Added: cs/benchmark/cigma/trunk/src/hold/HdfWriter.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/HdfWriter.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/HdfWriter.cpp 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,71 @@
+#include <cassert>
+#include <cstdlib>
+#include "HdfWriter.h"
+#include "h5io.h"
+
+using namespace std;
+using namespace cigma;
+
+
+// ---------------------------------------------------------------------------
+
+HdfWriter::HdfWriter()
+{
+ file_id = -1;
+}
+
+HdfWriter::~HdfWriter()
+{
+ close();
+}
+
+
+// ---------------------------------------------------------------------------
+
+int HdfWriter::open(string filename)
+{
+ file_id = h5io_file_open(filename.c_str(), "rw+");
+
+ if (file_id < 0)
+ {
+ return -1;
+ }
+ return 0;
+}
+
+void HdfWriter::close()
+{
+ herr_t status;
+ if (file_id != -1)
+ {
+ status = H5Fclose(file_id);
+ }
+ file_id = -1;
+}
+
+
+// ---------------------------------------------------------------------------
+
+int HdfWriter::write_coordinates(const char *loc, double *coordinates, int nno, int nsd)
+{
+ int ierr;
+ ierr = h5io_dset_write2(file_id, loc, "coordinates array", H5T_NATIVE_DOUBLE, coordinates, nno, nsd);
+ return ierr;
+}
+
+int HdfWriter::write_connectivity(const char *loc, int *connectivity, int nel, int ndofs)
+{
+ int ierr;
+ ierr = h5io_dset_write2(file_id, loc, "connectivity array", H5T_NATIVE_INT, connectivity, nel, ndofs);
+ return ierr;
+}
+
+// ---------------------------------------------------------------------------
+
+void cigma::HdfWriter::
+write_field(FE_Field *field)
+{
+ assert(field != 0);
+}
+
+// ---------------------------------------------------------------------------
Added: cs/benchmark/cigma/trunk/src/hold/HdfWriter.h
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/HdfWriter.h (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/HdfWriter.h 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,37 @@
+#ifndef __HDF_WRITER_H__
+#define __HDF_WRITER_H__
+
+#include "Writer.h"
+#include "hdf5.h"
+
+
+namespace cigma
+{
+ class HdfWriter;
+}
+
+
+class cigma::HdfWriter : public Writer
+{
+public:
+ HdfWriter();
+ ~HdfWriter();
+
+public:
+ WriterType getType() { return HDF_WRITER; }
+ int open(std::string filename);
+ void close();
+
+public:
+ int write_coordinates(const char *loc, double *coordinates, int nno, int nsd);
+ int write_connectivity(const char *loc, int *connectivity, int nel, int ndofs);
+
+public:
+ void write_field(FE_Field *field);
+
+public:
+ hid_t file_id;
+};
+
+
+#endif
Added: cs/benchmark/cigma/trunk/src/hold/Makefile.am
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/Makefile.am (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/Makefile.am 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,153 @@
+bin_PROGRAMS = cigma
+
+noinst_LIBRARIES = libh5io.a libann.a libcigma.a
+
+libh5io_a_CPPFLAGS = $(HDF5_INCLUDES)
+
+libann_a_CPPFLAGS = -I$(srcdir)/ann/include
+
+libcigma_a_CPPFLAGS = -Wno-deprecated $(HDF5_INCLUDES) $(PYTHON_INCLUDES) $(VTK_INCLUDES) -I$(srcdir)/h5io -I$(srcdir)/ann/include
+#libcigma_a_LIBADD = libh5io.a libann.a
+
+cigma_CPPFLAGS = -Wno-deprecated $(HDF5_INCLUDES) $(PYTHON_INCLUDES) $(VTK_INCLUDES) -I$(srcdir)/h5io -I$(srcdir)/ann/include
+cigma_LDFLAGS = -I$(srcdir)/h5io -I$(srcdir)/ann/include
+cigma_LDADD = $(HDF5_LIBS) $(HDF5_LDFLAGS) -L$(srcdir) $(PYTHON_LDFLAGS) $(PYTHON_LIBS) libcigma.a libh5io.a libann.a
+
+
+##############################################################################
+
+cigma_SOURCES = \
+ cigma.cpp
+
+libh5io_a_SOURCES = \
+ h5io/h5array.c \
+ h5io/h5attr.c \
+ h5io/h5dset.c \
+ h5io/h5file.c \
+ h5io/h5group.c \
+ h5io/h5io.h \
+ h5io/split.c \
+ h5io/split.h
+
+libann_a_SOURCES = \
+ ann/include/ANN/ANN.h \
+ ann/include/ANN/ANNperf.h \
+ ann/include/ANN/ANNx.h \
+ ann/src/ANN.cpp \
+ ann/src/bd_fix_rad_search.cpp \
+ ann/src/bd_pr_search.cpp \
+ ann/src/bd_search.cpp \
+ ann/src/bd_tree.cpp \
+ ann/src/bd_tree.h \
+ ann/src/brute.cpp \
+ ann/src/kd_dump.cpp \
+ ann/src/kd_fix_rad_search.cpp \
+ ann/src/kd_fix_rad_search.h \
+ ann/src/kd_pr_search.cpp \
+ ann/src/kd_pr_search.h \
+ ann/src/kd_search.cpp \
+ ann/src/kd_search.h \
+ ann/src/kd_split.cpp \
+ ann/src/kd_split.h \
+ ann/src/kd_tree.cpp \
+ ann/src/kd_tree.h \
+ ann/src/kd_util.cpp \
+ ann/src/kd_util.h \
+ ann/src/perf.cpp \
+ ann/src/pr_queue.h \
+ ann/src/pr_queue_k.h
+
+libcigma_a_SOURCES = \
+ AnnLocator.cpp \
+ AnnLocator.h \
+ AnyOption.cpp \
+ AnyOption.h \
+ Cell.cpp \
+ Cell.h \
+ Command.cpp \
+ Command.h \
+ CommandSet.cpp \
+ CommandSet.h \
+ CompareCmd.cpp \
+ CompareCmd.h \
+ convert.h \
+ CubeCmd.h \
+ CubeCmd.cpp \
+ CubeMeshPart.cpp \
+ CubeMeshPart.h \
+ DofHandler.cpp \
+ DofHandler.h \
+ ExtractCmd.h \
+ ExtractCmd.cpp \
+ FE.cpp \
+ FE_Field.cpp \
+ FE_Field.h \
+ FE.h \
+ FiatProxy.cpp \
+ FiatProxy.h \
+ Field.cpp \
+ Field.h \
+ FieldIO.cpp \
+ FieldIO.h \
+ HdfReader.cpp \
+ HdfReader.h \
+ HdfWriter.cpp \
+ HdfWriter.h \
+ HelpCmd.cpp \
+ HelpCmd.h \
+ Hex.cpp \
+ Hex.h \
+ ListCmd.cpp \
+ ListCmd.h \
+ Locator.cpp \
+ Locator.h \
+ MeshIO.cpp \
+ MeshIO.h \
+ MeshPart.cpp \
+ MeshPart.h \
+ Misc.cpp \
+ Misc.h \
+ Numeric.cpp \
+ Numeric.h \
+ Points.cpp \
+ Points.h \
+ PointField.cpp \
+ PointField.h \
+ Quad.cpp \
+ Quad.h \
+ QuadratureIO.cpp \
+ QuadratureIO.h \
+ QuadraturePoints.cpp \
+ QuadraturePoints.h \
+ QuadratureRule.cpp \
+ QuadratureRule.h \
+ Reader.cpp \
+ Reader.h \
+ ResidualField.cpp \
+ ResidualField.h \
+ SkelCmd.cpp \
+ SkelCmd.h \
+ StringUtils.cpp \
+ StringUtils.h \
+ Tet.cpp \
+ Tet.h \
+ TextReader.cpp \
+ TextReader.h \
+ TextWriter.cpp \
+ TextWriter.h \
+ Timer.h \
+ Tri.cpp \
+ Tri.h \
+ VtkList.cpp \
+ VtkList.h \
+ VtkReader.cpp \
+ VtkReader.h \
+ VtkXmlReader.cpp \
+ VtkXmlReader.h \
+ VtkWriter.cpp \
+ VtkWriter.h \
+ Writer.cpp \
+ Writer.h \
+ ZeroField.cpp \
+ ZeroField.h
+
Added: cs/benchmark/cigma/trunk/src/hold/MeshIO.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/MeshIO.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/MeshIO.cpp 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,443 @@
+#include <iostream>
+#include <cstdlib>
+#include "MeshIO.h"
+#include "StringUtils.h"
+
+using namespace std;
+using namespace cigma;
+
+// ---------------------------------------------------------------------------
+
+void load_args(AnyOption *opt, MeshIO *meshIO, const char *opt_prefix)
+{
+ assert(opt != 0);
+ assert(meshIO != 0);
+
+ char *in;
+ string optstr;
+ string mesh_name = opt_prefix;
+
+ in = opt->getValue(mesh_name.c_str());
+ if (in != 0)
+ {
+ meshIO->mesh_path = in;
+ }
+
+ optstr = mesh_name + "-coordinates";
+ in = opt->getValue(optstr.c_str());
+ if (in != 0)
+ {
+ meshIO->coords_path = in;
+ }
+
+ optstr = mesh_name + "-connectivity";
+ in = opt->getValue(optstr.c_str());
+ if (in != 0)
+ {
+ meshIO->connect_path = in;
+ }
+}
+
+void validate_args(MeshIO *meshIO, const char *cmd_name)
+{
+ assert(meshIO != 0);
+
+ if (meshIO->has_paths())
+ {
+ if (meshIO->mesh_path == "")
+ {
+ if (meshIO->coords_path == "")
+ {
+ cerr << cmd_name
+ << ": Detected missing option --mesh-coordinates"
+ << endl;
+
+ exit(1);
+ }
+ if (meshIO->connect_path == "")
+ {
+ cerr << cmd_name
+ << ": Detected missing option --mesh-connectivity"
+ << endl;
+
+ exit(1);
+ }
+ }
+ }
+}
+
+
+// ---------------------------------------------------------------------------
+
+MeshIO::MeshIO()
+{
+ meshPart = 0;
+ reader = coords_reader = connect_reader = 0;
+ writer = coords_writer = connect_writer = 0;
+}
+
+MeshIO::~MeshIO()
+{
+ if (meshPart != 0)
+ {
+ // XXX: traverse meshPart structure and delete everything
+ }
+}
+
+// ---------------------------------------------------------------------------
+
+
+bool MeshIO::has_paths()
+{
+ return (mesh_path != "") || (coords_path != "") || (connect_path != "");
+}
+
+bool MeshIO::has_coords_path()
+{
+ return (mesh_path != "") || (coords_path != "");
+}
+
+bool MeshIO::has_connect_path()
+{
+ return (mesh_path != "") || (connect_path != "");
+}
+
+bool MeshIO::has_valid_paths()
+{
+ // empty state is also valid
+ if (!has_paths())
+ {
+ return true;
+ }
+
+ // require both coords and connect paths if mesh path is absent
+ return (mesh_path != "") || ((coords_path != "") && (connect_path != ""));
+}
+
+// ---------------------------------------------------------------------------
+
+bool MeshIO::prepare()
+{
+ // This function initializes the required variables
+ // for use in MeshIO::load()
+ //
+ // precondition: all path variables need to be set
+ // appropriately by calling load_args()
+ //
+ // postcondition: if successful, location & file variables
+ // will be ready for use in load()
+ //
+
+ if (!(has_coords_path() && has_connect_path()))
+ {
+ // not enough information available to set location variables.
+ // top level program should use validate_args() to detect this
+ // condition ahead of time and bail out
+ return false;
+ }
+
+ //
+ // Note that each of coords_path and connect_path will override
+ // the corresponding component in mesh_path.
+ //
+ if (coords_path != "")
+ {
+ parse_dataset_path(coords_path, coords_loc, coords_file, coords_ext);
+ }
+ if (connect_path != "")
+ {
+ parse_dataset_path(connect_path, connect_loc, connect_file, connect_ext);
+ }
+ if (mesh_path != "")
+ {
+ parse_dataset_path(mesh_path, mesh_loc, mesh_file, mesh_ext);
+
+ if (coords_path == "")
+ {
+ // no coords path specified...
+ // read coords from mesh file
+ // coords_loc will have to be determined later in load()
+ coords_loc = "";
+ coords_file = mesh_file;
+ coords_ext = mesh_ext;
+ }
+ else
+ {
+ if (coords_loc == "")
+ {
+ // reinterpret the meaning of --mesh-coordinates
+ coords_loc = coords_file;
+ coords_file = mesh_file;
+ coords_ext = mesh_ext;
+ }
+ }
+
+ if (connect_path == "")
+ {
+ // no connect path specified...
+ // read connectivity from mesh file
+ // connect_loc will have to be determined later in load()
+ connect_loc = "";
+ connect_file = mesh_file;
+ connect_ext = mesh_ext;
+ }
+ else
+ {
+ if (connect_loc == "")
+ {
+ // reinterpret the meaning of --mesh-connectivity
+ connect_loc = connect_file;
+ connect_file = mesh_file;
+ connect_ext = mesh_ext;
+ }
+ }
+ }
+
+ return true;
+}
+
+void MeshIO::load()
+{
+ //
+ // The sole purpose of this function is to initialize
+ // the meshPart member variable. Thus, failure can be detected
+ // by checking whether meshPart is 0.
+ //
+
+
+ //cout << "Calling MeshIO::load()" << endl;
+
+
+ // Skip this method if paths are empty.
+ // This behavior is necessary because the --mesh option
+ // is not required on the command line
+ if (!has_paths())
+ {
+ return;
+ }
+
+ // bail out if prepare() fails
+ bool ready = prepare();
+ if (!ready)
+ {
+ return;
+ }
+
+
+ //
+ // Local declarations
+ //
+ int ierr;
+
+ int nno, nsd;
+ double *coords = 0;
+
+ int nel, ndofs;
+ int *connect = 0;
+
+ nno = nsd = 0;
+ nel = ndofs = 0;
+
+
+ // To implement at this point:
+ // XXX: check whether files exist if we are reading
+ // XXX: emit warnings or bail out if we are writing an existing file?
+ // XXX: check file's magic number and overwrite the extension if it doesn't match
+ //
+
+
+ //
+ // Instantiate readers
+ //
+ if (mesh_file != "")
+ {
+ new_reader(&reader, mesh_ext);
+ }
+ if (coords_file != "")
+ {
+ if (coords_file == mesh_file)
+ {
+ coords_reader = reader;
+ }
+ else
+ {
+ new_reader(&coords_reader, coords_ext);
+ }
+ }
+ if (connect_file != "")
+ {
+ if (connect_file == mesh_file)
+ {
+ connect_reader = reader;
+ }
+ else
+ {
+ new_reader(&connect_reader, connect_ext);
+ }
+ }
+ assert(coords_reader != 0);
+ assert(connect_reader != 0);
+
+
+ //
+ // Open files
+ //
+ if (reader != 0)
+ {
+ ierr = reader->open(mesh_file);
+ assert(ierr == 0);
+ }
+ if (coords_reader != reader)
+ {
+ ierr = coords_reader->open(coords_file);
+ assert(ierr == 0);
+ }
+ if (connect_reader != reader)
+ {
+ ierr = connect_reader->open(connect_file);
+ assert(ierr == 0);
+ }
+
+
+ //
+ // Last chance to determine coords_loc and connect_loc
+ // XXX: for HDF5 files, check dataset or attribute first?
+ // XXX: for HDF5 files, if mesh_loc is empty, does it *have* to be "/"?
+ //
+ if (coords_loc == "")
+ {
+ if (coords_reader->getType() == Reader::HDF_READER)
+ {
+ coords_loc = mesh_loc + "/coordinates";
+ // XXX: if at this point the coords_loc dataset doesn't exist,
+ // determine the following before bailing out:
+ // * assert that mesh_loc is a group
+ // * the group mesh_loc has a CoordLocation string attribute
+ // * load the CoordLocation string attribute into coords_loc
+ }
+ else if (coords_reader->getType() == Reader::TXT_READER)
+ {
+ // XXX: not yet implemented!
+ // XXX: bail out with a better message
+ // Here we would need support for reading multiple datasets
+ // from a single text file (low priority since we can use
+ // HDF5 for this very purpose).
+ assert(false);
+ }
+ else if (coords_reader->getType() == Reader::VTK_READER)
+ {
+ // XXX: it is not necessary to set coords_loc for VTK files
+ // since VtkReader::get_coordinates() ignores its location argument
+ }
+
+ }
+ if (connect_loc == "")
+ {
+ if (connect_reader->getType() == Reader::HDF_READER)
+ {
+ connect_loc = mesh_loc + "/connectivity";
+ // XXX: if at this point the connect_loc dataset doesn't exist,
+ // determine the following before bailing out:
+ // * assert that mesh_loc is a group
+ // * the group mesh_loc has a ConnectLocation string attribute
+ // * load the ConnectLocation string attribute into coords_loc
+ }
+ }
+
+ //
+ // Read datasets
+ // XXX: change these functions to return an error code
+ // XXX: need to check return values for failure
+ coords_reader->get_coordinates(coords_loc.c_str(), &coords, &nno, &nsd);
+ connect_reader->get_connectivity(connect_loc.c_str(), &connect, &nel, &ndofs);
+
+ /*
+ cout << "*** "
+ << "nno = " << nno << ", "
+ << "nsd = " << nsd << ", "
+ << "nel = " << nel << ", "
+ << "ndofs = " << ndofs << endl;
+ // */
+
+
+
+ /* XXX: use auto_ptr for the local readers, so we can throw exceptions
+ if (coords_path != "")
+ {
+ Reader *coords_reader;
+ parse_dataset_path(coords_path, coords_loc, coords_file, coords_ext);
+ new_reader(&coords_reader, coords_ext);
+ coords_reader->open(coords_file);
+ coords_reader->get_coordinates(coords_loc.c_str(), &coords, &nno, &nsd);
+ //coords_reader->close();
+ }
+
+ if (connect_path != "")
+ {
+ Reader *connect_reader;
+ parse_dataset_path(connect_path, connect_loc, connect_file, connect_ext);
+ new_reader(&connect_reader, connect_ext);
+ connect_reader->open(connect_file);
+ connect_reader->get_connectivity(connect_loc.c_str(), &connect, &nel, &ndofs);
+ //connect_reader->close();
+ }
+
+ if ((mesh_path != "") && ((coords == 0) || (connect == 0)))
+ {
+ Reader *mesh_reader;
+ parse_dataset_path(mesh_path, mesh_loc, mesh_file, mesh_ext);
+ new_reader(&mesh_reader, mesh_ext);
+ mesh_reader->open(mesh_file);
+ if (coords == 0)
+ {
+ if (mesh_reader->getType() == Reader::HDF_READER)
+ {
+ coords_loc = mesh_loc + "/coordinates";
+ }
+ mesh_reader->get_coordinates(coords_loc.c_str(), &coords, &nno, &nsd);
+ }
+ if (connect == 0)
+ {
+ if (mesh_reader->getType() == Reader::HDF_READER)
+ {
+ connect_loc = mesh_loc + "/connectivity";
+ }
+ mesh_reader->get_connectivity(connect_loc.c_str(), &connect, &nel, &ndofs);
+ }
+ //mesh_reader->close();
+ } // */
+
+
+
+
+ if ((coords != 0) && (connect != 0))
+ {
+ meshPart = new MeshPart();
+
+ meshPart->nno = nno;
+ meshPart->nsd = nsd;
+ meshPart->coords = coords;
+
+ meshPart->nel = nel;
+ meshPart->ndofs = ndofs;
+ meshPart->connect = connect;
+ }
+ else
+ {
+ if (coords == 0)
+ {
+ cerr << "MeshIO::load() error: "
+ << "Could not find mesh coordinates"
+ << endl;
+ }
+ if (connect == 0)
+ {
+ cerr << "MeshIO::load() error: "
+ << "Could not find mesh connectivity"
+ << endl;
+ }
+ }
+
+}
+
+// ---------------------------------------------------------------------------
Added: cs/benchmark/cigma/trunk/src/hold/MeshIO.h
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/MeshIO.h (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/MeshIO.h 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,53 @@
+#ifndef __MESH_IO_H__
+#define __MESH_IO_H__
+
+#include <string>
+#include "AnyOption.h"
+#include "Reader.h"
+#include "Writer.h"
+#include "MeshPart.h"
+
+class MeshIO
+{
+public:
+ std::string mesh_path;
+ std::string coords_path;
+ std::string connect_path;
+
+public:
+ std::string mesh_loc, mesh_file, mesh_ext;
+ std::string coords_loc, coords_file, coords_ext;
+ std::string connect_loc, connect_file, connect_ext;
+ bool prepare();
+
+public:
+ cigma::Reader *reader;
+ cigma::Reader *coords_reader;
+ cigma::Reader *connect_reader;
+ cigma::Writer *writer;
+ cigma::Writer *coords_writer;
+ cigma::Writer *connect_writer;
+ cigma::MeshPart *meshPart;
+
+public:
+ MeshIO();
+ ~MeshIO();
+
+public:
+ void load();
+ void save();
+
+public:
+ bool has_paths();
+ bool has_coords_path();
+ bool has_connect_path();
+ bool has_valid_paths();
+};
+
+
+
+void load_args(AnyOption *opt, MeshIO *meshIO, const char *opt_prefix);
+void validate_args(MeshIO *meshIO, const char *cmd_name);
+
+
+#endif
Added: cs/benchmark/cigma/trunk/src/hold/QuadratureIO.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/QuadratureIO.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/QuadratureIO.cpp 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,441 @@
+#include <iostream>
+#include "QuadratureIO.h"
+#include "StringUtils.h"
+
+using namespace std;
+using namespace cigma;
+
+// ---------------------------------------------------------------------------
+
+void load_args(AnyOption *opt, QuadratureIO *quadratureIO, const char *opt_prefix)
+{
+ assert(opt != 0);
+ assert(quadratureIO != 0);
+
+ char *in;
+
+ in = opt->getValue("rule-order");
+ if (in != 0)
+ {
+ quadratureIO->quadrature_order = in;
+ }
+
+ in = opt->getValue("rule");
+ if (in != 0)
+ {
+ quadratureIO->quadrature_path = in;
+ }
+
+ in = opt->getValue("rule-points");
+ if (in != 0)
+ {
+ quadratureIO->points_path = in;
+ }
+
+ in = opt->getValue("rule-weights");
+ if (in != 0)
+ {
+ quadratureIO->weights_path = in;
+ }
+}
+
+void validate_args(QuadratureIO *quadratureIO, const char *cmd_name)
+{
+ assert(quadratureIO != 0);
+
+ //
+ // Check for incompatible/inconsistent options
+ //
+
+ if ((quadratureIO->points_path == "") && (quadratureIO->weights_path != ""))
+ {
+ cerr << cmd_name
+ << ": You also need to set the option --rule-points"
+ << endl;
+
+ exit(1);
+ }
+
+ if ((quadratureIO->weights_path == "") && (quadratureIO->points_path != ""))
+ {
+ cerr << cmd_name
+ << ": You also need to set the option --rule-weights"
+ << endl;
+
+ exit(1);
+ }
+
+ if ((quadratureIO->weights_path != "") && (quadratureIO->points_path != ""))
+ {
+ if (quadratureIO->quadrature_path != "")
+ {
+ cerr << cmd_name
+ << ": Using explicit points and weights (--rule not necessary)"
+ << endl;
+
+ exit(1);
+ }
+ }
+
+ if ((quadratureIO->quadrature_order != "")
+ && ((quadratureIO->quadrature_path != "")
+ || ((quadratureIO->points_path != "") && (quadratureIO->weights_path != ""))))
+ {
+ cerr << cmd_name
+ << ": Cannot use --rule-order with an explicit quadrature rule"
+ << endl;
+
+ exit(1);
+ }
+
+
+}
+
+// ---------------------------------------------------------------------------
+
+QuadratureIO::QuadratureIO()
+{
+ quadrature = 0;
+}
+
+QuadratureIO::~QuadratureIO()
+{
+ if (quadrature != 0)
+ {
+ delete quadrature;
+ }
+}
+
+
+// ---------------------------------------------------------------------------
+
+void QuadratureIO::load(cigma::Cell *cell)
+{
+ //cout << "Calling QuadratureIO::load()" << endl;
+
+ assert(cell != 0);
+
+ int i,j;
+ int ierr;
+
+
+ //
+ // XXX: move these default rules into the appropriate Cell subclasses
+ //
+
+ // tri_qr(5)
+ const int tri_nno = 9;
+ const int tri_celldim = 2;
+ const int tri_nsd = 2; //XXX
+ double tri_qpts[tri_nno * tri_celldim] = {
+ -0.79456469, -0.82282408,
+ -0.86689186, -0.18106627,
+ -0.95213774, 0.57531892,
+ -0.08858796, -0.82282408,
+ -0.40946686, -0.18106627,
+ -0.78765946, 0.57531892,
+ 0.61738877, -0.82282408,
+ 0.04795814, -0.18106627,
+ -0.62318119, 0.57531892
+ };
+ double tri_qwts[tri_nno] = {
+ 0.22325768, 0.25471234, 0.07758553, 0.35721229, 0.40753974,
+ 0.12413685, 0.22325768, 0.25471234, 0.07758553
+ };
+ for (i = 0; i < tri_nno; i++)
+ {
+ for (j = 0; j < tri_celldim; j++)
+ {
+ int n = tri_celldim * i + j;
+ tri_qpts[n] += 1;
+ tri_qpts[n] *= 0.5;
+ }
+ }
+
+
+ // quad_qr(7)
+ const int quad_nno = 16;
+ const int quad_celldim = 2;
+ const int quad_nsd = 2; //XXX
+ double quad_qpts[quad_nno * quad_celldim] = {
+ -0.86113631, -0.86113631,
+ -0.33998104, -0.86113631,
+ 0.33998104, -0.86113631,
+ 0.86113631, -0.86113631,
+ 0.86113631, -0.33998104,
+ 0.86113631, 0.33998104,
+ 0.86113631, 0.86113631,
+ 0.33998104, 0.86113631,
+ -0.33998104, 0.86113631,
+ -0.86113631, 0.86113631,
+ -0.86113631, 0.33998104,
+ -0.86113631, -0.33998104,
+ -0.33998104, -0.33998104,
+ 0.33998104, -0.33998104,
+ -0.33998104, 0.33998104,
+ 0.33998104, 0.33998104
+ };
+ double quad_qwts[quad_nno] = {
+ 0.12100299, 0.22685185, 0.22685185, 0.12100299, 0.22685185,
+ 0.22685185, 0.12100299, 0.22685185, 0.22685185, 0.12100299,
+ 0.22685185, 0.22685185, 0.4252933 , 0.4252933 , 0.4252933 ,
+ 0.4252933
+ };
+
+
+ // tet_qr(3)
+ const int tet_nno = 8;
+ const int tet_celldim = 3;
+ const int tet_nsd = 3;
+ double tet_qpts[tet_nno * tet_celldim] = {
+ -0.68663473, -0.72789005, -0.75497035,
+ -0.83720867, -0.85864055, 0.08830369,
+ -0.86832263, 0.13186633, -0.75497035,
+ -0.93159441, -0.4120024 , 0.08830369,
+ 0.16949513, -0.72789005, -0.75497035,
+ -0.39245447, -0.85864055, 0.08830369,
+ -0.50857335, 0.13186633, -0.75497035,
+ -0.74470688, -0.4120024 , 0.08830369
+ };
+ double tet_qwts[tet_nno] = {
+ 0.29583885, 0.12821632, 0.16925605, 0.07335544, 0.29583885,
+ 0.12821632, 0.16925605, 0.07335544
+ };
+ for (i = 0; i < tet_nno; i++)
+ {
+ for (j = 0; j < tet_celldim; j++)
+ {
+ int n = tet_celldim * i + j;
+ tet_qpts[n] += 1;
+ tet_qpts[n] *= 0.5;
+ }
+ }
+
+ // hex_qr(3)
+ const int hex_nno = 8;
+ const int hex_celldim = 3;
+ const int hex_nsd = 3;
+ double hex_qpts[hex_nno * hex_celldim] = {
+ -0.57735027, -0.57735027, -0.57735027,
+ 0.57735027, -0.57735027, -0.57735027,
+ 0.57735027, 0.57735027, -0.57735027,
+ -0.57735027, 0.57735027, -0.57735027,
+ -0.57735027, -0.57735027, 0.57735027,
+ 0.57735027, -0.57735027, 0.57735027,
+ 0.57735027, 0.57735027, 0.57735027,
+ -0.57735027, 0.57735027, 0.57735027
+ };
+ double hex_qwts[hex_nno] = { 1., 1., 1., 1., 1., 1., 1., 1. };
+
+
+ int nq,nd;
+ double *qx, *qw;
+
+ nq = 0;
+ nd = 0;
+ qx = 0;
+ qw = 0;
+
+
+ if ((points_path != "") && (weights_path != ""))
+ {
+ Reader *points_reader;
+ string points_loc, points_file, points_ext;
+ parse_dataset_path(points_path, points_loc, points_file, points_ext);
+ new_reader(&points_reader, points_ext);
+ ierr = points_reader->open(points_file);
+ assert(ierr == 0);
+ cout << "quadrature points file = " << points_file << endl;
+ points_reader->get_dataset(points_loc.c_str(), &qx, &nq, &nd);
+ assert(nq > 0);
+ assert(nd > 0);
+ //points_reader->close();
+
+ int wtdims[2];
+ Reader *weights_reader;
+ string weights_loc, weights_file, weights_ext;
+ parse_dataset_path(weights_path, weights_loc, weights_file, weights_ext);
+ new_reader(&weights_reader, weights_ext);
+ ierr = weights_reader->open(weights_file);
+ assert(ierr == 0);
+ cout << "quadrature weights file = " << weights_file << endl;
+ weights_reader->get_dataset(weights_loc.c_str(), &qw, &wtdims[0], &wtdims[1]);
+ assert(wtdims[0] == nq);
+ assert(wtdims[1] == 1);
+ //weights_reader->close();
+ }
+ else if (quadrature_path != "")
+ {
+ int ierr;
+ string quadrature_loc, quadrature_file, quadrature_ext;
+ parse_dataset_path(quadrature_path, quadrature_loc, quadrature_file, quadrature_ext);
+ new_reader(&reader, quadrature_ext);
+ ierr = reader->open(quadrature_file);
+ assert(ierr == 0);
+ cout << "quadrature file = " << quadrature_file << endl;
+
+ if (reader->getType() == Reader::HDF_READER)
+ {
+ string points_loc = quadrature_loc + "/points";
+ cout << "quadrature points = " << points_loc << endl;
+ reader->get_dataset(points_loc.c_str(), &qx, &nq, &nd);
+
+ int wtdims[2];
+ string weights_loc = quadrature_loc + "/weights";
+ cout << "quadrature weights = " << weights_loc << endl;
+ reader->get_dataset(weights_loc.c_str(), &qw, &wtdims[0], &wtdims[1]);
+ assert(wtdims[0] == nq);
+ assert(wtdims[1] == 1);
+
+ }
+ else if (reader->getType() == Reader::TXT_READER)
+ {
+ int rows, cols;
+ double *wxdata;
+ reader->get_dataset(0, &wxdata, &rows, &cols);
+ nq = rows;
+ nd = cols - 1;
+ assert(nq >= 1);
+ assert(nd >= 1);
+ qx = new double[nq*nd];
+ qw = new double[nq];
+ for (i = 0; i < nq; i++)
+ {
+ qw[i] = wxdata[cols*i + 0];
+ for (j = 0; j < nd; j++)
+ {
+ qx[nd*i + j] = wxdata[cols*i + (j+1)];
+ }
+ }
+ delete [] wxdata;
+ }
+ else if (reader->getType() == Reader::VTK_READER)
+ {
+ assert(false); // XXX: VtkReader not supported!
+ }
+ }
+ else if (quadrature_order != "")
+ {
+ int order;
+ string_to_int(quadrature_order, order);
+
+ if (order > 0)
+ {
+ /* call FiatProxy
+ //fiat->set(quadrature);
+ int npts,dim;
+ double *qx,*qw;
+ fiat->quadrature(geometry, order, &qx, &qw, &npts, &dim);
+ quadrature->set_quadrature(qx, qw, nno, nsd);
+ delete [] qx;
+ delete [] qw;
+ // */
+ }
+ }
+ else
+ {
+ // assign reasonable defaults
+ switch (cell->geometry())
+ {
+ case Cell::TRIANGLE:
+ quadrature = new QuadraturePoints();
+ quadrature->set_quadrature(tri_qpts, tri_qwts, tri_nno, tri_celldim);
+ quadrature->set_globaldim(tri_nsd);
+ break;
+ case Cell::QUADRANGLE:
+ quadrature = new QuadraturePoints();
+ quadrature->set_quadrature(quad_qpts, quad_qwts, quad_nno, quad_celldim);
+ quadrature->set_globaldim(quad_nsd);
+ break;
+ case Cell::TETRAHEDRON:
+ quadrature = new QuadraturePoints();
+ quadrature->set_quadrature(tet_qpts, tet_qwts, tet_nno, tet_celldim);
+ quadrature->set_globaldim(tet_nsd);
+ break;
+ case Cell::HEXAHEDRON:
+ quadrature = new QuadraturePoints();
+ quadrature->set_quadrature(hex_qpts, hex_qwts, hex_nno, hex_celldim);
+ quadrature->set_globaldim(hex_nsd);
+ break;
+ default:
+ break;
+ }
+ }
+
+ if ((qx != 0) && (qw != 0))
+ {
+ quadrature = new QuadraturePoints();
+ quadrature->set_quadrature(qx, qw, nq, nd);
+ quadrature->set_globaldim(3); // XXX: how to treat 2D case?
+ }
+
+ assert(quadrature != 0);
+ assert(quadrature->n_points() > 0);
+ assert(quadrature->n_refdim() > 0);
+ assert(quadrature->n_globaldim() > 0);
+
+ cout << "quadrature rule = "
+ << quadrature->n_points() << " points, "
+ << "celldim " << quadrature->n_refdim() << ", "
+ << "dim " << quadrature->n_globaldim()
+ << endl;
+
+ //
+ // XXX: now we need to ensure that every point in our
+ // quadrature rule is indeed contained inside our cell.
+ //
+ bool all_contained = true;
+ for (i = 0; i < quadrature->n_points(); i++)
+ {
+ double refpt[3] = {0,0,0};
+ for (j = 0; j < quadrature->n_refdim(); j++)
+ {
+ refpt[j] = quadrature->point(i,j);
+ }
+ if (!cell->interior(refpt[0], refpt[1], refpt[2]))
+ {
+ all_contained = false;
+ break;
+ }
+
+ }
+ if (!all_contained)
+ {
+ // XXX: add more information here (e.g., actual reference cell used,
+ // such as tri3, quad4, tet4, hex8, ...)
+ cerr << "QuadratureIO::load() error: "
+ << "Not all points provided lie inside required reference cell"
+ << endl;
+ exit(1);
+ }
+
+ //
+ // XXX: emit warning if there are any negative quadrature weights?
+ //
+ bool some_neg_wts = false;
+ for (i = 0; i < quadrature->n_points(); i++)
+ {
+ if (quadrature->weight(i) < 0)
+ {
+ some_neg_wts = true;
+ break;
+ }
+ }
+ if (some_neg_wts)
+ {
+ //
+ // XXX: only emit warning. exit(1) is not necessary...
+ //
+ cerr << "QuadratureIO::load() warning: "
+ << "Some weights are negative"
+ << endl;
+ }
+
+ //
+ // Clean up any local allocations
+ //
+ if (qx != 0) delete [] qx;
+ if (qw != 0) delete [] qw;
+}
+
+// ---------------------------------------------------------------------------
Added: cs/benchmark/cigma/trunk/src/hold/QuadratureIO.h
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/QuadratureIO.h (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/QuadratureIO.h 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,36 @@
+#ifndef __QUADRATURE_IO_H__
+#define __QUADRATURE_IO_H__
+
+#include <string>
+#include "AnyOption.h"
+#include "Reader.h"
+#include "Writer.h"
+#include "QuadraturePoints.h"
+#include "Cell.h"
+
+class QuadratureIO
+{
+public:
+ cigma::Reader *reader;
+ cigma::Writer *writer;
+
+public:
+ std::string quadrature_order;
+ std::string quadrature_path;
+ std::string points_path;
+ std::string weights_path;
+ cigma::QuadraturePoints *quadrature;
+
+public:
+ QuadratureIO();
+ ~QuadratureIO();
+ void load(cigma::Cell *cell);
+ void save();
+};
+
+
+void load_args(AnyOption *opt, QuadratureIO *quadratureIO, const char *opt_prefix);
+void validate_args(QuadratureIO *quadratureIO, const char *cmd_name);
+
+
+#endif
Added: cs/benchmark/cigma/trunk/src/hold/Reader.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/Reader.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/Reader.cpp 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,15 @@
+#include "Reader.h"
+#include <cassert>
+#include <cstdlib>
+
+// ---------------------------------------------------------------------------
+cigma::Reader::Reader()
+{
+}
+
+cigma::Reader::~Reader()
+{
+}
+
+// ---------------------------------------------------------------------------
+
Added: cs/benchmark/cigma/trunk/src/hold/Reader.h
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/Reader.h (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/Reader.h 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,36 @@
+#ifndef __READER_H__
+#define __READER_H__
+
+#include <string>
+
+namespace cigma
+{
+ class Reader;
+}
+
+class cigma::Reader
+{
+public:
+ typedef enum {
+ NULL_READER,
+ HDF_READER,
+ VTK_READER,
+ TXT_READER
+ } ReaderType;
+
+public:
+ Reader();
+ ~Reader();
+
+public:
+ virtual ReaderType getType() = 0;
+ virtual void open(std::string filename) = 0;
+ virtual void close() = 0;
+
+public:
+ virtual void get_dataset(const char *loc, double **data, int *num, int *dim) = 0;
+ virtual void get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd) = 0;
+ virtual void get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs) = 0;
+};
+
+#endif
Added: cs/benchmark/cigma/trunk/src/hold/TextReader.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/TextReader.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/TextReader.cpp 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,154 @@
+#include "TextReader.h"
+#include <cassert>
+#include <cstdlib>
+#include <cstring>
+
+// ---------------------------------------------------------------------------
+cigma::TextReader::TextReader()
+{
+ fp = NULL;
+}
+
+cigma::TextReader::~TextReader()
+{
+ close();
+}
+
+// ---------------------------------------------------------------------------
+int cigma::TextReader::open(std::string filename)
+{
+ fp = NULL;
+
+ if (filename != "")
+ {
+ fp = fopen(filename.c_str(), "r");
+ }
+
+ // check for failure
+ if (fp == NULL)
+ {
+ return -1;
+ }
+
+ return 0;
+}
+
+void cigma::TextReader::close()
+{
+ if (fp != NULL)
+ {
+ fclose(fp);
+ }
+}
+
+
+// ---------------------------------------------------------------------------
+
+static FILE *get_fp(const char *loc, FILE *default_fp)
+{
+ FILE *fp = default_fp;
+ if ((loc != NULL) && (strcmp(loc, "") != 0))
+ {
+ fp = fopen(loc, "r");
+ }
+ return fp;
+}
+
+static bool read_dmat(FILE *fp, double **mat, int *rows, int *cols)
+{
+ assert(fp != NULL);
+
+ int nrows, ncols;
+ int ret;
+
+ ret = fscanf(fp, "%d", &nrows);
+ assert(ret == 1);
+ assert(nrows > 0);
+
+ ret = fscanf(fp, "%d", &ncols);
+ assert(ret == 1);
+ assert(ncols > 0);
+
+ double *m = new double[nrows * ncols];
+
+ for (int i = 0; i < nrows; i++)
+ {
+ for (int j = 0; j < ncols; j++)
+ {
+ int n = ncols*i + j;
+ ret = fscanf(fp, "%lf", &m[n]);
+ assert(ret == 1);
+ }
+ }
+
+ *mat = m;
+ *rows = nrows;
+ *cols = ncols;
+
+ return true;
+}
+
+static bool read_imat(FILE *fp, int **mat, int *rows, int *cols)
+{
+ assert(fp != NULL);
+
+ int nrows, ncols;
+ int ret;
+
+ ret = fscanf(fp, "%d", &nrows);
+ assert(ret == 1);
+ assert(nrows > 0);
+
+ ret = fscanf(fp, "%d", &ncols);
+ assert(ret == 1);
+ assert(ncols > 0);
+
+ int *m = new int[nrows * ncols];
+
+ for (int i = 0; i < nrows; i++)
+ {
+ for (int j = 0; j < ncols; j++)
+ {
+ int n = ncols*i + j;
+ ret = fscanf(fp, "%d", &m[n]);
+ assert(ret == 1);
+ }
+ }
+
+ *mat = m;
+ *rows = nrows;
+ *cols = ncols;
+
+ return true;
+}
+
+
+// ---------------------------------------------------------------------------
+
+void cigma::TextReader::get_dataset(const char *loc, double **data, int *num, int *dim)
+{
+ FILE *loc_fp = get_fp(loc, fp);
+ assert(loc_fp != NULL);
+ read_dmat(loc_fp, data, num, dim);
+}
+
+void cigma::TextReader::get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs)
+{
+ // XXX: add support for sections (c.f. gmsh format) so we can scan a single file
+ // for now, interpret loc argument to be an entirely new file
+ FILE *loc_fp = get_fp(loc, fp);
+ assert(loc_fp != NULL);
+ read_imat(loc_fp, connectivity, nel, ndofs);
+}
+
+void cigma::TextReader::get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd)
+{
+ // XXX: add support for sections (c.f. gmsh format) so we can scan a single file
+ // for now, interpret loc argument to be an entirely new file
+ FILE *loc_fp = get_fp(loc, fp);
+ assert(loc_fp != NULL);
+ read_dmat(loc_fp, coordinates, nno, nsd);
+}
+
+
+// ---------------------------------------------------------------------------
Added: cs/benchmark/cigma/trunk/src/hold/TextReader.h
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/TextReader.h (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/TextReader.h 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,38 @@
+#ifndef __TEXTREADER_H__
+#define __TEXTREADER_H__
+
+#include <cstdio>
+#include <string>
+
+#include "Reader.h"
+
+namespace cigma
+{
+ class TextReader;
+};
+
+class cigma::TextReader : public Reader
+{
+public:
+ TextReader();
+ ~TextReader();
+
+public:
+ ReaderType getType() { return TXT_READER; }
+
+public:
+ int open(std::string filename);
+ void close();
+
+public:
+ void get_dataset(const char *loc, double **data, int *num, int *dim);
+ void get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd);
+ void get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs);
+
+public:
+ FILE *fp; // default file pointer
+};
+
+// ---------------------------------------------------------------------------
+
+#endif
Added: cs/benchmark/cigma/trunk/src/hold/TextWriter.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/TextWriter.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/TextWriter.cpp 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,113 @@
+#include "TextWriter.h"
+#include <cassert>
+#include <cstdlib>
+
+// ---------------------------------------------------------------------------
+cigma::TextWriter::TextWriter()
+{
+ fp = NULL;
+}
+
+cigma::TextWriter::~TextWriter()
+{
+ close();
+}
+
+// ---------------------------------------------------------------------------
+
+int cigma::TextWriter::open(std::string filename)
+{
+ fp = NULL;
+
+ if (filename != "")
+ {
+ fp = fopen(filename.c_str(), "w");
+ }
+ else
+ {
+ fp = stdout;
+ }
+
+ if (fp == NULL)
+ {
+ return -1;
+ }
+ return 0;
+}
+
+void cigma::TextWriter::close()
+{
+ if (fp != NULL)
+ {
+ fclose(fp);
+ }
+}
+
+
+// ---------------------------------------------------------------------------
+
+void cigma::TextWriter::
+write_field(FE_Field *field)
+{
+ assert(field != 0);
+}
+
+
+// ---------------------------------------------------------------------------
+
+static bool write_dmat(FILE *fp, double *mat, int rows, int cols)
+{
+ assert(fp != NULL);
+ assert(rows > 0);
+ assert(cols > 0);
+
+ fprintf(fp, "%d %d\n", rows, cols);
+ for (int i = 0; i < rows; i++)
+ {
+ for (int j = 0; j < cols; j++)
+ {
+ fprintf(fp, " %g", mat[cols*i + j]);
+ }
+ fprintf(fp, "\n");
+ }
+
+ return true;
+}
+
+static bool write_imat(FILE *fp, int *mat, int rows, int cols)
+{
+ assert(fp != NULL);
+ assert(rows > 0);
+ assert(cols > 0);
+
+ fprintf(fp, "%d %d\n", rows, cols);
+ for (int i = 0; i < rows; i++)
+ {
+ for (int j = 0; j < cols; j++)
+ {
+ fprintf(fp, " %d", mat[cols*i + j]);
+ }
+ fprintf(fp, "\n");
+ }
+
+ return true;
+}
+
+// ---------------------------------------------------------------------------
+
+void cigma::TextWriter::write_connectivity(int *connectivity, int nel, int ndofs)
+{
+ write_imat(fp, connectivity, nel, ndofs);
+}
+
+void cigma::TextWriter::write_coordinates(double *coordinates, int nno, int nsd)
+{
+ write_dmat(fp, coordinates, nno, nsd);
+}
+
+void cigma::TextWriter::write_dofs(double *dofs, int num, int ndim)
+{
+ write_dmat(fp, dofs, num, ndim);
+}
+
+// ---------------------------------------------------------------------------
Added: cs/benchmark/cigma/trunk/src/hold/TextWriter.h
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/TextWriter.h (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/TextWriter.h 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,42 @@
+#ifndef __TEXTWRITER_H__
+#define __TEXTWRITER_H__
+
+#include <cstdio>
+#include <string>
+
+#include "Writer.h"
+
+
+namespace cigma
+{
+ class TextWriter;
+};
+
+
+class cigma::TextWriter : public Writer
+{
+public:
+ TextWriter();
+ ~TextWriter();
+
+public:
+ WriterType getType() { return TXT_WRITER; }
+ int open(std::string filename);
+ void close();
+
+public:
+ void write_field(FE_Field *field);
+
+public:
+ void write_connectivity(int *connectivity, int nel, int ndofs);
+ void write_coordinates(double *coordinates, int nno, int nsd);
+ void write_dofs(double *dofs, int nno, int ndim);
+
+public:
+ FILE *fp;
+};
+
+
+// ---------------------------------------------------------------------------
+
+#endif
Added: cs/benchmark/cigma/trunk/src/hold/VtkList.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/VtkList.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/VtkList.cpp 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,250 @@
+#include <cassert>
+#include "VtkList.h"
+
+#include "vtkDataSetReader.h"
+#include "vtkXMLStructuredGridReader.h"
+
+#include "vtkDataSet.h"
+#include "vtkStructuredGrid.h"
+
+#include "vtkPointData.h"
+#include "vtkCellData.h"
+#include "vtkDataArray.h"
+
+
+using namespace std;
+
+
+// ---------------------------------------------------------------------------
+
+/*
+VtkFileType getFileType(vtkDataSetReader *reader)
+{
+ // XXX: are these mutually exclusive?
+ if (reader->IsFilePolyData())
+ return VTK_FILE_POLYDATA;
+ if (reader->IsFileStructuredPoints())
+ return VTK_FILE_STRUCTURED_POINTS;
+ if (reader->IsFileStructuredGrid())
+ return VTK_FILE_STRUCTURED_GRID;
+ if (reader->IsFileUnstructuredGrid())
+ return VTK_FILE_UNSTRUCTURED_GRID;
+ if (reader->IsFileRectilinearGrid())
+ return VTK_FILE_RECTILINEAR_GRID;
+ return VTK_FILE_NONE;
+
+} // */
+
+
+const char *getFileTypeName(VtkFileType fileType)
+{
+ if (fileType == VTK_FILE_POLYDATA)
+ return "polydata";
+ if (fileType == VTK_FILE_STRUCTURED_POINTS)
+ return "structured points";
+ if (fileType == VTK_FILE_STRUCTURED_GRID)
+ return "structured grid";
+ if (fileType == VTK_FILE_UNSTRUCTURED_GRID)
+ return "unstructured grid";
+ if (fileType == VTK_FILE_RECTILINEAR_GRID)
+ return "rectilinear grid";
+ return "none";
+}
+
+
+// ---------------------------------------------------------------------------
+
+void vtkls(vtkDataSet *dataset)
+{
+ assert(dataset != 0);
+
+ vtkPointData *pointData = dataset->GetPointData();
+ //pointData->PrintSelf(cout, 0);
+
+ vtkCellData *cellData = dataset->GetCellData();
+ //cellData->PrintSelf(cout, 0);
+
+
+ // ---------------------------------------------------
+
+ int nno, nel;
+
+ nno = dataset->GetNumberOfPoints();
+ nel = dataset->GetNumberOfCells();
+
+ if (nno > 0)
+ {
+ cout << "Points = " << nno << endl;
+ }
+
+ if (nel > 0)
+ {
+ cout << "Cells = " << nel << endl;
+ }
+
+
+ // ---------------------------------------------------
+ int i, n;
+ int numArrays;
+ int numTuples, numComponents;
+
+ numArrays = pointData->GetNumberOfArrays();
+ if (numArrays > 0)
+ {
+ for (i = 0; i < numArrays; i++)
+ {
+ vtkDataArray *dataArray = pointData->GetArray(i);
+ const char *name = pointData->GetArrayName(i);
+ numTuples = dataArray->GetNumberOfTuples();
+ numComponents = dataArray->GetNumberOfComponents();
+ cout << "PointDataArray[" << i << "] = " << name << " ";
+ //cout << dataArray->GetClassName() << " ";
+ cout << "(" << numTuples << " x " << numComponents << ")";
+ cout << endl;
+ }
+ }
+
+ numArrays = cellData->GetNumberOfArrays();
+ if (numArrays > 0)
+ {
+ for (i = 0; i < numArrays; i++)
+ {
+ vtkDataArray *dataArray = cellData->GetArray(i);
+ const char *name = cellData->GetArrayName(i);
+ numTuples = dataArray->GetNumberOfTuples();
+ numComponents = dataArray->GetNumberOfComponents();
+ cout << "CellDataArray[" << i << "] = " << name << " ";
+ //cout << dataArray->GetClassName() << " ";
+ cout << "(" << numTuples << " x " << numComponents << ")";
+ cout << endl;
+ }
+ }
+
+ return;
+}
+
+
+// ---------------------------------------------------------------------------
+
+void list_vtk(const char *filename)
+{
+ int outputType;
+ vtkDataSetReader *reader = vtkDataSetReader::New();
+
+ reader->SetFileName(filename);
+ cout << "Reading " << filename << endl;
+
+ int ierr = -1;
+
+ ierr = reader->OpenVTKFile(); // does file exist?
+ if (ierr == 0)
+ {
+ cerr << "Could not open " << filename << endl;
+ exit(1);
+ }
+
+ ierr = reader->ReadHeader(); // is the vtk header present?
+ if (ierr == 0)
+ {
+ cerr << "Unrecognized file " << filename << endl;
+ exit(1);
+ }
+
+ outputType = reader->ReadOutputType();
+ if (outputType < 0)
+ {
+ cerr << "Invalid VTK file? " << filename << endl;
+ exit(1);
+ }
+
+
+ // ---------------------------------------------------
+ reader->Update();
+ //reader->PrintSelf(cout, 0);
+
+ vtkDataSet *dataset = reader->GetOutput();
+ //dataset->PrintSelf(cout, 0);
+
+ const bool group_by_arrays = true;
+ const bool group_by_class = false;
+
+ if (group_by_arrays)
+ {
+ vtkls(dataset);
+ }
+
+ if (group_by_class)
+ {
+ int n, i;
+
+ n = reader->GetNumberOfFieldDataInFile();
+ for (i = 0; i < n; i++)
+ {
+ const char *name = reader->GetFieldDataNameInFile(i);
+ //cout << "Number of field data = " << n << endl;
+ cout << "FieldData[" << i << "] = " << name;
+ cout << endl;
+ }
+
+ n = reader->GetNumberOfScalarsInFile();
+ //cout << "Number of scalars = " << n << endl;
+ for (i = 0; i < n; i++)
+ {
+ const char *name = reader->GetScalarsNameInFile(i);
+ cout << "Scalars[" << i << "] = " << name;
+ cout << endl;
+ }
+
+ n = reader->GetNumberOfVectorsInFile();
+ //cout << "Number of vectors = " << n << endl;
+ for (i = 0; i < n; i++)
+ {
+ const char *name = reader->GetVectorsNameInFile(i);
+ cout << "Vectors[" << i << "] = " << name;
+ cout << endl;
+ }
+
+ n = reader->GetNumberOfTensorsInFile();
+ //cout << "Number of tensors = " << n << endl;
+ for (i = 0; i < n; i++)
+ {
+ const char *name = reader->GetTensorsNameInFile(i);
+ cout << "Tensors[" << i << "] = " << name;
+ cout << endl;
+ }
+ }
+
+ reader->Delete();
+
+ return;
+}
+
+
+void list_vts(const char *filename)
+{
+ vtkXMLStructuredGridReader *reader = vtkXMLStructuredGridReader::New();
+
+ int canReadFile = 0;
+ canReadFile = reader->CanReadFile(filename);
+
+ if (!canReadFile)
+ {
+ cerr << "Could not read " << filename << endl;
+ exit(1);
+ }
+
+ reader->SetFileName(filename);
+ cout << "Reading " << filename << endl;
+
+ reader->Update();
+ //reader->PrintSelf(cout, 0);
+
+ vtkStructuredGrid *sgrid = reader->GetOutput();
+ vtkDataSet *dataset = static_cast<vtkDataSet*>(sgrid);
+
+ vtkls(dataset);
+
+ reader->Delete();
+}
+
+// ---------------------------------------------------------------------------
Added: cs/benchmark/cigma/trunk/src/hold/VtkList.h
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/VtkList.h (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/VtkList.h 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,21 @@
+#ifndef __VTK_LIST_H__
+#define __VTK_LIST_H__
+
+
+typedef enum {
+ VTK_FILE_NONE,
+ VTK_FILE_POLYDATA,
+ VTK_FILE_STRUCTURED_POINTS,
+ VTK_FILE_STRUCTURED_GRID,
+ VTK_FILE_UNSTRUCTURED_GRID,
+ VTK_FILE_RECTILINEAR_GRID
+} VtkFileType;
+
+//VtkFileType getFileType(vtkDataSetReader *reader);
+const char *getFileTypeName(VtkFileType fileType);
+
+void list_vtk(const char *filename);
+void list_vts(const char *filename);
+
+
+#endif
Added: cs/benchmark/cigma/trunk/src/hold/VtkReader.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/VtkReader.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/VtkReader.cpp 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,286 @@
+#include <iostream>
+#include <cstdlib>
+#include <cassert>
+
+#include "VtkReader.h"
+
+#include "vtkCellArray.h"
+#include "vtkPointData.h"
+#include "vtkDoubleArray.h"
+#include "vtkFloatArray.h"
+
+
+// ---------------------------------------------------------------------------
+
+cigma::VtkReader::VtkReader()
+{
+ reader = 0;
+ grid = 0;
+}
+
+cigma::VtkReader::~VtkReader()
+{
+ //XXX: avoid calling close() for now
+ //close();
+}
+
+
+
+// ---------------------------------------------------------------------------
+
+int cigma::VtkReader::open(std::string filename)
+{
+ int ierr;
+
+ /* XXX: throw exception if file doesn't exist */
+ reader = vtkUnstructuredGridReader::New();
+ reader->SetFileName(filename.c_str());
+
+
+ //* detect bad vtk file
+ ierr = reader->OpenVTKFile();
+ if (ierr == 0)
+ {
+ //reader->Delete(); // XXX
+ reader = 0;
+ return -1;
+ }
+
+ ierr = reader->ReadHeader();
+ if (ierr == 0)
+ {
+ //reader->Delete(); // XXX
+ reader = 0;
+ return -1;
+ } // */
+
+
+ reader->Update();
+ //reader->PrintSelf(std::cout, 0);
+
+ grid = reader->GetOutput();
+ //grid->PrintSelf(std::cout, 4);
+
+ return 0;
+}
+
+void cigma::VtkReader::close()
+{
+ /* XXX: fix this!
+ if (grid != 0)
+ {
+ grid->Delete();
+ grid = 0;
+ }
+ if (reader != 0)
+ {
+ reader->Delete();
+ reader = 0;
+ }
+ // */
+}
+
+
+
+// ---------------------------------------------------------------------------
+
+void cigma::VtkReader::
+get_dataset(const char *loc, double **data, int *num, int *dim)
+{
+ assert(grid != 0);
+ assert(loc != 0);
+
+ int size = 0;
+ int dims[2] = {0, 0};
+ double *array = 0;
+
+ vtkPointData *pointData = grid->GetPointData();
+
+ if (pointData->HasArray(loc) == 1)
+ {
+ vtkDataArray *dataArray = pointData->GetArray(loc);
+ assert(dataArray != 0);
+
+ dims[0] = dataArray->GetNumberOfTuples();
+ dims[1] = dataArray->GetNumberOfComponents();
+
+ size = dims[0] * dims[1];
+ array = new double[size];
+
+ int dataType = dataArray->GetDataType();
+
+ if (dataType == VTK_DOUBLE)
+ {
+ double *ptr = static_cast<double*>(dataArray->GetVoidPointer(0));
+ for (int i = 0; i < size; i++)
+ {
+ array[i] = ptr[i];
+ }
+ }
+ else if (dataType == VTK_FLOAT)
+ {
+ float *ptr = static_cast<float*>(dataArray->GetVoidPointer(0));
+ for (int i = 0; i < size; i++)
+ {
+ array[i] = ptr[i];
+ }
+ }
+ else
+ {
+ assert(false); // XXX: throw exception instead
+ }
+ }
+
+ *data = array;
+ *num = dims[0];
+ *dim = dims[1];
+}
+
+void cigma::VtkReader::
+get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd)
+{
+ get_coordinates(coordinates, nno, nsd);
+}
+
+void cigma::VtkReader::
+get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs)
+{
+ get_connectivity(connectivity, nel, ndofs);
+}
+
+
+// ---------------------------------------------------------------------------
+
+void cigma::VtkReader::
+get_coordinates(double **coordinates, int *nno, int *nsd)
+{
+ assert(grid != 0);
+
+ vtkPoints *points = grid->GetPoints();
+ //points->PrintSelf(std::cout, 0);
+
+ int dims[2];
+ dims[0] = points->GetNumberOfPoints();
+ dims[1] = 3;
+
+ int size = dims[0] * dims[1];
+ int dataType = points->GetDataType();
+ double *coords = new double[size];
+
+ if (dataType == VTK_DOUBLE)
+ {
+ double *ptr = static_cast<double*>(points->GetVoidPointer(0));
+ for (int i = 0; i < size; i++)
+ {
+ coords[i] = ptr[i];
+ }
+ }
+ else if (dataType == VTK_FLOAT)
+ {
+ float *ptr = static_cast<float*>(points->GetVoidPointer(0));
+ for (int i = 0; i < size; i++)
+ {
+ coords[i] = ptr[i];
+ }
+ }
+ else
+ {
+ assert(false); // XXX: throw exception instead
+ }
+
+ *coordinates = coords;
+ *nno = dims[0];
+ *nsd = dims[1];
+}
+
+void cigma::VtkReader::
+get_connectivity(int **connectivity, int *nel, int *ndofs)
+{
+ assert(grid != 0);
+
+ vtkCellArray *cellArray = grid->GetCells();
+ //cellArray->PrintSelf(std::cout, 0);
+
+ vtkIdType numCells = grid->GetNumberOfCells();
+ vtkIdType *cellArrayPtr = cellArray->GetPointer();
+
+ int dofs_per_elt = cellArrayPtr[0];
+ int offset = dofs_per_elt + 1;
+ int *connect = new int[numCells * dofs_per_elt];
+
+ for (int e = 0; e < numCells; ++e)
+ {
+ for (int i = 1; i <= dofs_per_elt; ++i)
+ {
+ connect[dofs_per_elt * e + (i-1)] = cellArrayPtr[offset*e + i];
+ }
+ }
+
+ *connectivity = connect;
+ *nel = numCells;
+ *ndofs = dofs_per_elt;
+}
+
+void cigma::VtkReader::
+get_vector_point_data(const char *name, double **vectors, int *num, int *dim)
+{
+ assert(grid != 0);
+
+ vtkPointData *pointData = grid->GetPointData();
+ //pointData->PrintSelf(std::cout, 0);
+
+ vtkDataArray *dataArray;
+ if (name != 0) {
+ assert(pointData->HasArray(name) == 1);
+ dataArray = pointData->GetVectors(name);
+ } else {
+ dataArray = pointData->GetVectors();
+ }
+ //dataArray->PrintSelf(std::cout, 0);
+
+ int dataType = dataArray->GetDataType();
+ assert(dataType == VTK_DOUBLE);
+ double *ptr = static_cast<double*>(dataArray->GetVoidPointer(0));
+
+ // XXX: copy the data, or keep the reference?
+ // if dataType from the file is float, then we'd need to copy anyway
+ // perhaps we need a void pointer and a type
+ // new function signature would be
+ // void get_vector_point_data(const char *name, void **vectors, int *num, int *dim, int *type)
+
+ *vectors = ptr;
+ *num = dataArray->GetNumberOfTuples();
+ *dim = dataArray->GetNumberOfComponents();
+}
+
+void cigma::VtkReader::
+get_scalar_point_data(const char *name, double **scalars, int *num, int *dim)
+{
+ assert(grid != 0);
+
+ vtkPointData *pointData = grid->GetPointData();
+ //pointData->PrintSelf(std::cout, 0);
+
+ vtkDataArray *dataArray;
+ if (name != 0) {
+ assert(pointData->HasArray(name) == 1);
+ dataArray = pointData->GetScalars(name);
+ } else {
+ dataArray = pointData->GetScalars();
+ }
+ //dataArray->PrintSelf(std::cout, 0);
+
+
+ int dataType = dataArray->GetDataType();
+ assert(dataType == VTK_DOUBLE);
+ double *ptr = static_cast<double*>(dataArray->GetVoidPointer(0));
+
+ // XXX: see comment in get_vector_point_data()
+
+ *scalars = ptr;
+ *num = dataArray->GetNumberOfTuples();
+ *dim = dataArray->GetNumberOfComponents();
+}
+
+
+// ---------------------------------------------------------------------------
Added: cs/benchmark/cigma/trunk/src/hold/VtkReader.h
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/VtkReader.h (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/VtkReader.h 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,80 @@
+#ifndef __VTK_READER_H__
+#define __VTK_READER_H__
+
+#include <string>
+
+#include "Reader.h"
+
+#include "vtkDataSetReader.h"
+#include "vtkUnstructuredGridReader.h"
+#include "vtkUnstructuredGrid.h"
+
+namespace cigma
+{
+ class VtkReader;
+}
+
+class cigma::VtkReader : public cigma::Reader
+{
+public:
+ typedef enum {
+ VTK_FILE_NONE,
+ VTK_FILE_POLYDATA,
+ VTK_FILE_STRUCTURED_POINTS,
+ VTK_FILE_STRUCTURED_GRID,
+ VTK_FILE_UNSTRUCTURED_GRID,
+ VTK_FILE_RECTILINEAR_GRID
+ } VtkFileType;
+
+ ReaderType getType()
+ {
+ return VTK_READER;
+ }
+
+ VtkFileType getFileType()
+ {
+ if (reader != 0)
+ {
+ // XXX: are these mutually exclusive?
+ if (reader->IsFilePolyData())
+ return VTK_FILE_POLYDATA;
+ if (reader->IsFileStructuredPoints())
+ return VTK_FILE_STRUCTURED_POINTS;
+ if (reader->IsFileStructuredGrid())
+ return VTK_FILE_STRUCTURED_GRID;
+ if (reader->IsFileUnstructuredGrid())
+ return VTK_FILE_UNSTRUCTURED_GRID;
+ if (reader->IsFileRectilinearGrid())
+ return VTK_FILE_RECTILINEAR_GRID;
+ }
+ return VTK_FILE_NONE;
+ }
+
+public:
+ VtkReader();
+ ~VtkReader();
+
+public:
+ int open(std::string filename);
+ void close();
+
+public:
+ void get_dataset(const char *loc, double **data, int *num, int *dim);
+ void get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd);
+ void get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs);
+
+public:
+ void get_coordinates(double **coordinates, int *nno, int *nsd);
+ void get_connectivity(int **connectivity, int *nel, int *ndofs);
+ void get_vector_point_data(const char *name, double **vectors, int *num, int *dim);
+ void get_scalar_point_data(const char *name, double **scalars, int *num, int *dim);
+
+public:
+ vtkUnstructuredGridReader *reader;
+ //vtkDataSetReader *reader;
+ vtkUnstructuredGrid *grid;
+};
+
+// ---------------------------------------------------------------------------
+
+#endif
Added: cs/benchmark/cigma/trunk/src/hold/VtkWriter.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/VtkWriter.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/VtkWriter.cpp 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,221 @@
+#include "VtkWriter.h"
+#include <cassert>
+#include <cstdlib>
+
+// ---------------------------------------------------------------------------
+cigma::VtkWriter::VtkWriter()
+{
+ fp = NULL;
+}
+
+cigma::VtkWriter::~VtkWriter()
+{
+ close();
+}
+
+// ---------------------------------------------------------------------------
+
+int cigma::VtkWriter::
+open(std::string filename)
+{
+ fp = fopen(filename.c_str(), "w");
+ if (fp == NULL)
+ {
+ return -1;
+ }
+ return 0;
+}
+
+void cigma::VtkWriter::
+close()
+{
+ if (fp != NULL)
+ {
+ fclose(fp);
+ fp = NULL;
+ }
+}
+
+
+// ---------------------------------------------------------------------------
+
+void cigma::VtkWriter::
+write_field(FE_Field *field)
+{
+ assert(field != 0);
+}
+
+// ---------------------------------------------------------------------------
+
+void cigma::VtkWriter::
+write_header()
+{
+ assert(fp != NULL);
+ fprintf(fp, "# vtk DataFile Version 3.0\n");
+ fprintf(fp, "This line is a comment\n");
+ fprintf(fp, "ASCII\n");
+ fprintf(fp, "DATASET UNSTRUCTURED_GRID\n");
+}
+
+
+void cigma::VtkWriter::
+write_points(double *points, int npts, int ndim)
+{
+ assert(fp != NULL);
+ assert(ndim > 1);
+
+ fprintf(fp, "POINTS %d double\n", npts);
+ for (int i = 0; i < npts; i++)
+ {
+ fprintf(fp, " %g", points[ndim*i + 0]);
+ fprintf(fp, " %g", points[ndim*i + 1]);
+ if (ndim == 3)
+ fprintf(fp, " %g", points[ndim*i + 2]);
+ else
+ fprintf(fp, " 0.0");
+ fprintf(fp, "\n");
+ }
+}
+
+void cigma::VtkWriter::
+write_cells(int *cells, int nel, int ndofs)
+{
+ assert(fp != NULL);
+
+ fprintf(fp, "CELLS %d %d\n", nel, nel*(1 + ndofs));
+ for (int i = 0; i < nel; i++)
+ {
+ fprintf(fp, " %d", ndofs);
+ for (int j = 0; j < ndofs; j++)
+ {
+ fprintf(fp, " %d", cells[ndofs*i + j]);
+ }
+ fprintf(fp, "\n");
+ }
+}
+
+
+void cigma::VtkWriter::
+write_cell_types(int nsd, int nel, int ndofs)
+{
+ assert(fp != NULL);
+
+ fprintf(fp, "CELL_TYPES %d\n", nel);
+
+ int vtkType = 0;
+ if (nsd == 3)
+ {
+ switch (ndofs)
+ {
+ case 4: vtkType = 10; break; // VTK_TETRA=10
+ case 8: vtkType = 12; break; // VTK_HEXAHEDRON=12
+ case 6: vtkType = 13; break; // VTK_WEDGE=13
+ case 5: vtkType = 14; break; // VTK_PYRAMID=14
+ case 10: vtkType = 24; break; // VTK_QUADRATIC_TETRA=24
+ case 20: vtkType = 25; break; // VTK_QUADRATIC_HEXAHEDRON=25
+ }
+ }
+ else if (nsd == 2)
+ {
+ switch (ndofs)
+ {
+ case 3: vtkType = 5; break; // VTK_TRIANGLE=5
+ case 4: vtkType = 9; break; // VTK_QUAD=9
+ case 6: vtkType = 22; break; // VTK_QUADRATIC_TRIANGLE=22
+ case 8: vtkType = 23; break; // VTK_QUADRATIC_QUAD=23
+ }
+ }
+ assert(vtkType > 0);
+
+ for (int i = 0; i < nel; i++)
+ {
+ fprintf(fp, "%d\n", vtkType);
+ }
+}
+
+void cigma::VtkWriter::
+write_point_data(const char *name, double *data, int nno, int ndim)
+{
+ assert(fp != NULL);
+
+ fprintf(fp, "POINT_DATA %d\n", nno);
+
+ if (ndim == 1)
+ {
+ fprintf(fp, "SCALARS %s double 1\n", name);
+ fprintf(fp, "LOOKUP_TABLE default\n");
+ for (int i = 0; i < nno; i++)
+ {
+ fprintf(fp, "%g\n", data[i]);
+ }
+ return;
+ }
+ else if ((ndim == 2) || (ndim == 3))
+ {
+ fprintf(fp, "VECTORS %s double\n", name);
+ for (int i = 0; i < nno; i++)
+ {
+ fprintf(fp, " %g", data[ndim*i + 0]);
+ fprintf(fp, " %g", data[ndim*i + 1]);
+ if (ndim == 3)
+ fprintf(fp, " %g\n", data[ndim*i + 2]);
+ else
+ fprintf(fp, " 0.0\n");
+ }
+ }
+ else if (ndim == 9)
+ {
+ fprintf(fp, "TENSORS %s double\n", name);
+ for (int i = 0; i < nno; i++)
+ {
+ for (int j = 0; j < ndim; j++)
+ {
+ fprintf(fp, " %g", data[ndim*i + j]);
+ }
+ fprintf(fp, "\n");
+ }
+ }
+}
+
+void cigma::VtkWriter::
+write_cell_data(const char *name, double *data, int nel, int ndim)
+{
+ assert(fp != NULL);
+
+ fprintf(fp, "CELL_DATA %d\n", nel);
+
+ if (ndim == 1)
+ {
+ fprintf(fp, "SCALARS %s float 1\n", name);
+ fprintf(fp, "LOOKUP_TABLE default\n");
+ for (int i = 0; i < nel; i++)
+ fprintf(fp, "%g\n", data[i]);
+ }
+ else if ((ndim == 2) || (ndim == 3))
+ {
+ fprintf(fp, "VECTORS %s float\n", name);
+ for (int i = 0; i < nel; i++)
+ {
+ fprintf(fp, " %g", data[ndim*i + 0]);
+ fprintf(fp, " %g", data[ndim*i + 1]);
+ if (ndim == 3)
+ fprintf(fp, " %g\n", data[ndim*i + 2]);
+ else
+ fprintf(fp, " 0.0\n");
+ }
+ }
+ else if (ndim == 9)
+ {
+ fprintf(fp, "TENSORS %s float\n", name);
+ for (int i = 0; i < nel; i++)
+ {
+ for (int j = 0; j < ndim; j++)
+ {
+ fprintf(fp, " %g", data[ndim*i + j]);
+ }
+ fprintf(fp, "\n");
+ }
+ }
+}
+
+// ---------------------------------------------------------------------------
Added: cs/benchmark/cigma/trunk/src/hold/VtkWriter.h
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/VtkWriter.h (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/VtkWriter.h 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,41 @@
+#ifndef __VTK_WRITER_H__
+#define __VTK_WRITER_H__
+
+#include <cstdio>
+
+#include "Writer.h"
+
+
+namespace cigma
+{
+ class VtkWriter;
+}
+
+
+class cigma::VtkWriter : public Writer
+{
+public:
+ VtkWriter();
+ ~VtkWriter();
+
+public:
+ WriterType getType() { return VTK_WRITER; }
+ int open(std::string filename);
+ void close();
+
+public:
+ void write_field(FE_Field *field);
+
+public:
+ void write_header();
+ void write_points(double *points, int npts, int ndim);
+ void write_cells(int *cells, int nel, int ndofs);
+ void write_cell_types(int nsd, int nel, int ndofs);
+ void write_point_data(const char *name, double *data, int nno, int ndim);
+ void write_cell_data(const char *name, double *data, int nel, int ndim);
+
+public:
+ FILE *fp;
+};
+
+#endif
Added: cs/benchmark/cigma/trunk/src/hold/VtkXmlReader.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/VtkXmlReader.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/VtkXmlReader.cpp 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,223 @@
+#include <iostream>
+#include <cassert>
+#include "VtkXmlReader.h"
+
+#include "vtkCellArray.h"
+#include "vtkPointData.h"
+#include "vtkDoubleArray.h"
+#include "vtkFloatArray.h"
+
+
+using namespace std;
+
+
+// ---------------------------------------------------------------------------
+
+cigma::VtkXmlReader::VtkXmlReader()
+{
+ reader = 0;
+ grid = 0;
+}
+
+cigma::VtkXmlReader::~VtkXmlReader()
+{
+ //close(); // XXX: make sure this call resets both reader and grid
+}
+
+// ---------------------------------------------------------------------------
+
+int cigma::VtkXmlReader::
+open(std::string filename)
+{
+ reader = vtkXMLStructuredGridReader::New();
+ reader->SetFileName(filename.c_str());
+ reader->Update();
+ //cout << "Reading " << filename << endl;
+ //reader->PrintSelf(cout, 4)
+
+ grid = reader->GetOutput();
+ //cout << endl << endl;
+ //grid->PrintSelf(cout, 4);
+
+ return 0;
+}
+
+void cigma::VtkXmlReader::
+close()
+{
+ // XXX: don't forget to reset values of reader and grid
+}
+
+
+// ---------------------------------------------------------------------------
+
+void cigma::VtkXmlReader::
+get_dataset(const char *loc, double **data, int *num, int *dim)
+{
+ assert(grid != 0);
+ assert(loc != 0);
+
+ int size = 0;
+ int dims[2] = {0, 0};
+ double *array = 0;
+
+ vtkPointData *pointData = grid->GetPointData();
+
+ if (pointData->HasArray(loc) == 1)
+ {
+ vtkDataArray *dataArray = pointData->GetArray(loc);
+ assert(dataArray != 0);
+
+ dims[0] = dataArray->GetNumberOfTuples();
+ dims[1] = dataArray->GetNumberOfComponents();
+
+ size = dims[0] * dims[1];
+ array = new double[size];
+
+ int dataType = dataArray->GetDataType();
+
+ if (dataType == VTK_DOUBLE)
+ {
+ double *ptr = static_cast<double*>(dataArray->GetVoidPointer(0));
+ for (int i = 0; i < size; i++)
+ {
+ array[i] = ptr[i];
+ }
+ }
+ else if (dataType == VTK_FLOAT)
+ {
+ float *ptr = static_cast<float*>(dataArray->GetVoidPointer(0));
+ for (int i = 0; i < size; i++)
+ {
+ array[i] = ptr[i];
+ }
+ }
+ else
+ {
+ assert(false); // XXX: throw exception instead
+ }
+ }
+
+ *data = array;
+ *num = dims[0];
+ *dim = dims[1];
+}
+
+void cigma::VtkXmlReader::
+get_coordinates(const char *loc, double **coordinates, int *nno, int *nsd)
+{
+ get_coordinates(coordinates, nno, nsd);
+}
+
+void cigma::VtkXmlReader::
+get_connectivity(const char *loc, int **connectivity, int *nel, int *ndofs)
+{
+ get_connectivity(connectivity, nel, ndofs);
+}
+
+
+// ---------------------------------------------------------------------------
+
+void cigma::VtkXmlReader::
+get_coordinates(double **coordinates, int *nno, int *nsd)
+{
+ assert(grid != 0);
+
+ vtkPoints *points = grid->GetPoints();
+ //points->PrintSelf(cout, 0);
+
+
+ // XXX: need to know if we're loading 2D cells, so we can
+ // recalculate size & nsd....if 3D, we leave things alone
+
+
+ int dims[2];
+ dims[0] = points->GetNumberOfPoints();
+ //dims[1] = 3;
+ dims[1] = 2; // XXX: for now, assume 2D grid of quadrangles
+
+ int size = dims[0] * dims[1];
+ double *coords = new double[size];
+
+
+ int dataType = points->GetDataType();
+
+ if (dataType == VTK_DOUBLE)
+ {
+ double *ptr = static_cast<double*>(points->GetVoidPointer(0));
+ for (int i = 0; i < dims[0]; i++)
+ {
+ coords[2*i + 0] = ptr[3*i + 0];
+ coords[2*i + 1] = ptr[3*i + 1];
+ } /*
+ for (int i = 0; i < size; i++)
+ {
+ coords[i] = ptr[i];
+ } // */
+ }
+ else if (dataType == VTK_FLOAT)
+ {
+ float *ptr = static_cast<float*>(points->GetVoidPointer(0));
+ for (int i = 0; i < dims[0]; i++)
+ {
+ coords[2*i + 0] = ptr[3*i + 0];
+ coords[2*i + 1] = ptr[3*i + 1];
+ } /*
+ for (int i = 0; i < size; i++)
+ {
+ coords[i] = ptr[i];
+ } // */
+ }
+ else
+ {
+ assert(false); // XXX: throw exception instead
+ }
+
+ *coordinates = coords;
+ *nno = dims[0];
+ *nsd = dims[1];
+}
+
+void cigma::VtkXmlReader::
+get_connectivity(int **connectivity, int *nel, int *ndofs)
+{
+ //
+ // XXX: for now, assume 2D grid of quadrangles
+ //
+ int nx, ny;
+ int ex, ey;
+ int *connect = 0;
+ int ndofs_per_cell;
+ int num_cells;
+ int extent[6];
+
+ num_cells = grid->GetNumberOfCells();
+ ndofs_per_cell = 4;
+
+ grid->GetExtent(extent);
+
+ nx = extent[1] + 1;
+ ny = extent[3] + 1;
+
+ ex = nx - 1;
+ ey = ny - 1;
+
+ connect = new int[ex * ey * ndofs_per_cell];
+ for (int j = 0; j < ny-1; j++)
+ {
+ for (int i = 0; i < nx-1; i++)
+ {
+ int e = i + j*ex;
+ connect[4*e + 0] = i + j*nx;
+ connect[4*e + 1] = (i+1) + j*nx;
+ connect[4*e + 2] = (i+1) + (j+1)*nx;
+ connect[4*e + 3] = i + (j+1)*nx;
+ }
+ }
+
+ *connectivity = connect;
+ *nel = num_cells;
+ *ndofs = ndofs_per_cell;
+}
+
+// ---------------------------------------------------------------------------
Added: cs/benchmark/cigma/trunk/src/hold/VtkXmlReader.h
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/VtkXmlReader.h (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/VtkXmlReader.h 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,64 @@
+#ifndef __VTK_XML_READER_H__
+#define __VTK_XML_READER_H__
+
+#include <string>
+#include "Reader.h"
+
+#include "vtkXMLStructuredGridReader.h"
+#include "vtkStructuredGrid.h"
+
+namespace cigma
+{
+ class VtkXmlReader;
+}
+
+class cigma::VtkXmlReader : public Reader
+{
+public:
+ typedef enum {
+ VTK_FILE_NONE,
+ VTK_FILE_STRUCTURED_GRID
+ // XXX: others...
+ } VtkXmlFileType;
+
+ ReaderType getType()
+ {
+ return VTK_READER;
+ }
+
+ VtkXmlFileType getFileType()
+ {
+ if (reader != 0)
+ {
+ // XXX: assume structured grid for now
+ return VTK_FILE_STRUCTURED_GRID;
+ }
+ return VTK_FILE_NONE;
+ }
+
+public:
+ VtkXmlReader();
+ ~VtkXmlReader();
+
+public:
+ int open(std::string filename);
+ void close();
+
+public:
+ void get_dataset(const char *loc, double **data, int *num, int *dim);
+ void get_coordinates(const char *log, double **coordinates, int *nno, int *nsd);
+ void get_connectivity(const char *log, int **connectivity, int *nel, int *ndofs);
+
+
+public:
+ void get_coordinates(double **coordinates, int *nno, int *nsd);
+ void get_connectivity(int **connectivity, int *nel, int *ndofs);
+ void get_vector_point_data(const char *name, double **vectors, int *num, int *dim);
+ void get_scalar_point_data(const char *name, double **scalars, int *num, int *dim);
+
+public:
+ vtkXMLStructuredGridReader *reader;
+ vtkStructuredGrid *grid;
+};
+
+#endif
Added: cs/benchmark/cigma/trunk/src/hold/Writer.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/Writer.cpp (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/Writer.cpp 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,46 @@
+#include <cassert>
+#include <cstdlib>
+#include "Writer.h"
+#include "HdfWriter.h"
+#include "TextWriter.h"
+#include "VtkWriter.h"
+
+
+using namespace cigma;
+
+
+// ---------------------------------------------------------------------------
+
+void new_writer(cigma::Writer **writer, std::string ext)
+{
+ if (ext == ".h5")
+ {
+ *writer = new HdfWriter();
+ return;
+ }
+
+ if (ext == ".txt")
+ {
+ *writer = new TextWriter();
+ return;
+ }
+
+ if (ext == ".vtk")
+ {
+ *writer = new VtkWriter();
+ return;
+ }
+
+}
+
+// ---------------------------------------------------------------------------
+
+cigma::Writer::Writer()
+{
+}
+
+cigma::Writer::~Writer()
+{
+}
+
+// ---------------------------------------------------------------------------
Added: cs/benchmark/cigma/trunk/src/hold/Writer.h
===================================================================
--- cs/benchmark/cigma/trunk/src/hold/Writer.h (rev 0)
+++ cs/benchmark/cigma/trunk/src/hold/Writer.h 2008-03-19 19:00:15 UTC (rev 11480)
@@ -0,0 +1,43 @@
+#ifndef __WRITER_H__
+#define __WRITER_H__
+
+#include <cstdio>
+#include <string>
+
+#include "FE_Field.h"
+
+
+namespace cigma
+{
+ class Writer;
+};
+
+
+class cigma::Writer
+{
+public:
+ typedef enum {
+ NULL_WRITER,
+ HDF_WRITER,
+ VTK_WRITER,
+ TXT_WRITER
+ } WriterType;
+
+public:
+ Writer();
+ ~Writer();
+
+public:
+ virtual WriterType getType() = 0;
+ virtual int open(std::string filename) = 0;
+ virtual void close() = 0;
+
+public:
+ virtual void write_field(FE_Field *field) = 0;
+};
+
+
+void new_writer(cigma::Writer **writer, std::string ext);
+
+
+#endif
More information about the cig-commits
mailing list