[cig-commits] r11721 - cs/benchmark/cigma/trunk/src
luis at geodynamics.org
luis at geodynamics.org
Wed Apr 2 11:01:38 PDT 2008
Author: luis
Date: 2008-04-02 11:01:38 -0700 (Wed, 02 Apr 2008)
New Revision: 11721
Modified:
cs/benchmark/cigma/trunk/src/QuadratureReader.cpp
cs/benchmark/cigma/trunk/src/QuadratureReader.h
Log:
Improvements to QuadratureReader
Modified: cs/benchmark/cigma/trunk/src/QuadratureReader.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/QuadratureReader.cpp 2008-04-02 18:01:36 UTC (rev 11720)
+++ cs/benchmark/cigma/trunk/src/QuadratureReader.cpp 2008-04-02 18:01:38 UTC (rev 11721)
@@ -2,6 +2,7 @@
#include <cassert>
#include <cstdlib>
#include "QuadratureReader.h"
+#include "Reader.h"
#include "StringUtils.h"
using namespace std;
@@ -11,8 +12,9 @@
QuadratureReader::QuadratureReader()
{
+ quadrature = 0;
+ meshPart = 0;
verbose = false;
- quadrature = 0;
}
QuadratureReader::~QuadratureReader()
@@ -24,6 +26,7 @@
}
}
+
// ---------------------------------------------------------------------------
void QuadratureReader::load_args(AnyOption *opt, const char *opt_prefix)
@@ -31,36 +34,40 @@
assert(opt != 0);
string prefix = opt_prefix;
- string opt_name;
+ string optstr;
char *in;
- opt_name = prefix + "-order";
- in = opt->getValue(opt_name.c_str());
+ optstr = prefix + "-order";
+ in = opt->getValue(optstr.c_str());
if (in != 0)
{
this->quadratureOrder = in;
}
- opt_name = prefix;
- in = opt->getValue(opt_name.c_str());
+ optstr = prefix;
+ in = opt->getValue(optstr.c_str());
if (in != 0)
{
this->quadraturePath = in;
}
- opt_name = prefix + "-points";
- in = opt->getValue(opt_name.c_str());
+ optstr = prefix + "-points";
+ in = opt->getValue(optstr.c_str());
if (in != 0)
{
this->pointsPath = in;
}
- opt_name = prefix + "-weights";
- in = opt->getValue(opt_name.c_str());
+ optstr = prefix + "-weights";
+ in = opt->getValue(optstr.c_str());
if (in != 0)
{
this->weightsPath = in;
}
+
+ optstr = prefix + "-mesh";
+ meshPartReader.load_args(opt, optstr.c_str());
+
}
void QuadratureReader::validate_args(const char *cmd_name)
@@ -105,16 +112,55 @@
<< endl;
exit(1);
}
+
+ meshPartReader.validate_args(cmd_name);
}
// ---------------------------------------------------------------------------
-void QuadratureReader::load_quadrature(Cell *cell)
+void QuadratureReader::set_mesh(MeshPart *meshPart)
{
- assert(cell != 0);
+ //
+ // Provide a default mesh just in case none is explicitly given.
+ // Note that load_mesh() will overwrite meshPart, and will quit
+ // if no default has been provided.
+ //
+ assert(this->meshPartReader.meshPart == 0);
+ this->meshPart = meshPart;
+}
+void QuadratureReader::load_mesh()
+{
+ //
+ // Load integration mesh, or quit.
+ //
+ meshPartReader.load_mesh();
+ if (meshPartReader.meshPart != 0)
+ {
+ // integration mesh loaded, so override default
+ meshPart = meshPartReader.meshPart;
+ }
+ if (meshPart == 0)
+ {
+ // no mesh found, so quit
+ cerr << "QuadratureReader: Could not load integration mesh" << endl;
+ exit(1);
+ }
+}
+
+void QuadratureReader::load_quadrature()
+{
+ //
+ // Load quadrature mesh
+ //
+ this->load_mesh();
+ assert(meshPart != 0);
+ assert(meshPart->cell != 0);
+
+ //
+ // Now, load quadrature points and weights on reference cell
+ //
int i,j;
-
int nq, nd;
double *qx, *qw;
@@ -131,6 +177,7 @@
ierr = pointsReader->open(pointsFile.c_str());
if (ierr < 0)
{
+ cerr << "QuadratureReader: ";
cerr << "Could not open quadrature points file " << pointsFile << endl;
exit(1);
}
@@ -138,6 +185,7 @@
ierr = pointsReader->get_dataset(pointsLoc.c_str(), &qx, &nq, &nd);
if (ierr < 0)
{
+ cerr << "QuadratureReader: ";
cerr << "Could not access quadrature points in file " << pointsFile << endl;
exit(1);
}
@@ -154,6 +202,7 @@
ierr = weightsReader->open(weightsFile.c_str());
if (ierr < 0)
{
+ cerr << "QuadratureReader: ";
cerr << "Could not open quadrature weights file " << weightsFile << endl;
exit(1);
}
@@ -161,16 +210,19 @@
ierr = weightsReader->get_dataset(weightsLoc.c_str(), &qw, &wtdims[0], &wtdims[1]);
if (ierr < 0)
{
+ cerr << "QuadratureReader: ";
cerr << "Could not access quadrature weights in file " << weightsFile << endl;
exit(1);
}
if (wtdims[0] != nq)
{
+ cerr << "QuadratureReader: ";
cerr << "Number of points and weights do not match!" << endl;
exit(1);
}
if (wtdims[1] != 1)
{
+ cerr << "QuadratureReader: ";
cerr << "Array of weights is not one-dimensionsal!" << endl;
exit(1);
}
@@ -182,7 +234,8 @@
int ierr;
string quadratureLoc, quadratureFile, quadratureExt;
parse_dataset_path(quadraturePath, quadratureLoc, quadratureFile, quadratureExt);
- reader = NewReader(quadratureExt.c_str());
+
+ Reader *reader = NewReader(quadratureExt.c_str());
ierr = reader->open(quadratureFile.c_str());
if (ierr < 0)
{
@@ -198,6 +251,7 @@
ierr = reader->get_dataset(pointsLoc.c_str(), &qx, &nq, &nd);
if (ierr < 0)
{
+ cerr << "QuadratureReader: ";
cerr << "Could not read quadrature points from location " << pointsLoc << endl;
exit(1);
}
@@ -208,16 +262,19 @@
ierr = reader->get_dataset(weightsLoc.c_str(), &qw, &wtdims[0], &wtdims[1]);
if (ierr < 0)
{
+ cerr << "QuadratureReader: ";
cerr << "Could not read quadrature weights from location " << weightsLoc << endl;
exit(1);
}
if (wtdims[0] != nq)
{
+ cerr << "QuadratureReader: ";
cerr << "Number of points and weights do not match!" << endl;
exit(1);
}
if (wtdims[1] != 1)
{
+ cerr << "QuadratureReader: ";
cerr << "Array of weights is not one-dimensional!" << endl;
exit(1);
}
@@ -227,6 +284,9 @@
cerr << "Cannot read explicit quadrature rule from a (" << quadratureExt << ") file." << endl;
exit(1);
}
+
+ reader->close();
+ delete reader;
}
else if (quadratureOrder != "")
{
@@ -244,59 +304,84 @@
delete [] qx;
delete [] qw;
// */
+ assert(false);
}
}
else
{
- // assign reasonable defaults
- quadrature = new QuadraturePoints();
- cell->qr_default(&qw, &qx, &nq, &nd);
- quadrature->set_quadrature(qw, qx, nq, nd);
+ // load default quadrature rule
+ meshPart->cell->qr_default(&qw, &qx, &nq, &nd);
}
if ((qx != 0) && (qw != 0))
{
- quadrature = new QuadraturePoints();
+ quadrature = new QuadratureRule();
quadrature->set_quadrature(qw, qx, nq, nd);
+ quadrature->set_mesh(meshPart);
}
assert(quadrature != 0);
- assert(quadrature->n_points() > 0);
- assert(quadrature->n_dim() > 0);
+ assert(quadrature->points != 0);
+ assert(quadrature->points->n_points() > 0);
+ assert(quadrature->points->n_dim() > 0);
- cout << "quadrature rule = "
- << quadrature->n_points()
- << " points on a "
- << quadrature->n_dim()
- << "-d cell"
- << endl;
if (verbose)
{
- for (i = 0; i < quadrature->n_points(); i++)
+ QuadraturePoints *qpoints = quadrature->points;
+
+ cout << endl
+ << "quadrature rule = "
+ << qpoints->n_points()
+ << " points on a "
+ << qpoints->n_dim()
+ << "-d cell"
+ << endl;
+
+ for (i = 0; i < qpoints->n_points(); i++)
{
cout << "\t";
- cout << quadrature->weight(i);
+ cout << qpoints->weight(i);
cout << "\t";
- for (j = 0; j < quadrature->n_dim(); j++)
+ for (j = 0; j < qpoints->n_dim(); j++)
{
- cout << quadrature->refpoint(i,j) << " ";
+ cout << qpoints->refpoint(i,j) << " ";
}
cout << endl;
}
}
+ this->check_interior();
+ this->warn_on_negative_weights();
+
+ //
+ // Clean up any local allocations
+ //
+ if (qx != 0) delete [] qx;
+ if (qw != 0) delete [] qw;
+}
+
+void QuadratureReader::check_interior()
+{
+ assert(quadrature != 0);
+ assert(meshPart != 0);
+ assert(meshPart->cell != 0);
+
//
// now we need to ensure that every point in our
// quadrature rule is indeed contained inside our cell
//
+ int i,j;
bool all_contained = true;
- for (i = 0; i < quadrature->n_points(); i++)
+ QuadraturePoints *qpoints = quadrature->points;
+ MeshPart *meshPart = meshPart;
+ Cell *cell = meshPart->cell;
+ for (int i = 0; i < qpoints->n_points(); i++)
{
double refpt[3] = {0,0,0};
- for (j = 0; j < quadrature->n_dim(); j++)
+ for (int j = 0; j < qpoints->n_dim(); j++)
{
- refpt[j] = quadrature->refpoint(i,j);
+ refpt[j] = qpoints->refpoint(i,j);
}
if (!cell->interior(refpt[0], refpt[1], refpt[2]))
{
@@ -308,19 +393,28 @@
{
// XXX: add more information here (e.g., actual reference cell used,
// such as tri3, quad4, tet4, hex8, ...)
- cerr << "QuadratureReader error: "
- << "Some quadrature points lie outside required reference cell"
+ cerr << "QuadratureReader: ";
+ cerr << "Some quadrature points lie outside required reference cell"
<< endl;
exit(1);
}
+}
+
+void QuadratureReader::warn_on_negative_weights()
+{
+ assert(quadrature != 0);
+ assert(quadrature->points != 0);
+
//
// emit warning if there are any negative quadrature weights?
//
+ int i;
bool some_neg_wts = false;
- for (i = 0; i < quadrature->n_points(); i++)
+ QuadraturePoints *qpoints = quadrature->points;
+ for (i = 0; i < qpoints->n_points(); i++)
{
- if (quadrature->weight(i) < 0)
+ if (qpoints->weight(i) < 0)
{
some_neg_wts = true;
break;
@@ -329,18 +423,12 @@
if (some_neg_wts)
{
//
- // XXX: only emit warning. exit(1) is not necessary...
+ // only emit warning. exit(1) is not necessary...
//
- cerr << "QuadratureReader warning: "
- << "Some weights are negative"
+ cerr << "QuadratureReader: ";
+ cerr << "Warning, some weights are negative"
<< endl;
}
-
- //
- // Clean up any local allocations
- //
- if (qx != 0) delete [] qx;
- if (qw != 0) delete [] qw;
}
// ---------------------------------------------------------------------------
Modified: cs/benchmark/cigma/trunk/src/QuadratureReader.h
===================================================================
--- cs/benchmark/cigma/trunk/src/QuadratureReader.h 2008-04-02 18:01:36 UTC (rev 11720)
+++ cs/benchmark/cigma/trunk/src/QuadratureReader.h 2008-04-02 18:01:38 UTC (rev 11721)
@@ -3,10 +3,11 @@
#include <string>
#include "AnyOption.h"
-#include "Reader.h"
-#include "QuadraturePoints.h"
-#include "Cell.h"
+//#include "Cell.h"
+#include "MeshPartReader.h"
+#include "QuadratureRule.h"
+
class QuadratureReader
{
public:
@@ -18,9 +19,15 @@
void validate_args(const char *cmd_name);
public:
- void load_quadrature(cigma::Cell *cell);
+ void set_mesh(cigma::MeshPart *meshPart);
+ void load_mesh();
public:
+ void load_quadrature();
+ void check_interior();
+ void warn_on_negative_weights();
+
+public:
std::string quadratureOrder;
std::string quadraturePath;
std::string pointsPath;
@@ -28,8 +35,10 @@
bool verbose;
public:
- cigma::Reader *reader;
- cigma::QuadraturePoints *quadrature;
+ cigma::QuadratureRule *quadrature;
+ cigma::MeshPart *meshPart;
+ MeshPartReader meshPartReader;
};
+
#endif
More information about the cig-commits
mailing list