[cig-commits] r13544 - cs/cigma/trunk/src

[luis at geodynamics.org]

Author: luis
Date: 2008-12-09 18:13:52 -0800 (Tue, 09 Dec 2008)
New Revision: 13544

Modified:
   cs/cigma/trunk/src/Quadrature.cpp
   cs/cigma/trunk/src/Quadrature.h
Log:
Factory and setter methods for Quadrature class

Modified: cs/cigma/trunk/src/Quadrature.cpp
===================================================================
--- cs/cigma/trunk/src/Quadrature.cpp	2008-12-10 02:13:50 UTC (rev 13543)
+++ cs/cigma/trunk/src/Quadrature.cpp	2008-12-10 02:13:52 UTC (rev 13544)
@@ -1,6 +1,13 @@
+#include "Quadrature.h"
+#include "io_file_reader.h"
+#include "core_array.h"
+#include "core_readers.h"
 #include <iostream>
 #include <cassert>
-#include "Quadrature.h"
+#include <string>
+#include <sstream>
+
+using namespace std;
 using namespace cigma;
 using boost::shared_ptr;
 
@@ -133,6 +140,46 @@
     return Q;
 }
 
+shared_ptr<Quadrature> Quadrature::New(const QuadratureInfo& q_info)
+{
+    //TRI_LOG_STR("Quadrature::New()");
+    shared_ptr<Quadrature> Q;
+
+    if (q_info.p_quadrature)
+    {
+        //TRI_LOG(q_info.p_quadrature);
+        Q = shared_ptr<Quadrature>(new Quadrature);
+        Q->setPath(q_info.p_quadrature);
+        if (Q->cell_type == Cell::NONE)
+        {
+            Q->cell_type = Cell::string2type(q_info.cell_type_name);
+        }
+        return Q;
+    }
+
+    if (q_info.p_weights && q_info.p_points)
+    {
+        //TRI_LOG(q_info.p_weights);
+        //TRI_LOG(q_info.p_points);
+        Q = shared_ptr<Quadrature>(new Quadrature);
+        Q->setPath2(q_info.p_weights, q_info.p_points);
+        Q->cell_type = Cell::string2type(q_info.cell_type_name);
+        return Q;
+    }
+
+    if (q_info.cell_type_name != "")
+    {
+        //TRI_LOG(q_info.cell_type_name);
+        Cell::type cell_type = Cell::string2type(q_info.cell_type_name);
+        if (cell_type != Cell::NONE)
+        {
+            Q = Quadrature::default_rule(cell_type);
+        }
+    }
+
+    return Q;
+}
+
 // ----------------------------------------------------------------------------
 
 Quadrature::Quadrature()
@@ -230,3 +277,85 @@
 // ----------------------------------------------------------------------------
 
 
+void Quadrature::setPath(const DataPath& p_quad)
+{
+    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);
+
+    if (wx == 0)
+    {
+        std::ostringstream stream;
+        stream << "Failed to read quadrature from path '" << p_quad << "'" << std::ends;
+        throw cigma::Exception("Quadrature::setPath", stream.str());
+    }
+
+    const int npts = wx->n_points();
+    const int ndim = -1 + wx->n_dim();
+
+    this->reinit(npts, ndim);
+
+    for (int q = 0; q < npts; q++)
+    {
+        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;
+}
+
+void Quadrature::setPath2(const DataPath& p_weights, const DataPath& p_points)
+{
+    TRI_LOG_STR("Quadrature::setPath2");
+    TRI_LOG(p_weights);
+    TRI_LOG(p_points);
+
+    // load weights array from p_weights
+    cigma::array<double> *wts = ReadArray(p_weights);
+    if (wts == 0)
+    {
+        std::ostringstream stream;
+        stream << "Failed to read quadrature weights from path '" << p_weights << "'" << std::ends;
+        throw cigma::Exception("Quadrature::setPath2", stream.str());
+    }
+    this->weights = wts->_data;
+    wts->_data = 0;
+    delete wts;
+
+    // load points array from p_points
+    cigma::array<double> *pts = ReadArray(p_points);
+    if (pts == 0)
+    {
+        std::ostringstream stream;
+        stream << "Failed to read quadrature points from path '" << p_points << "'" << std::ends;
+        throw cigma::Exception("Quadrature::setPath2", stream.str());
+    }
+    this->points = pts->_data;
+    pts->_data = 0;
+    delete pts;
+}
+
+void Quadrature::setCell(string cell_name)
+{
+    TRI_LOG_STR("Quadrature::setCell");
+    TRI_LOG(cell_name);
+
+    this->cell_type = Cell::string2type(cell_name);
+
+    if (cell_type != Cell::NONE)
+    {
+        shared_ptr<Quadrature> Q = Quadrature::default_rule(cell_type);
+        this->reinit(Q->npts, Q->ndim);
+        this->setData(Q->points, Q->weights);
+    }
+}
+
+// ----------------------------------------------------------------------------

Modified: cs/cigma/trunk/src/Quadrature.h
===================================================================
--- cs/cigma/trunk/src/Quadrature.h	2008-12-10 02:13:50 UTC (rev 13543)
+++ cs/cigma/trunk/src/Quadrature.h	2008-12-10 02:13:52 UTC (rev 13544)
@@ -1,10 +1,13 @@
 #ifndef __CIGMA_QUADRATURE_H__
 #define __CIGMA_QUADRATURE_H__
 
-#include "Cell.h"
 #include <boost/noncopyable.hpp>
 #include <boost/shared_ptr.hpp>
 
+#include "Cell.h"
+#include "DataPath.h"
+#include "core_args.h"
+
 namespace cigma
 {
     class Quadrature;
@@ -24,17 +27,24 @@
     int n_points() const;
     int n_dim() const;
 
-    Cell::type getCellType() const;
     double getPoint(int q, int j) const;
     double getWeight(int q) const;
+    Cell::type getCellType() const;
 
-    void setCellType(Cell::type cell_type);
+    void reinit(int npts, int ndim);
+    void setData(double *points, double *weights); // XXX: rename to init(double*,double*)
+
     void setPoint(int q, int j, double val);
     void setWeight(int q, double val);
+    void setCellType(Cell::type cell_type);
 
-    void reinit(int npts, int ndim);
-    void setData(double *points, double *weights); // XXX: rename to init(double*,double*)
+    /* incremental setters */
+    void setPath(const DataPath& p_quad);
+    void setPath2(const DataPath& p_weights, const DataPath& p_points);
+    void setCell(std::string cell_name);
 
+    /* static factory methods */
+    static boost::shared_ptr<Quadrature> New(const QuadratureInfo& q_info);
     static boost::shared_ptr<Quadrature> default_rule(Cell::type cell_type);
 
 public: