[cig-commits] r9211 - cs/benchmark/cigma/trunk/src
luis at geodynamics.org
luis at geodynamics.org
Fri Feb 1 12:52:08 PST 2008
Author: luis
Date: 2008-02-01 12:52:06 -0800 (Fri, 01 Feb 2008)
New Revision: 9211
Modified:
cs/benchmark/cigma/trunk/src/CompareCmd.cpp
cs/benchmark/cigma/trunk/src/Misc.cpp
cs/benchmark/cigma/trunk/src/Misc.h
Log:
Improvements to CompareCmd.cpp and Misc.cpp
* Now using generic IO loaders from Misc.cpp (VTK case works)
* Took advantage of virtual methods in Reader
* Removed more obsoleted code (previously commented out)
Modified: cs/benchmark/cigma/trunk/src/CompareCmd.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/CompareCmd.cpp 2008-02-01 20:52:04 UTC (rev 9210)
+++ cs/benchmark/cigma/trunk/src/CompareCmd.cpp 2008-02-01 20:52:06 UTC (rev 9211)
@@ -106,26 +106,31 @@
/*
* Initialization order:
- * Load Integration mesh
- * Load Quadrature rule
* Load First field
* If FE_Field
- * Load MeshA
- * Load DofsB
+ * Load DofsA
+ * Load MeshA if req'd
+ * Load RuleA if req'd
* Else
* Load Analytic Field
* Load Second field
* If FE_Field
- * Load MeshB
* 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 */
+ //XXX: add "rule" to last arg
+ //XXX: rename configure_quadrature to configure_rule
+
configure_mesh(opt, &meshIO, "mesh");
- configure_quadrature(opt, &quadratureIO);
+ configure_quadrature(opt, &quadratureIO, "rule");
configure_field(opt, &firstIO, "first");
configure_field(opt, &secondIO, "second");
configure_field(opt, &residualsIO, "output");
@@ -151,26 +156,22 @@
if (secondIO.field_path == "")
{
cerr << "compare: Please specify the option --second" << endl;
- exit(2);
+ exit(1);
}
if (residualsIO.field_path == "")
{
cerr << "compare: Please specify the option --output" << endl;
+ exit(1);
}
- /* XXX: if no mesh specified, get it from fieldA
- * if fieldA has no mesh (e.g. specified by analytic soln), then
- * swap fieldA and fieldB, and try again.
- * if fieldA still has no mesh, then produce error if no mesh
- * was specified to begin with.
- */
- meshIO.load();
- mesh = meshIO.meshPart;
+ /* 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, "
@@ -178,16 +179,11 @@
<< field_a->fe->cell->n_nodes() << " dofs/cell, "
<< "rank " << field_a->n_rank() << endl;
- if (mesh == 0)
- {
- mesh = firstIO.field->meshPart;
- }
- assert(mesh != 0);
-
-
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, "
@@ -196,6 +192,22 @@
<< "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);
+
+
quadratureIO.load(mesh->cell);
quadrature = quadratureIO.quadrature;
if (quadrature == 0)
@@ -204,7 +216,11 @@
}
assert(quadrature != 0);
+ // XXX: perhaps quadrature data should be loaded first!
+ // move this line back into FieldIO load() method
+ field_a->fe->set_quadrature(quadrature);
+
/* Determine the output-frequency option */
verbose = opt->getFlag("verbose");
@@ -290,6 +306,7 @@
field_a->tabulate();
// ... calculate phi_a[]
+ // XXX: use tabulation instead of calling eval repeatedly!
for (q = 0; q < nq; q++)
{
field_a->eval((*quadrature)[q], &phi_a[valdim*q]);
Modified: cs/benchmark/cigma/trunk/src/Misc.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/Misc.cpp 2008-02-01 20:52:04 UTC (rev 9210)
+++ cs/benchmark/cigma/trunk/src/Misc.cpp 2008-02-01 20:52:06 UTC (rev 9211)
@@ -87,7 +87,7 @@
// ---------------------------------------------------------------------------
-void configure_quadrature(AnyOption *opt, QuadratureIO *quadratureIO)
+void configure_quadrature(AnyOption *opt, QuadratureIO *quadratureIO, const char *opt_prefix)
{
assert(opt != 0);
assert(quadratureIO != 0);
@@ -171,113 +171,6 @@
// ---------------------------------------------------------------------------
-/*
-void read_double_dset(HdfReader *reader, double **data, int *num, int *dim)
-{
-}
-
-void read_double_dset(TextReader *reader, double **data, int *num, int *dim)
-{
-}
-
-void read_double_dset(VtkReader *reader, double **data, int *num, int *dim)
-{
-}
-
-void read_double_dset(Reader *reader, double **data, int *num, int *dim)
-{
-
- if (reader->getType() == Reader::HDF_READER)
- {
- meshPart = new MeshPart();
- HdfReader *hdfReader = static_cast<HdfReader*>(reader);
- read_double_dset(hdfReader, data, num, dim);
-
- }
- else if (reader->getType() == Reader::TXT_READER)
- {
- meshPart = new MeshPart();
- TextReader *textReader = static_cast<TextReader*>(reader);
- read_double_dset(textReader, data, num, dim);
-
- }
- else if (reader->getType() == Reader::VTK_READER)
- {
- meshPart = new MeshPart();
- VtkReader *vtkReader = static_cast<VtkReader*>(reader);
- read_double_dset(vtkReader, data, num, dim);
- }
-
-}
-*/
-
-// ---------------------------------------------------------------------------
-
-/*
-void read_int_dset(HdfReader *reader, int **data, int *num, int *dim)
-{
-}
-
-void read_int_dset(TextReader *reader, int **data, int *num, int *dim)
-{
-}
-
-void read_int_dset(VtkReader *reader, int **data, int *num, int *dim)
-{
-}
-
-void read_int_dset(Reader *reader, int **data, int *num, int *dim)
-{
-
- if (reader->getType() == Reader::HDF_READER)
- {
- meshPart = new MeshPart();
- HdfReader *hdfReader = static_cast<HdfReader*>(reader);
- read_int_dset(hdfReader, data, num, dim);
-
- }
- else if (reader->getType() == Reader::TXT_READER)
- {
- meshPart = new MeshPart();
- TextReader *textReader = static_cast<TextReader*>(reader);
- read_int_dset(textReader, data, num, dim);
-
- }
- else if (reader->getType() == Reader::VTK_READER)
- {
- meshPart = new MeshPart();
- VtkReader *vtkReader = static_cast<VtkReader*>(reader);
- read_int_dset(vtkReader, data, num, dim);
- }
-
-}
-*/
-
-
-// ---------------------------------------------------------------------------
-
-
-/*
-void read_coords(Reader *reader, double **coords, int *nno, int *nsd)
-{
- read_double_dset(reader, coords, nno, nsd);
-}
-
-void read_connect(Reader *reader, int **connect, int *nel, int *ndofs)
-{
- read_int_dset(reader, connect, nel, ndofs);
-}
-
-void read_dofs(Reader *reader, double **dofs, int *nno, int *valdim)
-{
- read_double_dset(reader, dofs, nno, valdim);
-}
-*/
-
-
-
-// ---------------------------------------------------------------------------
-
MeshIO::MeshIO()
{
meshPart = 0;
@@ -293,6 +186,8 @@
void MeshIO::load()
{
+ cout << "Calling MeshIO::load()" << endl;
+
string mesh_loc, mesh_file, mesh_ext;
string coords_loc, coords_file, coords_ext;
string connect_loc, connect_file, connect_ext;
@@ -306,133 +201,50 @@
nno = nsd = 0;
nel = ndofs = 0;
-
+ // XXX: use auto_ptr for the local readers, so we can throw exceptions
if (coords_path != "")
{
- // XXX: use auto_ptr for this local reader, so we can throw exceptions
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);
-
- //read_coords(coords_reader, &coords, &nno, &nsd);
-
- if (coords_reader->getType() == Reader::HDF_READER)
- {
- HdfReader *hdfReader = static_cast<HdfReader*>(coords_reader);
- hdfReader->get_coordinates(coords_loc.c_str(), &coords, &nno, &nsd);
- hdfReader->close();
- }
- else if (coords_reader->getType() == Reader::TXT_READER)
- {
- TextReader *textReader = static_cast<TextReader*>(coords_reader);
- textReader->get_coordinates(&coords, &nno, &nsd);
- textReader->close();
-
- }
- else if (coords_reader->getType() == Reader::VTK_READER)
- {
- VtkReader *vtkReader = static_cast<VtkReader*>(coords_reader);
- vtkReader->get_coordinates(&coords, &nno, &nsd);
- vtkReader->close();
- }
+ coords_reader->get_coordinates(coords_loc.c_str(), &coords, &nno, &nsd);
+ coords_reader->close();
}
if (connect_path != "")
{
- // XXX: use auto_ptr for this local reader, so we can throw exceptions
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);
-
- //read_connect(connect_reader, &nel, &ndofs);
-
- if (connect_reader->getType() == Reader::HDF_READER)
- {
- HdfReader *hdfReader = static_cast<HdfReader*>(connect_reader);
- hdfReader->get_coordinates(coords_loc.c_str(), &coords, &nno, &nsd);
- hdfReader->close();
- }
- else if (connect_reader->getType() == Reader::TXT_READER)
- {
- TextReader *textReader = static_cast<TextReader*>(connect_reader);
- textReader->get_coordinates(&coords, &nno, &nsd);
- textReader->close();
- }
- else if (connect_reader->getType() == Reader::VTK_READER)
- {
- VtkReader *vtkReader = static_cast<VtkReader*>(connect_reader);
- vtkReader->get_coordinates(&coords, &nno, &nsd);
- vtkReader->close();
- }
+ connect_reader->get_connectivity(connect_loc.c_str(), &connect, &nel, &ndofs);
+ connect_reader->close();
}
if ((mesh_path != "") && ((coords == 0) || (connect == 0)))
{
- // XXX: use auto_ptr for this local reader, so we can throw exceptions
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 (mesh_reader->getType() == Reader::HDF_READER)
+ if (coords == 0)
{
- coords_loc = mesh_loc + "/coordinates";
- connect_loc = mesh_loc + "/connectivity";
- HdfReader *hdfReader = static_cast<HdfReader*>(mesh_reader);
- hdfReader->get_coordinates(coords_loc.c_str(), &coords, &nno, &nsd);
- hdfReader->get_connectivity(connect_loc.c_str(), &connect, &nel, &ndofs);
+ if (mesh_reader->getType() == Reader::HDF_READER)
+ {
+ coords_loc = mesh_loc + "/coordinates";
+ }
+ mesh_reader->get_coordinates(coords_loc.c_str(), &coords, &nno, &nsd);
}
- else if (mesh_reader->getType() == Reader::TXT_READER)
+ if (connect == 0)
{
- TextReader *textReader = static_cast<TextReader*>(mesh_reader);
- textReader->get_coordinates(&coords, &nno, &nsd);
- textReader->get_connectivity(&connect, &nel, &ndofs);
+ if (mesh_reader->getType() == Reader::HDF_READER)
+ {
+ connect_loc = mesh_loc + "/connectivity";
+ }
+ mesh_reader->get_connectivity(connect_loc.c_str(), &connect, &nel, &ndofs);
}
- else if (mesh_reader->getType() == Reader::VTK_READER)
- {
- // read mesh from single vtk file
- VtkReader *vtkReader = static_cast<VtkReader*>(mesh_reader);
- vtkReader->get_coordinates(&coords, &nno, &nsd);
- vtkReader->get_connectivity(&connect, &nel, &ndofs);
- }
-
- reader->close();
-
- /*
- if (mesh_ext == ".h5")
- {
- assert(mesh_reader->getType() == Reader::HDF_READER);
-
- HdfReader *hdfReader = static_cast<HdfReader*>(mesh_reader);
- // open mesh_file in read only mode
-
- // assert that mesh_path points to an HDF5 group
-
- // read metadata from that group
-
- // if no metadata is available, look inside group
- // for two datasets named coordinates & connectivity
-
- }
- else if (mesh_ext == ".txt")
- {
- assert(reader->getType() == Reader::TXT_READER);
-
- // read mesh from single text file
- }
- else if (mesh_ext == ".vtk")
- {
- assert(reader->getType() == Reader::VTK_READER);
-
- // read mesh from single vtk file
- VtkReader *vtkReader = static_cast<VtkReader*>(mesh_reader);
- vtkReader->open(mesh_file);
- vtkReader->get_coordinates(&coords, &nno, &nsd);
- vtkReader->get_connectivity(&connect, &nel, &ndofs);
- vtkReader->close();
- } */
+ mesh_reader->close();
}
if ((coords != 0) && (connect != 0))
@@ -446,7 +258,6 @@
meshPart->nel = nel;
meshPart->ndofs = ndofs;
meshPart->connect = connect;
-
}
else
{
@@ -454,7 +265,6 @@
{
cerr << "MeshIO::load() error: Could not find mesh coordinates";
cerr << endl;
-
}
if (connect == 0)
{
@@ -463,7 +273,6 @@
}
}
-
}
@@ -485,6 +294,8 @@
void QuadratureIO::load(cigma::Cell *cell)
{
+ cout << "Calling QuadratureIO::load()" << endl;
+
assert(cell != 0);
// XXX: change *_nsd to *_celldim since we are
@@ -626,7 +437,7 @@
{
quadrature = new Quadrature();
quadrature->set_quadrature(qx, qw, nq, nd);
- quadrature->set_globaldim(3); // XXX
+ quadrature->set_globaldim(3); // XXX: how to treat 2D case?
}
assert(quadrature != 0);
@@ -654,163 +465,48 @@
void FieldIO::load()
{
+ cout << "Calling FieldIO::load()" << endl;
- int nno, nsd;
- int nel, ndofs;
- double *coords;
- int *connect;
-
int dofs_nno, dofs_valdim;
double *dofs;
+ string dofs_loc, dofs_file, dofs_ext;
if (field_path != "")
{
- assert(false);
- }
+ 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();
- /* XXX: For the following two cases, I need a static initializer on
- * Reader class that instantiates the right subclass based on the
- * detected filetype,
- *
- * After that, here I also need to use a switch statement on the value
- * of reader->getReaderType() in order to downcast the Reader
- * object properly...
- *
- * Having done that, I can rely on function polymorphism to call the right
- * routine: e.g choose between LoadField(VtkReader *reader, ...),
- * LoadField(HdfReader *reader, ...), LoadField(TextReader *reader)
- *
- * For detecting the filetype, one could rely only on the extension,
- * or possibly check for a magic number at the beginning of the file.
- */
- /*
- switch (reader->getType())
- {
- case Reader::HDF_READER:
- //XXX: cast to HdfReader
- break;
-
- case Reader::TXT_READER:
- break;
-
- case Reader::VTK_READER:
- // XXX: cast to VtkReader
- //fieldReader.open(inputfile);
- fieldReader->get_coordinates(&coords, &nno, &nsd);
- fieldReader->get_connectivity(&connect, &nel, &ndofs);
- //fieldReader->get_dofs(location.c_str(), &dofs, &dofs_nno, &dofs_valdim);
- //fieldReader->get_vector_point_data(location.c_str(), &dofs, &dofs_nno, &dofs_valdim);
- fieldReader->get_scalar_point_data(location.c_str(), &dofs, &dofs_nno, &dofs_valdim);
- //fieldReader->get_point_data(location.c_str(), &dofs, &dofs_nno, &dofs_valdim);
- break;
-
- default:
- break;
- }
-
-
-
- field->dim = nsd;
- field->rank = dofs_valdim;
- //field->meshPart = new cigma::VtkUgMeshPart();
- field->meshPart = new cigma::MeshPart();
- field->meshPart->set_coordinates(coords, nno, nsd);
- field->meshPart->set_connectivity(connect, nel, ndofs);
-
-
- // move to set_mesh()
- field->meshPart->set_cell();
- assert(field->meshPart->cell != 0);
-
- // XXX: Create locator only when necessary
- cigma::AnnLocator *locator = new cigma::AnnLocator();
- field->meshPart->set_locator(locator);
- ////
-
-
- //
- switch (field->dim)
- {
- case 2:
- // 2D case
- switch (ndofs)
+ if (meshIO.mesh_path == "")
{
- case 3:
- field->fe = new cigma::FE();
- //field->fe->cell = new cigma::Tri();
- field->fe->quadrature = new cigma::Quadrature();
- //field->fe->quadrature->set_quadrature(tri_qpts, tri_qwts, tri_nno, tri_nsd);
- //field->fe->quadrature->set_globaldim(tri_nsd);
- //field->fe->jxw = new double[tri_nno];
- //field->fe->basis_tab = new double[tri_nno * ndofs];
- //field->fe->basis_jet = new double[tri_nno * ndofs * tri_nsd];
- //field->meshPart->cell = field->fe->cell;
- break;
- case 4:
- field->fe = new cigma::FE();
- //field->fe->cell = new cigma::Quad();
- field->fe->quadrature = new cigma::Quadrature();
- //field->fe->quadrature->set_quadrature(quad_qpts, quad_qwts, quad_nno, quad_nsd);
- //field->fe->quadrature->set_globaldim(quad_nsd);
- //field->fe->jxw = new double[quad_nno];
- //field->fe->basis_tab = new double[quad_nno * ndofs];
- //field->fe->basis_jet = new double[quad_nno * ndofs * quad_nsd];
- //field->meshPart->cell = field->fe->cell;
- break;
+ meshIO.mesh_path = dofs_file;
}
- break;
- case 3:
- // 3D case
- switch (ndofs)
- {
- case 4:
- field->fe = new cigma::FE();
- //field->fe->cell = new cigma::Tet();
- field->fe->quadrature = new cigma::Quadrature();
- //field->fe->quadrature->set_quadrature(tet_qpts, tet_qwts, tet_nno, tet_nsd);
- //field->fe->quadrature->set_globaldim(3);
- //field->fe->jxw = new double[tet_nno];
- //field->fe->basis_tab = new double[tet_nno * ndofs];
- //field->fe->basis_jet = new double[tet_nno * ndofs * tet_nsd];
- //field->meshPart->cell = field->fe->cell;
- break;
- case 8:
- field->fe = new cigma::FE();
- //field->fe->cell = new cigma::Hex();
- field->fe->quadrature = new cigma::Quadrature();
- //field->fe->quadrature->set_quadrature(hex_qpts, hex_qwts, hex_nno, hex_nsd);
- //field->fe->quadrature->set_globaldim(hex_nsd);
- //field->fe->jxw = new double[hex_nno];
- //field->fe->basis_tab = new double[hex_nno * ndofs];
- //field->fe->basis_jet = new double[hex_nno * ndofs * hex_nsd];
- //field->meshPart->cell = field->fe->cell;
- break;
- }
- break;
- }
- field->fe->cell = field->meshPart->cell;
- ////
+ meshIO.load();
+ assert(meshIO.meshPart != 0);
- //XXX: move to field->set_quadrature(...)
- Quadrature *Q = new cigma::Quadrature();
- //load_quadrature(field->meshPart->cell, Q); XXX
- field->fe = new cigma::FE();
- field->fe->set_cell_quadrature(field->meshPart->cell, Q);
+ 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);
- assert(field->fe != 0);
- assert(field->fe->cell != 0);
+ // XXX: instantiate this only when necessary!
+ AnnLocator *locator = new AnnLocator();
+ field->meshPart->set_locator(locator);
+ field->dofHandler = new DofHandler();
+ field->dofHandler->set_data(dofs, dofs_nno, dofs_valdim);
+ }
- field->dofHandler = new cigma::DofHandler();
- field->dofHandler->meshPart = field->meshPart;
- field->dofHandler->set_data(dofs, dofs_nno, dofs_valdim);
-
- // */
return;
}
@@ -818,7 +514,9 @@
void FieldIO::save()
{
+ cout << "Calling FieldIO::save()" << endl;
assert(field != 0);
+ assert(false);
}
Modified: cs/benchmark/cigma/trunk/src/Misc.h
===================================================================
--- cs/benchmark/cigma/trunk/src/Misc.h 2008-02-01 20:52:04 UTC (rev 9210)
+++ cs/benchmark/cigma/trunk/src/Misc.h 2008-02-01 20:52:06 UTC (rev 9211)
@@ -98,7 +98,7 @@
double pick_from_interval(double a, double b);
void bbox_random_point(double minpt[3], double maxpt[3], double x[3]);
-void configure_quadrature(AnyOption *opt, QuadratureIO *quadratureIO);
+void configure_quadrature(AnyOption *opt, QuadratureIO *quadratureIO, const char *opt_prefix);
void configure_mesh(AnyOption *opt, MeshIO *meshIO, const char *opt_prefix);
void configure_field(AnyOption *opt, FieldIO *fieldIO, const char *opt_prefix);
More information about the cig-commits
mailing list