[cig-commits] r14095 - cs/cigma/trunk/src
luis at geodynamics.org
luis at geodynamics.org
Wed Feb 18 08:15:00 PST 2009
Author: luis
Date: 2009-02-18 08:15:00 -0800 (Wed, 18 Feb 2009)
New Revision: 14095
Modified:
cs/cigma/trunk/src/Locator.cpp
cs/cigma/trunk/src/MeshPart.cpp
cs/cigma/trunk/src/MeshPart.h
cs/cigma/trunk/src/cli_compare_cmd.cpp
cs/cigma/trunk/src/cli_compare_cmd.h
cs/cigma/trunk/src/core_compare_op.cpp
cs/cigma/trunk/src/io_exo_reader.cpp
cs/cigma/trunk/src/io_hdf5_reader.cpp
cs/cigma/trunk/src/io_vtk_reader.cpp
cs/cigma/trunk/src/loc_bbox.cpp
cs/cigma/trunk/src/loc_kdtree.cpp
Log:
Add options to control depth of the cell search algorithm
This only applies when function is associated with a finite element field.
Also, updated "cigma compare" to always include the summary
(unless a --quiet flag is given).
Modified: cs/cigma/trunk/src/Locator.cpp
===================================================================
--- cs/cigma/trunk/src/Locator.cpp 2009-02-18 16:14:58 UTC (rev 14094)
+++ cs/cigma/trunk/src/Locator.cpp 2009-02-18 16:15:00 UTC (rev 14095)
@@ -2,6 +2,7 @@
#include "tri_logger.hpp"
#include "loc_bbox.h"
#include "loc_kdtree.h"
+#include <sstream>
using namespace std;
using namespace cigma;
@@ -16,7 +17,16 @@
{
return Locator::KdtreeType;
}
- return Locator::BruteForceType;
+ else if (name == "none")
+ {
+ return Locator::BruteForceType;
+ }
+
+ // XXX: throw an exception for now. we should add a method
+ // to query for valid names instead.
+ std::ostringstream stream;
+ stream << "Unknown locator type '" << name << "'" << std::ends;
+ throw cigma::Exception("Locator::string2type", stream.str());
}
boost::shared_ptr<Locator> Locator::New(Locator::LocatorType loc_type)
Modified: cs/cigma/trunk/src/MeshPart.cpp
===================================================================
--- cs/cigma/trunk/src/MeshPart.cpp 2009-02-18 16:14:58 UTC (rev 14094)
+++ cs/cigma/trunk/src/MeshPart.cpp 2009-02-18 16:15:00 UTC (rev 14095)
@@ -8,6 +8,7 @@
#include <cassert>
#include <string>
#include <sstream>
+#include <algorithm>
using namespace std;
using namespace cigma;
using boost::shared_ptr;
@@ -79,6 +80,21 @@
}
}
+void MeshPart::setLocator(Locator::LocatorType loc_type)
+{
+ if (loc_type != Locator::BruteForceType)
+ {
+ // XXX: how do we handle the case when ncells is very small?
+ const int ncells = this->n_cells();
+ const int min_nnk = std::min(BoundingBoxLocator::default_nnk, KdtreeLocator::default_nnk);
+ if (ncells > min_nnk)
+ {
+ locator = Locator::New(loc_type);
+ locator->init(*this);
+ }
+ }
+}
+
void MeshPart::setCellType(Cell::type cell_type)
{
this->cell_type = cell_type;
Modified: cs/cigma/trunk/src/MeshPart.h
===================================================================
--- cs/cigma/trunk/src/MeshPart.h 2009-02-18 16:14:58 UTC (rev 14094)
+++ cs/cigma/trunk/src/MeshPart.h 2009-02-18 16:15:00 UTC (rev 14095)
@@ -51,6 +51,7 @@
void setElementBlock(const boost::shared_ptr<ElementBlock>& eb);
void setLocator(const boost::shared_ptr<Locator>& loc);
void setLocator();
+ void setLocator(Locator::LocatorType loc_type);
void setCellType(Cell::type cell_type);
/* incremental setters */
Modified: cs/cigma/trunk/src/cli_compare_cmd.cpp
===================================================================
--- cs/cigma/trunk/src/cli_compare_cmd.cpp 2009-02-18 16:14:58 UTC (rev 14094)
+++ cs/cigma/trunk/src/cli_compare_cmd.cpp 2009-02-18 16:15:00 UTC (rev 14095)
@@ -19,6 +19,16 @@
brief = "Calculate local and global residuals over two fields";
usage = "Usage: cigma compare <FIRST-FIELD> <SECOND-FIELD> -o <RESIDUAL-FILE>";
appendix = "";
+
+ nnk[0] = 0;
+ nnk[1] = 0;
+ nnk[2] = 0;
+
+ h[0] = 0;
+ h[1] = 0;
+ h[2] = 0;
+
+ threshold = std::numeric_limits<double>::infinity();
}
CompareCmd::~CompareCmd()
@@ -31,11 +41,11 @@
{
BaseCmd::defineOptions();
- opts_section = "Compare Options: ";
opts.add_options()
("first,a", po::value<string>(&first), "Path to first field")
("second,b", po::value<string>(&second), "Path to second field")
("output,o", po::value<string>(&outputfile), "Output file (for residuals)")
+ ("timer-frequency,t", po::value<int>(&(op.freq)), "Output frequency of timer")
;
opts2.add_options()
@@ -54,11 +64,21 @@
("rule,r", po::value<string>(&rule), "Integration rule (points and weights)")
("rule-weights", po::value<string>(&qw), "...integration weights on reference cell")
("rule-points", po::value<string>(&qx), "...integration points on reference cell")
+ ("global-threshold,g", po::value<double>(&threshold), "Threshold value for global residual")
("do-not-overwrite", "Do NOT overwrite existing data")
- ("timer-frequency,t", po::value<int>(&(op.freq)), "Output frequency of timer")
- ("global-threshold,g", po::value<double>(&threshold), "Threshold value for global residual")
;
+ other.add_options()
+ ("quiet", "Suppress all output")
+ ;
+
+ hidden.add_options()
+ ("first-mesh-locator", po::value<string>(&loc1))
+ ("first-mesh-nnk", po::value<int>(&nnk[1]))
+ ("second-mesh-locator", po::value<string>(&loc2))
+ ("second-mesh-nnk", po::value<int>(&nnk[2]))
+ ;
+
// define positional arguments
args.add("first", 1);
args.add("second", 1);
@@ -70,16 +90,6 @@
TRI_LOG_STR("CompareCmd::configure()");
- if (vm.count("verbose"))
- {
- op.verbose = true;
- }
-
- if (!vm.count("timer-frequency"))
- {
- op.freq = 1000;
- }
-
if (!vm.count("first"))
{
string msg("No first field was specified (use --first option)");
@@ -92,12 +102,6 @@
throw cigma::Exception("CompareCmd::configure", msg, 1);
}
- //if (!vm.count("output"))
- //{
- // string msg("No output file for residuals was specified (use --output option)");
- // throw cigma::Exception("CompareCmd::configure", msg, 1);
- //}
-
if (!vm.count("global-threshold"))
{
threshold = std::numeric_limits<double>::infinity();
@@ -111,13 +115,31 @@
}
}
+ if (!vm.count("timer-frequency"))
+ {
+ op.freq = 1000;
+ }
+
+ if (vm.count("quiet"))
+ {
+ op.verbose = false;
+ }
+ else if (vm.count("verbose"))
+ {
+ op.verbose = true;
+ }
+ else
+ {
+ op.verbose = false;
+ }
+
op.first_info.fn_name = first;
op.first_info.field_info.p_field = DataPath(first);
op.first_info.field_info.mesh_info.p_mesh = DataPath(mesh1);
op.first_info.field_info.mesh_info.p_nc = DataPath(nc1);
op.first_info.field_info.mesh_info.p_eb = DataPath(eb1);
op.first_info.field_info.mesh_info.cell_type_name = cell1;
- op.first_info.field_info.mesh_info.use_locator = true;
+ op.first_info.field_info.mesh_info.use_locator = false;
op.first_info.field_info.fe_info.q_info.p_quadrature = DataPath(qr1);
op.first_info.field_info.fe_info.q_info.p_weights = DataPath(qw1);
op.first_info.field_info.fe_info.q_info.p_points = DataPath(qx1);
@@ -129,7 +151,7 @@
op.second_info.field_info.mesh_info.p_nc = DataPath(nc2);
op.second_info.field_info.mesh_info.p_eb = DataPath(eb2);
op.second_info.field_info.mesh_info.cell_type_name = cell2;
- op.second_info.field_info.mesh_info.use_locator = true;
+ op.second_info.field_info.mesh_info.use_locator = false;
op.second_info.field_info.fe_info.q_info.p_quadrature = DataPath(qr2);
op.second_info.field_info.fe_info.q_info.p_weights = DataPath(qw2);
op.second_info.field_info.fe_info.q_info.p_points = DataPath(qx2);
@@ -146,6 +168,83 @@
op.domain_info.fe_info.q_info.cell_type_name = cell;
op.configure();
+
+ if (op.domain)
+ {
+ if ((op.domain)->mesh)
+ {
+ h[0] = (op.domain)->mesh->getMaxCellDiameter();
+ }
+ }
+ else
+ {
+ const string msg("Error: Could not configure integration mesh!");
+ throw cigma::Exception("CompareCmd::configure", msg);
+ }
+
+ if (op.first)
+ {
+ if ((op.first)->getType() == Function::FieldType)
+ {
+ Field *field = static_cast<Field*>(&(*(op.first)));
+
+ if (field && field->mesh)
+ {
+ if (vm.count("first-mesh-locator"))
+ {
+ field->mesh->setLocator(Locator::string2type(loc1));
+ }
+ else
+ {
+ field->mesh->setLocator();
+ }
+
+ if (field->mesh->locator && (nnk[1] > 0))
+ {
+ field->mesh->locator->setMaxNeighbors(nnk[1]);
+ }
+
+ h[1] = field->mesh->getMaxCellDiameter();
+ }
+ }
+ }
+ else
+ {
+ const string msg("Error: Could not configure first function!");
+ throw cigma::Exception("CompareCmd::configure", msg);
+ }
+
+ if (op.second)
+ {
+ if ((op.second)->getType() == Function::FieldType)
+ {
+ Field *field = static_cast<Field*>(&(*(op.second)));
+
+ if (field && field->mesh)
+ {
+ if (vm.count("second-mesh-locator"))
+ {
+ field->mesh->setLocator(Locator::string2type(loc2));
+ }
+ else
+ {
+ field->mesh->setLocator();
+ }
+
+ if (field->mesh->locator && (nnk[2] > 0))
+ {
+ field->mesh->locator->setMaxNeighbors(nnk[2]);
+ }
+
+ h[2] = field->mesh->getMaxCellDiameter();
+ }
+ }
+ }
+ else
+ {
+ const string msg("Error: Could not configure second function!");
+ throw cigma::Exception("CompareCmd::configure", msg);
+ }
}
int CompareCmd::run()
@@ -157,16 +256,15 @@
DataPath path(outputfile);
+ /* run the comparison */
int status = op.run();
+ /* gather results */
+ double L2 = op.residuals->L2();
double Linf = op.residuals->infinity_norm();
- double L2 = op.residuals->L2();
- double relative_L2 = op.residuals->relative_L2();
double volume = op.residuals->mesh->getVolume();
double L2_vol = L2/volume;
- double h = op.residuals->mesh->getMaxCellDiameter();
- //double h1 = op.first->getMaxCellDiameter();
- //double h2 = op.second->getMaxCellDiameter();
+ double L2_sqrt_vol = L2/sqrt(volume);
if (L2 > threshold)
{
@@ -177,31 +275,51 @@
status = 1;
}
- if (op.verbose)
+ if (!vm.count("quiet"))
{
+ // report global error metrics
+
+ const string indent(" ");
cout << setprecision(12);
-
- // report global error metrics
+
cout << endl;
cout << "Summary of comparison: " << endl;
- cout << " L2 = " << L2 << endl;
- cout << " Linf = " << Linf << endl;
- cout << " L2/volume = " << L2_vol << endl;
- cout << " sqrt((L2)^2/volume) = " << relative_L2 << endl;
- cout << " total volume of domain = " << volume << endl;
+ cout << indent << "L2 = " << L2 << endl;
+ cout << indent << "Linf = " << Linf << endl;
+ cout << indent << "volume = " << volume << endl;
+ cout << indent << "L2/volume = " << L2_vol << endl;
+ cout << indent << "L2/sqrt(volume) = " << L2_sqrt_vol << endl;
+
+ if (h[1] > 0)
+ {
+ cout << indent << "h1 = " << h[1] << endl;
+ }
+
+ if (h[2] > 0)
+ {
+ cout << indent << "h2 = " << h[2] << endl;
+ }
+
cout << endl;
-
}
if (vm.count("output"))
{
if (op.verbose)
{
- // write out data
- cout << "Creating file " << path.filename() << endl;
+ if (!path.exists())
+ {
+ cout << "Creating ";
+ }
+ else
+ {
+ cout << "Updating ";
+ }
+ cout << "'" << path << "'" << endl;
}
const bool normalize = false;
- WriteResiduals(path, op.residuals, normalize, !vm.count("do-not-overwrite"));
+ const bool overwrite = !vm.count("do-not-overwrite");
+ WriteResiduals(path, op.residuals, normalize, overwrite);
}
return status;
Modified: cs/cigma/trunk/src/cli_compare_cmd.h
===================================================================
--- cs/cigma/trunk/src/cli_compare_cmd.h 2009-02-18 16:14:58 UTC (rev 14094)
+++ cs/cigma/trunk/src/cli_compare_cmd.h 2009-02-18 16:15:00 UTC (rev 14095)
@@ -10,7 +10,7 @@
class CompareCmd;
}
-class cigma::CompareCmd : public BaseCmd
+class cigma::CompareCmd : public cigma::BaseCmd
{
public:
@@ -24,9 +24,13 @@
public:
std::string mesh, nc, eb, rule, qw, qx, cell;
- std::string first, mesh1, nc1, eb1, qr1, qw1, qx1, cell1;
- std::string second, mesh2, nc2, eb2, qr2, qw2, qx2, cell2;
+ std::string first, mesh1, loc1, nc1, eb1, qr1, qw1, qx1, cell1;
+ std::string second, mesh2, loc2, nc2, eb2, qr2, qw2, qx2, cell2;
+
+ int nnk[3];
+ double h[3];
double threshold;
+
std::string outputfile;
CompareOp op;
Modified: cs/cigma/trunk/src/core_compare_op.cpp
===================================================================
--- cs/cigma/trunk/src/core_compare_op.cpp 2009-02-18 16:14:58 UTC (rev 14094)
+++ cs/cigma/trunk/src/core_compare_op.cpp 2009-02-18 16:15:00 UTC (rev 14095)
@@ -120,13 +120,13 @@
throw cigma::Exception("CompareOp::configure", msg);
}
}
- assert(domain);
- assert(domain->mesh);
- assert(domain->mesh->coords);
- assert(domain->mesh->connect);
- assert(domain->fe);
- assert(domain->fe->cell);
- assert(domain->fe->quadrature);
+ //assert(domain);
+ //assert(domain->mesh);
+ //assert(domain->mesh->coords);
+ //assert(domain->mesh->connect);
+ //assert(domain->fe);
+ //assert(domain->fe->cell);
+ //assert(domain->fe->quadrature);
// Special cases for first function
if (dynamic_cast<cigma::ZeroFn*>(&(*first)))
@@ -215,12 +215,6 @@
residuals->setMesh(domain->mesh);
residuals->reset();
- //if (mesh->volume.size() == 0)
- //{
- // mesh->computeCellVolumes();
- //}
- //valarray<double>& vol = mesh->volume;
-
if (verbose)
{
cout << timer.header("cells");
Modified: cs/cigma/trunk/src/io_exo_reader.cpp
===================================================================
--- cs/cigma/trunk/src/io_exo_reader.cpp 2009-02-18 16:14:58 UTC (rev 14094)
+++ cs/cigma/trunk/src/io_exo_reader.cpp 2009-02-18 16:15:00 UTC (rev 14095)
@@ -237,7 +237,7 @@
mesh->coords = this->getNodeCoordinates("coords");
mesh->connect = this->getElementBlock("connect1");
mesh->cell_type = this->getCellType("connect1");
- mesh->setLocator();
+ //mesh->setLocator();
return mesh;
}
Modified: cs/cigma/trunk/src/io_hdf5_reader.cpp
===================================================================
--- cs/cigma/trunk/src/io_hdf5_reader.cpp 2009-02-18 16:14:58 UTC (rev 14094)
+++ cs/cigma/trunk/src/io_hdf5_reader.cpp 2009-02-18 16:15:00 UTC (rev 14095)
@@ -238,7 +238,7 @@
}
}
- mesh->setLocator();
+ //mesh->setLocator();
}
else
{
Modified: cs/cigma/trunk/src/io_vtk_reader.cpp
===================================================================
--- cs/cigma/trunk/src/io_vtk_reader.cpp 2009-02-18 16:14:58 UTC (rev 14094)
+++ cs/cigma/trunk/src/io_vtk_reader.cpp 2009-02-18 16:15:00 UTC (rev 14095)
@@ -287,7 +287,7 @@
mesh->connect = this->getElementBlock(0);
mesh->cell_type = this->getCellType(0);
- mesh->setLocator();
+ //mesh->setLocator();
return mesh;
}
Modified: cs/cigma/trunk/src/loc_bbox.cpp
===================================================================
--- cs/cigma/trunk/src/loc_bbox.cpp 2009-02-18 16:14:58 UTC (rev 14094)
+++ cs/cigma/trunk/src/loc_bbox.cpp 2009-02-18 16:15:00 UTC (rev 14095)
@@ -123,6 +123,8 @@
void BoundingBoxLocator::setMaxNeighbors(int k)
{
+ TRI_LOG_STR("BoundingBoxLocator::setMaxNeighbors()");
+
assert(k > 0);
if (nnIdx != 0) { delete [] nnIdx; }
@@ -131,6 +133,8 @@
nnk = k;
nnIdx = new ANNidx[nnk];
nnDists = new ANNdist[nnk];
+
+ TRI_LOG(nnk);
}
void BoundingBoxLocator::search(double *point)
Modified: cs/cigma/trunk/src/loc_kdtree.cpp
===================================================================
--- cs/cigma/trunk/src/loc_kdtree.cpp 2009-02-18 16:14:58 UTC (rev 14094)
+++ cs/cigma/trunk/src/loc_kdtree.cpp 2009-02-18 16:15:00 UTC (rev 14095)
@@ -81,7 +81,7 @@
kdtree = new ANNkd_tree(dataPoints, npts, ndim);
- searchRadius = 3 * mesh.getMaxCellDiameter();
+ searchRadius = 2.5 * mesh.getMaxCellDiameter();
TRI_LOG(npts);
TRI_LOG(ndim);
@@ -91,10 +91,12 @@
void KdtreeLocator::setMaxNeighbors(int k)
{
+ TRI_LOG_STR("KdtreeLocator::setMaxNeighbors()");
assert(k > 0);
nnk = k;
nnIdx.resize(nnk);
nnDists.resize(nnk);
+ TRI_LOG(nnk);
}
void KdtreeLocator::search(double *point)
More information about the CIG-COMMITS
mailing list