[cig-commits] r14088 - cs/cigma/trunk/src
luis at geodynamics.org
luis at geodynamics.org
Wed Feb 18 08:14:45 PST 2009
Author: luis
Date: 2009-02-18 08:14:44 -0800 (Wed, 18 Feb 2009)
New Revision: 14088
Added:
cs/cigma/trunk/src/cli_element_info_cmd.cpp
cs/cigma/trunk/src/cli_element_info_cmd.h
cs/cigma/trunk/src/cli_function_info_cmd.cpp
cs/cigma/trunk/src/cli_function_info_cmd.h
Modified:
cs/cigma/trunk/src/cli_application.cpp
Log:
Added element-info and function-info commands
Modified: cs/cigma/trunk/src/cli_application.cpp
===================================================================
--- cs/cigma/trunk/src/cli_application.cpp 2009-02-18 16:14:43 UTC (rev 14087)
+++ cs/cigma/trunk/src/cli_application.cpp 2009-02-18 16:14:44 UTC (rev 14088)
@@ -18,8 +18,9 @@
#include "cli_eval_cmd.h"
#include "cli_extract_cmd.h"
#include "cli_list_cmd.h"
-#include "cli_info_cmd.h"
#include "cli_mesh_info_cmd.h"
+#include "cli_function_info_cmd.h"
+#include "cli_element_info_cmd.h"
using namespace std;
using namespace cigma;
@@ -33,8 +34,9 @@
this->addCommand(new EvalCmd());
this->addCommand(new ExtractCmd());
this->addCommand(new ListCmd());
- this->addCommand(new InfoCmd());
this->addCommand(new MeshInfoCmd());
+ this->addCommand(new FunctionInfoCmd());
+ this->addCommand(new ElementInfoCmd());
}
@@ -211,13 +213,14 @@
* more elegant, since we only ever add new commands at compile time.
*/
string prefix = " ";
+ const int max_name_width = 13; // XXX: calculate max_width from list of command names
for (CmdNames::iterator i = names.begin(); i != names.end(); ++i)
{
Command *cmd = getCommand(*i);
if (cmd != 0)
{
cout << prefix
- << std::setw(12) << std::left
+ << std::setw(max_name_width + 2) << std::left
<< cmd->name
<< cmd->brief
<< endl;
Added: cs/cigma/trunk/src/cli_element_info_cmd.cpp
===================================================================
--- cs/cigma/trunk/src/cli_element_info_cmd.cpp (rev 0)
+++ cs/cigma/trunk/src/cli_element_info_cmd.cpp 2009-02-18 16:14:44 UTC (rev 14088)
@@ -0,0 +1,280 @@
+#include "cli_element_info_cmd.h"
+#include "Common.h"
+#include <iostream>
+#include <iomanip>
+#include "fe_hex8.h"
+#include "fe_tet4.h"
+#include "fe_tet10.h"
+#include "fe_quad4.h"
+#include "fe_tri3.h"
+#include "fe_tri6.h"
+#include "Quadrature.h"
+
+using namespace cigma;
+using namespace std;
+
+using boost::shared_ptr;
+
+namespace po = boost::program_options;
+
+ElementInfoCmd::ElementInfoCmd()
+{
+ name = "element-info";
+ brief = "Show information about available element types";
+ usage = "Usage: cigma element-info [ <CellType> ]";
+ appendix = "";
+}
+
+ElementInfoCmd::~ElementInfoCmd()
+{
+}
+
+void ElementInfoCmd::defineOptions()
+{
+ BaseCmd::defineOptions();
+ opts_section = "Options: ";
+ opts.add_options()
+ ("cell-type,c", po::value<string>(&celltype), "Element cell type")
+ ;
+ args.add("cell-type", -1);
+}
+
+void ElementInfoCmd::configure()
+{
+ BaseCmd::configure();
+ TRI_LOG_STR("ElementInfoCmd::configure()");
+}
+
+static void print_header(const std::string& celltype)
+{
+ cout << "Information on cell '" << celltype << "'" << endl;
+ cout << endl;
+}
+
+static void print_cell_nodes(double *refverts, int nno, int ndim)
+{
+ const char *indent = " ";
+ cout << "Reference-Cell Nodes:" << endl;
+ for (int i = 0; i < nno; i++)
+ {
+ cout << indent;
+ cout << noshowpos << i << showpos << fixed;
+ cout << " ";
+ for (int j = 0; j < ndim; j++)
+ {
+ cout << refverts[ndim*i + j] << " ";
+ }
+ cout << endl;
+ }
+ cout << noshowpos;
+ cout << endl;
+}
+
+static void print_hex8_shape()
+{
+ const char *indent = " ";
+ cout << "Shape functions:" << endl;
+ cout << indent << "N[0] = 0.125 * (1.0 - u) * (1.0 - v) * (1.0 - w)\n";
+ cout << indent << "N[1] = 0.125 * (1.0 + u) * (1.0 - v) * (1.0 - w)\n";
+ cout << indent << "N[2] = 0.125 * (1.0 + u) * (1.0 + v) * (1.0 - w)\n";
+ cout << indent << "N[3] = 0.125 * (1.0 - u) * (1.0 + v) * (1.0 - w)\n";
+ cout << indent << "N[4] = 0.125 * (1.0 - u) * (1.0 - v) * (1.0 + w)\n";
+ cout << indent << "N[5] = 0.125 * (1.0 + u) * (1.0 - v) * (1.0 + w)\n";
+ cout << indent << "N[6] = 0.125 * (1.0 + u) * (1.0 + v) * (1.0 + w)\n";
+ cout << indent << "N[7] = 0.125 * (1.0 - u) * (1.0 + v) * (1.0 + w)\n";
+ cout << endl;
+}
+
+static void print_tet4_shape()
+{
+ const char *indent = " ";
+ cout << "Shape functions:" << endl;
+ cout << indent << "N[0] = 1.0 - u - v - w\n";
+ cout << indent << "N[1] = u \n";
+ cout << indent << "N[2] = v \n";
+ cout << indent << "N[3] = w\n";
+ cout << endl;
+}
+
+static void print_tet10_shape()
+{
+ const char *indent = " ";
+ cout << "Shape functions:" << endl;
+ cout << indent << "N[0] = 1 - 3*u - 3*v + 4*u*v + 4*u*w + 4*v*w + 2*u*u + 2*v*v + 2*w*w\n";
+ cout << indent << "N[1] = 2*u*u - u\n";
+ cout << indent << "N[2] = 2*v*v - v\n";
+ cout << indent << "N[3] = 2*w*w - w\n";
+ cout << indent << "N[4] = 4*u*v\n";
+ cout << indent << "N[5] = 4*v*w\n";
+ cout << indent << "N[6] = 4*u*w\n";
+ cout << indent << "N[7] = 4*(w - u*w - v*w - w*w)\n";
+ cout << indent << "N[8] = 4*(v - u*v - v*w - v*v)\n";
+ cout << indent << "N[9] = 4*(u - u*v - u*w - u*u)\n";
+ cout << endl;
+}
+
+static void print_quad4_shape()
+{
+ const char *indent = " ";
+ cout << "Shape functions:" << endl;
+ cout << indent << "N[0] = 0.25 * (1.0 - u) * (1.0 - v)\n";
+ cout << indent << "N[1] = 0.25 * (1.0 + u) * (1.0 - v)\n";
+ cout << indent << "N[2] = 0.25 * (1.0 + u) * (1.0 + v)\n";
+ cout << indent << "N[3] = 0.25 * (1.0 - u) * (1.0 + v)\n";
+ cout << endl;
+}
+
+static void print_tri3_shape()
+{
+ const char *indent = " ";
+ cout << "Shape functions:" << endl;
+ cout << indent << "N[0] = 1.0 - u - v\n";
+ cout << indent << "N[1] = u \n";
+ cout << indent << "N[2] = v\n";
+ cout << endl;
+}
+
+static void print_tri6_shape()
+{
+ const char *indent = " ";
+ cout << "Shape functions:" << endl;
+ cout << indent << "N[0] = 1 - 3*u - 3*v + 4*u*v + 2*u*u + 2*v*v\n";
+ cout << indent << "N[1] = 2*u*u - u\n";
+ cout << indent << "N[2] = 2*v*v - v\n";
+ cout << indent << "N[3] = 4*u*v\n";
+ cout << indent << "N[4] = 4*(v - u*v - v*v)\n";
+ cout << indent << "N[5] = 4*(u - u*v - u*u)\n";
+ cout << endl;
+}
+
+static void print_default_quadrature(Cell::type celltype)
+{
+ const char *indent = " ";
+ cout << "Default integration rule (weights & points):" << endl;
+ shared_ptr<Quadrature> Q = Quadrature::default_rule(celltype);
+ if (Q)
+ {
+ cout << noshowpos << fixed;
+ for (int q = 0; q < Q->n_points(); q++)
+ {
+ cout << indent;
+ cout << noshowpos << Q->getWeight(q);
+ cout << " " << showpos;
+ for (int j = 0; j < Q->n_dim(); j++)
+ {
+ cout << Q->getPoint(q, j) << " ";
+ }
+ cout << endl;
+ }
+ }
+ else
+ {
+ cout << indent << "No rule found for cell '" << Cell::type2string(celltype) << endl;
+ }
+ cout << endl;
+}
+
+int ElementInfoCmd::run()
+{
+ BaseCmd::run();
+ TRI_LOG_STR("ElementInfoCmd::run()");
+
+ const char *indent = " ";
+
+ if (celltype == "")
+ {
+ TRI_LOG_STR(">> Listing available elements");
+
+ //
+ // XXX: define static method Cell::listElements()
+ //
+ cout << "List of elements available for use in a Field:" << endl;
+ cout << indent << "tet4, tet10" << endl;
+ cout << indent << "hex8" << endl;
+ cout << indent << "tri3, tri6" << endl;
+ cout << indent << "quad4" << endl;
+
+ //
+ // XXX: restore full list one by one.
+ //
+ //cout << indent << "hex8, hex20, hex27" << endl;
+ //cout << indent << "tet4, tet10" << endl;
+ //cout << indent << "prism6, prism15, prism18" << endl;
+ //cout << indent << "pyramid5" << endl;
+ //cout << indent << "quad4, quad8, quad9" << endl;
+ //cout << indent << "tri3, tri6" << endl;
+ //cout << indent << "edge2, edge3, edge4" << endl;
+ }
+ else
+ {
+ TRI_LOG_STR(">> Showing information on cell type");
+ TRI_LOG(celltype);
+
+ // print name, dimension, number of nodes
+ // print reference nodes and coordinates
+ // print shape functions
+ // print default quadrature pts & wts
+
+ // XXX: Since the list of elements is small, we hardcode information
+ // about them here. As soon as we start adding more element types,
+ // this information should be moved to its corresponding Cell
+ // subclass.
+
+ if (celltype == "hex8")
+ {
+ hex8 cell;
+ print_header(celltype);
+ print_cell_nodes(hex8::refverts, cell.n_nodes(), cell.n_celldim());
+ print_hex8_shape();
+ print_default_quadrature(cell.cell_type());
+ }
+ else if (celltype == "tet4")
+ {
+ tet4 cell;
+ print_header(celltype);
+ print_cell_nodes(tet4::refverts, cell.n_nodes(), cell.n_celldim());
+ print_tet4_shape();
+ print_default_quadrature(cell.cell_type());
+ }
+ else if (celltype == "tet10")
+ {
+ tet10 cell;
+ print_header(celltype);
+ print_cell_nodes(tet10::refverts, cell.n_nodes(), cell.n_celldim());
+ print_tet10_shape();
+ print_default_quadrature(cell.cell_type());
+ }
+ else if (celltype == "quad4")
+ {
+ quad4 cell;
+ print_header(celltype);
+ print_cell_nodes(quad4::refverts, cell.n_nodes(), cell.n_celldim());
+ print_quad4_shape();
+ print_default_quadrature(cell.cell_type());
+ }
+ else if (celltype == "tri3")
+ {
+ tri3 cell;
+ print_header(celltype);
+ print_cell_nodes(tri3::refverts, cell.n_nodes(), cell.n_celldim());
+ print_tri3_shape();
+ print_default_quadrature(cell.cell_type());
+ }
+ else if (celltype == "tri6")
+ {
+ tri6 cell;
+ print_header(celltype);
+ print_cell_nodes(tri6::refverts, cell.n_nodes(), cell.n_celldim());
+ print_tri6_shape();
+ print_default_quadrature(cell.cell_type());
+ }
+ else
+ {
+ cout << "No information available on '" << celltype << "'" << endl;
+ return 1;
+ }
+ }
+ return 0;
+}
+
+
Added: cs/cigma/trunk/src/cli_element_info_cmd.h
===================================================================
--- cs/cigma/trunk/src/cli_element_info_cmd.h (rev 0)
+++ cs/cigma/trunk/src/cli_element_info_cmd.h 2009-02-18 16:14:44 UTC (rev 14088)
@@ -0,0 +1,27 @@
+#ifndef CIGMA_ELEMENT_INFO_CMD_H
+#define CIGMA_ELEMENT_INFO_CMD_H
+
+#include "cli_base_cmd.h"
+#include "Cell.h"
+
+namespace cigma
+{
+ class ElementInfoCmd;
+}
+
+class cigma::ElementInfoCmd : public BaseCmd
+{
+public:
+ ElementInfoCmd();
+ ~ElementInfoCmd();
+
+ void defineOptions();
+ void configure();
+ int run();
+
+public:
+ std::string celltype;
+
+};
+
+#endif
Added: cs/cigma/trunk/src/cli_function_info_cmd.cpp
===================================================================
--- cs/cigma/trunk/src/cli_function_info_cmd.cpp (rev 0)
+++ cs/cigma/trunk/src/cli_function_info_cmd.cpp 2009-02-18 16:14:44 UTC (rev 14088)
@@ -0,0 +1,363 @@
+#include "cli_function_info_cmd.h"
+#include "Common.h"
+#include "Function.h"
+#include "Field.h"
+#include "core_args.h"
+#include "core_info.h"
+#include "core_writers.h"
+#include <iostream>
+#include <iomanip>
+#include <sstream>
+
+using namespace cigma;
+using namespace std;
+using namespace boost;
+
+namespace po = boost::program_options;
+
+const string indent = " ";
+
+// ----------------------------------------------------------------------------
+
+FunctionInfoCmd::FunctionInfoCmd()
+{
+ name = "function-info";
+ brief = "Show list of pre-defined functions";
+ usage = "Usage: cigma function-info";
+ //usage = "Usage: cigma function-info <FunctionObject>";
+ appendix = "";
+}
+
+FunctionInfoCmd::~FunctionInfoCmd()
+{
+}
+
+void FunctionInfoCmd::defineOptions()
+{
+ BaseCmd::defineOptions();
+ opts_section = "Options: ";
+ opts.add_options()
+ ("function,f", po::value<string>(&fn_name), "Input function")
+ ("mesh,m", po::value<string>(&meshfile), "Mesh associated with function, if present")
+ ("query-point,p", po::value<string>(&pointstr), "Evaluate function at given point")
+ ;
+ args.add("function", -1);
+}
+
+void FunctionInfoCmd::configure()
+{
+ BaseCmd::configure();
+ TRI_LOG_STR("FunctionInfoCmd::configure()");
+
+ if (vm.count("function"))
+ {
+ FunctionInfo fn_info;
+ fn_info.fn_name = fn_name;
+ fn_info.field_info.p_field = DataPath(fn_name);
+
+ if (vm.count("mesh"))
+ {
+ fn_info.field_info.mesh_info.p_mesh = DataPath(meshfile);
+ }
+
+ fn = Function::NewFunction(fn_info);
+ }
+}
+
+static bool is_registered_function(const string& fn_name)
+{
+ return (fn_name == "zero")
+ || (fn_name == "one")
+ || (fn_name == "test.cube")
+ || (fn_name == "test.square")
+ || (fn_name == "bm.ssnog.displacement")
+ || (fn_name == "bm.gale.circular_inclusion.pressure")
+ ;
+}
+static void print_registered_function_info(const string& fn_name)
+{
+ // XXX: Since the list of initially available functions is
+ // currently small, we hardcode the information here.
+ // As the list grows, this kind of information should be moved
+ // to the constructor of the respective Function subclasses,
+ // and queried here by using an API.
+
+ if (fn_name == "zero")
+ {
+ cout << "Function Info:" << endl;
+ cout << indent << "Name = " << fn_name << endl;
+ }
+ else if (fn_name == "one")
+ {
+ cout << "Function Info:" << endl;
+ cout << indent << "Name = " << fn_name << endl;
+ }
+ else if (fn_name == "test.cube")
+ {
+ cout << "Function Info:" << endl;
+ cout << indent << "Name = " << fn_name << endl;
+ }
+ else if (fn_name == "test.square")
+ {
+ cout << "Function Info:" << endl;
+ cout << indent << "Name = " << fn_name << endl;
+ }
+ else if (fn_name == "bm.ssnog.displacement")
+ {
+ cout << "Function Info:" << endl;
+ cout << indent << "Name = " << fn_name << endl;
+ }
+ else if (fn_name == "bm.gale.circular_inclusion.pressure")
+ {
+ cout << "Function Info:" << endl;
+ cout << indent << "Name = " << fn_name << endl;
+ }
+}
+
+
+int FunctionInfoCmd::run()
+{
+ BaseCmd::run();
+
+ if (fn_name == "")
+ {
+ TRI_LOG_STR(">> Listing available functions");
+ cout << "List of registered functions:" << endl;
+ global_function_registry.listFunctions();
+ return 0;
+ }
+
+ TRI_LOG_STR(">> Showing information on function");
+ // print name, short description
+ // print domain description (including dimension)
+ // print description of function values (including rank dim)
+
+ bool registered = is_registered_function(fn_name);
+
+ if (!fn)
+ {
+ cout << "No information available on '" << fn_name << "'" << endl;
+ return 1;
+ }
+
+ int dim, rank;
+ dim = fn->n_dim();
+ rank = fn->n_rank();
+
+ TRI_LOG(dim);
+ TRI_LOG(rank);
+
+ if (registered)
+ {
+ print_registered_function_info(fn_name);
+ }
+
+ Field *field = 0;
+ if (dynamic_cast<Field*>(&(*fn)))
+ {
+ field = static_cast<Field*>(&(*fn));
+ }
+
+ cout << setprecision(8);
+
+ if (vm.count("query-point"))
+ {
+ TRI_LOG_STR(">> Running point query");
+
+ //
+ // Handle case when query-point is specified
+ //
+ valarray<double> point;
+ get_point_from_string(pointstr, point);
+
+ if ((fn_name == "zero") || (fn_name == "one"))
+ {
+ dim = point.size();
+ }
+
+ if (!registered)
+ {
+ cout << "Function Info:" << endl;
+ cout << indent << "Name = " << fn_name << endl;
+ }
+
+ cout << indent << "Dimension of domain = " << dim << endl;
+ cout << indent << "Dimension of range = " << rank << endl;
+
+ if (field != 0)
+ {
+ // Mesh Info
+ print_mesh_info(cout, *(field->mesh), meshfile);
+ }
+
+ cout << endl;
+ cout << "Point Query:" << endl;
+
+ if (point.size() != dim)
+ {
+ cout << "Error: Dimension of point (" << point.size()
+ << ") does not match dimension of function domain."
+ << endl;
+ return 1;
+ }
+
+ /*
+ int cellid = -1;
+ bool found = false;
+ if (field != 0)
+ {
+ found = field->mesh->findCell(&point[0], &cellid);
+ if (found)
+ {
+ cout << indent << "Found point in cell " << cellid << endl;
+ }
+ }
+ // */
+
+ valarray<double> value;
+ value.resize(rank);
+ bool found2 = fn->eval(&point[0], &value[0]);
+
+ // Print point
+ cout << indent << "Point ";
+ print_point(cout, point);
+ cout << endl;
+
+ // Print function value
+ cout << indent << "Value ";
+ if (found2)
+ {
+ print_point(cout, value);
+ }
+ else
+ {
+ cout << "undefined";
+ }
+ cout << endl;
+
+
+ // Print node coordinates in cell, as well as associated degrees of freedom
+ if (field != 0)
+ {
+ int ndim = field->mesh->n_dim();
+ int ncells = field->mesh->n_cells();
+ assert(dim == ndim);
+
+ int cellid = -1;
+ bool found = field->mesh->findCell(&point[0], &cellid);
+
+ if ((0 <= cellid) && (cellid < ncells))
+ {
+ int i,j;
+ shared_ptr<Cell> cell = field->mesh->getCell();
+ field->mesh->getCell(cellid, *cell);
+ const int nno_per_cell = cell->n_nodes();
+ const int celldim = cell->n_celldim();
+ const double *globverts = cell->globverts;
+
+ valarray<double> local_dofs;
+ local_dofs.resize(nno_per_cell * rank);
+ field->getCellDofs(cellid, &local_dofs[0]);
+
+ valarray<double> refpoint;
+ //refpoint.resize(celldim);
+ refpoint.resize(3);
+ cell->xyz2uvw(&point[0], &refpoint[0]);
+
+ valarray<double> N;
+ N.resize(nno_per_cell);
+ cell->shape(1, &refpoint[0], &N[0]);
+
+ double jac[3][3];
+ const double jacobian = cell->jacobian(&refpoint[0], jac);
+ const double volume = cell->volume();
+
+ //
+ // Found point in cell <CellID>
+ //
+ cout << endl;
+ cout << "Found point in cell " << cellid << endl;
+
+ cout << indent << "Global Point = ";
+ print_point(cout, point);
+ cout << endl;
+
+ cout << indent << "Reference Point = ";
+ print_point(cout, refpoint, celldim);
+ cout << endl;
+
+ cout << indent << "Cell Jacobian = " << jacobian << endl;
+ cout << indent << "Cell Volume = " << volume << endl;
+
+ //
+ // Node coordinates for cell <CellID>
+ //
+ print_node_coords_of_cell(cout, *(field->mesh), cellid);
+
+ //
+ // Degrees of freedom for cell <CellID>
+ //
+ cout << endl;
+ cout << "Degrees of freedom for cell " << cellid << endl;
+ for (i = 0; i < nno_per_cell; i++)
+ {
+ int n = field->mesh->connect->getId(cellid, i);
+
+ cout << " ";
+ cout << noshowpos;
+ cout << setw(7) << right << n;
+ cout << showpos;
+ cout << " ";
+
+ for (j = 0; j < rank; j++)
+ {
+ cout << local_dofs[rank*i + j] << " ";
+ }
+
+ cout << endl;
+ }
+
+ cout << noshowpos;
+
+ //
+ // Shape function values at <RefPoint>
+ //
+ cout << endl;
+ cout << "Shape function values at reference point" << endl;
+ for (i = 0; i < nno_per_cell; i++)
+ {
+ cout << indent << noshowpos
+ << "N[" << i << "] = " << showpos
+ << N[i]
+ << endl;
+ }
+
+ }
+ else
+ {
+ cout << "Error: Point not found in any cell!" << endl;
+ }
+
+ } // if (field != 0)
+
+ }
+ else
+ {
+ // Case when no query-point is given
+ if (!registered)
+ {
+ cout << "Function Info:" << endl;
+ cout << indent << "Name = " << fn_name << endl;
+ cout << indent << "Dimension of domain = " << dim << endl;
+ cout << indent << "Dimension of range = " << rank << endl;
+
+ if (field != 0)
+ {
+ print_mesh_info(cout, *(field->mesh), meshfile);
+ }
+ }
+ }
+
+ return 0;
+}
+
Added: cs/cigma/trunk/src/cli_function_info_cmd.h
===================================================================
--- cs/cigma/trunk/src/cli_function_info_cmd.h (rev 0)
+++ cs/cigma/trunk/src/cli_function_info_cmd.h 2009-02-18 16:14:44 UTC (rev 14088)
@@ -0,0 +1,28 @@
+#ifndef CIGMA_FUNCTION_INFO_H
+#define CIGMA_FUNCTION_INFO_H
+
+#include "cli_base_cmd.h"
+#include "FunctionRegistry.h"
+
+namespace cigma
+{
+ class FunctionInfoCmd;
+}
+
+class cigma::FunctionInfoCmd : public BaseCmd
+{
+public:
+ FunctionInfoCmd();
+ ~FunctionInfoCmd();
+
+ void defineOptions();
+ void configure();
+ int run();
+
+public:
+ std::string fn_name;
+ std::string meshfile;
+ std::string pointstr;
+ boost::shared_ptr<Function> fn;
+};
+#endif
More information about the CIG-COMMITS
mailing list