scalars.hpp

Go to the documentation of this file.

#ifndef __FDL_SCALARS_HPP
#define __FDL_SCALARS_HPP
#include <cmath>
#include <cassert>
#include "ftnsbase.hpp"
#include "io/xmlparser.h"

namespace fdl {
        
template <typename scalar_t>
class Gaussian : public ScalarBase<scalar_t>
{
private:
        scalar_t alpha_, sigma2_;
        scalar_t x0_[3]; // no reason for phyics in dim >= 4. 
        
public:
        Gaussian(unsigned dim, scalar_t x0[], 
                                scalar_t sigma =scalar_t(1), scalar_t alpha = scalar_t(1)) 
        : ScalarBase<scalar_t>(dim), alpha_(alpha)
        {
                assert (dim <= 3);
                for (unsigned j = 0; j < dim; j++) 
                        x0_[j] = x0[j];
                sigma2_ = sigma * sigma;
        }
        
        Gaussian(TiXmlNode* node)
        : ScalarBase<scalar_t>(3), alpha_(1), sigma2_(1)
        {
                XmlParser Params(node);
                scalar_t sigma;
                x0_[0] = x0_[1] = x0_[2] = scalar_t(0);
                const char* coordNames[3] = {"x0", "y0", "z0"};
                Params.getAttribArray("Parameters", 3, coordNames, x0_);
                Params.getAttrib("Parameters", "alpha", alpha_);
                if (Params.getAttrib("Parameters", "sigma", sigma))
                        sigma2_ = sigma * sigma;
        }
        
        ~Gaussian() {}
        
        virtual scalar_t operator()(scalar_t t, const scalar_t x[])
        {
                scalar_t r2(0), s;
                for (unsigned j = 0; j < this->dim_; j++) {
                        s = x[j] - x0_[j];
                        r2 += s * s;
                }
                return alpha_ * exp( -r2 / sigma2_); 
        }
};

template <typename scalar_t>
class Spherical : public ScalarBase<scalar_t>
{
private:
        scalar_t alpha_, r2_;
        scalar_t x0_[3]; // no reason for phyics in dim >= 4. 
        
public:
        Spherical(unsigned dim, scalar_t x0[], 
                                scalar_t sigma =scalar_t(1), scalar_t alpha = scalar_t(1)) 
        : ScalarBase<scalar_t>(dim), alpha_(alpha) 
        {
                assert (dim <= 3);
                for (unsigned j = 0; j < dim; j++) 
                        x0_[j] = x0[j];
                r2_ = sigma * sigma;
        }

        Spherical(TiXmlNode* node)
        : ScalarBase<scalar_t>(3), alpha_(1), r2_(1)
        {
                XmlParser Params(node);
                scalar_t r;
                x0_[0] = x0_[1] = x0_[2] = scalar_t(0);
                const char* coordNames[3] = {"x0", "y0", "z0"};
                Params.getAttribArray("Parameters", 3, coordNames, x0_);
                Params.getAttrib("Parameters", "alpha", alpha_);
                if (Params.getAttrib("Parameters", "radius", r)) 
                        r2_ = r * r;
        }
        
        ~Spherical() {}
        
        virtual scalar_t operator()(scalar_t t, const scalar_t x[])
        {
                scalar_t r2(0), s;
                for (unsigned j = 0; j < this->dim_; j++) {
                        s = x[j] - x0_[j];
                        r2 += s * s;
                }
                return (r2 < r2_) ? alpha_ : scalar_t(0); 
        }
};

template <typename scalar_t>
class Cylindrical : public ScalarBase<scalar_t>
{
private:
        scalar_t alpha_;
        scalar_t z0_; 
        scalar_t z1_; 
        scalar_t xy0_[2]; 
        scalar_t r2_; 
        
public:
        Cylindrical(unsigned dim, scalar_t xy0[], scalar_t radius, 
                                                        scalar_t z0, scalar_t z1, scalar_t alpha = scalar_t(1))
        : ScalarBase<scalar_t>(dim), z0_(z0), z1_(z1), alpha_(alpha)
        {
                assert (dim <= 3);
                xy0_[0] = xy0[0]; xy0_[1] = xy0[1];
                r2_ = radius * radius;
        }
        
        Cylindrical(TiXmlNode* node)
        : ScalarBase<scalar_t>(3), z0_(0), z1_(0), alpha_(1), r2_(0)
        {
                XmlParser Params(node);
                scalar_t r;
                xy0_[0] = xy0_[1] = scalar_t(0);
                const char* coordNames[2] = {"x0", "y0"};
                Params.getAttribArray("Parameters", 2, coordNames, xy0_);
                Params.getAttrib("Parameters", "zmin", z0_);
                Params.getAttrib("Parameters", "zmax", z1_);
                Params.getAttrib("Parameters", "alpha", alpha_);
                if (Params.getAttrib("Parameters", "radius", r))
                        r2_ = r * r;
        }
        
        ~Cylindrical() {}
        
        virtual scalar_t operator()(scalar_t t, const scalar_t x[])
        {
                scalar_t s0 = x[0] - xy0_[0];
                scalar_t s1 = x[1] - xy0_[1];
                scalar_t r2 = s0 * s0 + s1 * s1;
                
                return (r2 < r2_ && z0_ <= x[2] && x[2] <= z1_) ? alpha_ : scalar_t(0);
        }
};

template <typename scalar_t>
static ScalarBase<scalar_t>* ScalarCatalog(TiXmlNode* node)
{
        std::string defn;
        std::stringstream errBffr;
        XmlParser Data(node);
        if (!Data.getField("Definition", defn)) {
                errBffr << Data.getRootName() << " missing <Definition> tag.\n";
                throw std::invalid_argument(errBffr.str());
        }
        if (!Data.hasChild("Parameters")) {
                errBffr << Data.getRootName() << " missing <Parameters> tag.\n";
                throw std::invalid_argument(errBffr.str());
        }
        if (defn == "Gaussian") {
                return new Gaussian<scalar_t>(node);
        }
        else if (defn == "Spherical") {
                return new Spherical<scalar_t>(node);
        }
        else if (defn == "Cylindrical") {
                return new Cylindrical<scalar_t>(node);
        }
        else {
                errBffr << defn << " not found in fdl::ScalarCatalog.\n";
                throw std::invalid_argument(errBffr.str());
        }
        //std::cout << "ScalarCatalog\t scalar=" << scalar << std::endl;
}
        
} // namespace fdl

#endif // __FDL_SCALARS_HPP