[cig-commits] r13720 - cs/cigma/trunk/src
luis at geodynamics.org
luis at geodynamics.org
Wed Dec 17 02:32:05 PST 2008
Author: luis
Date: 2008-12-17 02:32:05 -0800 (Wed, 17 Dec 2008)
New Revision: 13720
Modified:
cs/cigma/trunk/src/Quadrature.cpp
Log:
Generalized Quadrature::setPath() init method
* Look for a CellType attribute in the right place
Modified: cs/cigma/trunk/src/Quadrature.cpp
===================================================================
--- cs/cigma/trunk/src/Quadrature.cpp 2008-12-17 10:32:03 UTC (rev 13719)
+++ cs/cigma/trunk/src/Quadrature.cpp 2008-12-17 10:32:05 UTC (rev 13720)
@@ -1,5 +1,6 @@
#include "Quadrature.h"
#include "io_file_reader.h"
+#include "io_hdf5_reader.h"
#include "core_array.h"
#include "core_readers.h"
#include <iostream>
@@ -142,7 +143,7 @@
shared_ptr<Quadrature> Quadrature::New(const QuadratureInfo& q_info)
{
- TRI_LOG_STR("Quadrature::New");
+ TRI_LOG_STR("Quadrature::New()");
TRI_LOG(q_info);
shared_ptr<Quadrature> Q;
@@ -167,6 +168,7 @@
{
Cell::type cell_type = Cell::string2type(q_info.cell_type_name);
Q = Quadrature::default_rule(cell_type);
+ return Q;
}
return Q;
@@ -271,55 +273,132 @@
void Quadrature::setPath(const DataPath& p_quad)
{
- TRI_LOG_STR("Quadrature::setPath");
+ TRI_LOG_STR("Quadrature::setPath()");
TRI_LOG(p_quad);
// XXX: check that location points to HDF5 group
// XXX: check for cell type attribute
- // load weights & points simultaneously
- cigma::array<double> *wx = ReadArray(p_quad);
+ string filename = p_quad.filename();
+ string location = p_quad.location();
+ shared_ptr<FileReader> reader = FileReader::New(filename, "r");
+ const FileReader::ReaderType rt = reader->getReaderType();
- if (wx == 0)
+ bool combined_weights_and_points = false;
+
+ if (rt == FileReader::HDF5_FILE_READER)
{
+ HDF5_Reader* h5 = static_cast<HDF5_Reader*>(&(*reader));
+
+ bool locIsGroup = h5->locationIsGroup(location.c_str());
+ TRI_LOG(locIsGroup);
+ if (locIsGroup)
+ {
+ DataPath p_wts, p_pts;
+
+ p_wts.set_filename(filename);
+ p_wts.set_location(location + "/weights");
+
+ p_pts.set_filename(filename);
+ p_pts.set_location(location + "/points");
+
+ // read CellType attribute, if necessary
+ if (cell_type == Cell::NONE)
+ {
+ string cell_name;
+ int status = h5->readAttrString(location.c_str(), "CellType", cell_name);
+ TRI_LOG(status);
+ if (status == 0)
+ {
+ TRI_LOG(cell_name);
+ this->cell_type = Cell::string2type(cell_name);
+ }
+ }
+
+ // load weights & points from separate arrays
+ this->setPath2(p_wts, p_pts);
+ }
+ else
+ {
+ bool locIsDataset = h5->locationIsDataset(location.c_str());
+ TRI_LOG(locIsDataset);
+ if (locIsDataset)
+ {
+ // load weights & points simultaneously
+ combined_weights_and_points = true;
+ }
+ else
+ {
+ std::ostringstream stream;
+ stream << "Could not read quadrature points from data path '"
+ << p_quad << "'" << std::ends;
+ throw cigma::Exception("Quadrature::setPath", stream.str());
+ }
+ }
+ }
+ else if ((rt == FileReader::TEXT_FILE_READER) || (rt == FileReader::VTK_FILE_READER))
+ {
+ // load weights & points simultaneously
+ combined_weights_and_points = true;
+ }
+ else
+ {
std::ostringstream stream;
- stream << "Failed to read quadrature from path '"
- << p_quad << "'" << std::ends;
+ stream << "Could not read file '" << filename << "'" << std::ends;
throw cigma::Exception("Quadrature::setPath", stream.str());
}
- const int npts = wx->n_points();
- const int ndim = -1 + wx->n_dim();
+ if (combined_weights_and_points)
+ {
+ cigma::array<double> *wx = ReadArray(p_quad);
- this->reinit(npts, ndim);
+ if (wx == 0)
+ {
+ std::ostringstream stream;
+ stream << "Failed to read quadrature from path '"
+ << p_quad << "'" << std::ends;
+ throw cigma::Exception("Quadrature::setPath", stream.str());
+ }
- for (int q = 0; q < npts; q++)
- {
- this->setWeight(q, wx->data(q, 0));
- for (int j = 0; j < ndim; j++)
+ const int q_npts = wx->n_points();
+ const int q_ndim = -1 + wx->n_dim();
+
+ this->reinit(q_npts, q_ndim);
+
+ for (int q = 0; q < npts; q++)
{
- this->setPoint(q, j, wx->data(q, 1+j));
+ this->setWeight(q, wx->data(q, 0));
+ for (int j = 0; j < ndim; j++)
+ {
+ this->setPoint(q, j, wx->data(q, 1+j));
+ }
}
+ delete wx;
+
+ // check points
+ bool all_points_in_cell = points_in_cell();
+ TRI_LOG(all_points_in_cell);
+ if (!all_points_in_cell)
+ {
+ std::ostringstream stream;
+ stream << "Quadrature points not compatible with cell '";
+ stream << Cell::type2string(cell_type);
+ stream << "'" << std::ends;
+ throw cigma::Exception("Quadrature::setPath", stream.str());
+ }
}
- delete wx;
- // check points
- if (!points_in_cell())
- {
- std::ostringstream stream;
- stream << "Quadrature points not compatible with cell '";
- stream << Cell::type2string(cell_type);
- stream << "'" << std::ends;
- throw cigma::Exception("Quadrature::setPath", stream.str());
- }
}
+
void Quadrature::setPath2(const DataPath& p_weights, const DataPath& p_points)
{
TRI_LOG_STR("Quadrature::setPath2");
TRI_LOG(p_weights);
TRI_LOG(p_points);
+ int q_nwts, q_npts, q_ndim;
+
// load weights array from p_weights
cigma::array<double> *wts = ReadArray(p_weights);
if (wts == 0)
@@ -329,6 +408,7 @@
<< p_weights << "'" << std::ends;
throw cigma::Exception("Quadrature::setPath2", stream.str());
}
+ q_nwts = wts->n_points();
this->weights = wts->_data;
wts->_data = 0;
delete wts;
@@ -342,12 +422,25 @@
<< p_points << "'" << std::ends;
throw cigma::Exception("Quadrature::setPath2", stream.str());
}
+ this->npts = pts->n_points();
+ this->ndim = pts->n_dim();
this->points = pts->_data;
pts->_data = 0;
delete pts;
+ if (q_nwts != this->npts)
+ {
+ std::ostringstream stream;
+ stream << "Number of weights (" << q_nwts << ") "
+ << "does not match the number of quadrature "
+ << "points (" << npts << ")" << std::ends;
+ throw cigma::Exception("Quadrature::setPath2", stream.str());
+ }
+
// check points
- if (!points_in_cell())
+ bool all_points_in_cell = points_in_cell();
+ TRI_LOG(all_points_in_cell);
+ if (!all_points_in_cell)
{
std::ostringstream stream;
stream << "Quadrature points not compatible with cell '";
@@ -392,12 +485,19 @@
bool Quadrature::points_in_cell() const
{
- int i,j;
- if (cell_type == Cell::NONE) { return false; }
+ if (cell_type == Cell::NONE)
+ {
+ string msg("Quadrature cell type is not set (cannot check interior)");
+ throw cigma::Exception("Quadrature::points_in_cell", msg);
+ }
shared_ptr<Cell> cell = Cell::New(cell_type);
- if (!cell) { return false; }
+ if (!cell)
+ {
+ return false;
+ }
+ int i,j;
for (i = 0; i < npts; i++)
{
double refpt[3] = {0,0,0};
More information about the CIG-COMMITS
mailing list