sdf_domainmask.h

Go to the documentation of this file.

/*
  This file is part of MADNESS.

  Copyright (C) 2007,2010 Oak Ridge National Laboratory

  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 of the License, or
  (at your option) any later version.

  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

  For more information please contact:

  Robert J. Harrison
  Oak Ridge National Laboratory
  One Bethel Valley Road
  P.O. Box 2008, MS-6367

  email: harrisonrj@ornl.gov
  tel:   865-241-3937
  fax:   865-572-0680

  $Id: sdf_shape.h 1792 2010-01-26 12:58:18Z rjharrison $
*/

#ifndef MADNESS_MRA_SDF_DOMAINMASK_H__INCLUDED
#define MADNESS_MRA_SDF_DOMAINMASK_H__INCLUDED

#include <madness/mra/mra.h>

namespace madness {

    template <std::size_t NDIM>
    class SignedDFInterface {
    public:
        virtual double sdf(const Vector<double,NDIM>& x) const = 0;

        virtual Vector<double,NDIM> grad_sdf(const Vector<double,NDIM>& x)
            const = 0;

        virtual ~SignedDFInterface() {}
    };

    class DomainMaskInterface {
    public:
        virtual double mask(double d) const = 0;

        virtual double dmask(double d) const = 0;

        virtual double surface(double d) const = 0;

        virtual double dsurface(double d) const = 0;

        virtual ~DomainMaskInterface() {}
    };

    template <std::size_t NDIM>
    class DomainMaskSDFFunctor : public FunctionFunctorInterface<double,NDIM> {

    private:
        DomainMaskSDFFunctor() {}

    private: // protected may be better if this becomes highly inherited
        std::shared_ptr<DomainMaskInterface> mask;

        std::shared_ptr<SignedDFInterface<NDIM> > sdf;

        int mswitch; 

    public:
        // switch values
        static const int MASK; 
        static const int MASK_COMPLEMENT; 
        static const int DMASK; 
        static const int SURFACE; 
        static const int DSURFACE; 

        DomainMaskSDFFunctor(std::shared_ptr<DomainMaskInterface> mask,
                std::shared_ptr<SignedDFInterface<NDIM> > sdf)
            : mask(mask), sdf(sdf), mswitch(MASK)
        {}

        DomainMaskSDFFunctor(std::shared_ptr<DomainMaskInterface> mask,
                std::shared_ptr<SignedDFInterface<NDIM> > sdf, int _mswitch)
            : mask(mask), sdf(sdf), mswitch(MASK) {
            if(_mswitch == MASK || _mswitch == DMASK || _mswitch == SURFACE ||
               _mswitch == DSURFACE || _mswitch == MASK_COMPLEMENT) {
                mswitch = _mswitch;
            }
            else {
                error("Unrecognized function option in DomainMaskSDFFunctor" \
                      "::DomainMaskSDFFunctor()");
            }
        }

        double operator()(const Vector<double,NDIM>& x) const {
            if(mswitch == DomainMaskSDFFunctor<NDIM>::MASK)
                return mask->mask(sdf->sdf(x));
            else if(mswitch == DomainMaskSDFFunctor<NDIM>::MASK_COMPLEMENT)
                return 1.0 - mask->mask(sdf->sdf(x));
            else if(mswitch == DomainMaskSDFFunctor<NDIM>::DMASK)
                return mask->dmask(sdf->sdf(x));
            else if(mswitch == DomainMaskSDFFunctor<NDIM>::SURFACE)
                return mask->surface(sdf->sdf(x));
            else if(mswitch == DomainMaskSDFFunctor<NDIM>::DSURFACE)
                return mask->dsurface(sdf->sdf(x));
            else {
                error("Unknown function from DomainMaskInterface in " \
                      "DomainMaskSDFFunctor::operator()");
                return 0.0;
            }
        }

        void setMaskFunction(int _mswitch) {
            if(_mswitch == MASK || _mswitch == DMASK || _mswitch == SURFACE ||
               _mswitch == DSURFACE || _mswitch == MASK_COMPLEMENT) {
                mswitch = _mswitch;
            }
            else {
                error("Unrecognized function option in DomainMaskSDFFunctor" \
                      "::setMaskFunction()");
            }
        }

        virtual ~DomainMaskSDFFunctor() {}
    };

    template<std::size_t NDIM>
    const int DomainMaskSDFFunctor<NDIM>::MASK = 1;

    template<std::size_t NDIM>
    const int DomainMaskSDFFunctor<NDIM>::MASK_COMPLEMENT = 2;

    template<std::size_t NDIM>
    const int DomainMaskSDFFunctor<NDIM>::DMASK = 3;

    template<std::size_t NDIM>
    const int DomainMaskSDFFunctor<NDIM>::SURFACE = 4;

    template<std::size_t NDIM>
    const int DomainMaskSDFFunctor<NDIM>::DSURFACE = 5;

    class LLRVDomainMask : public DomainMaskInterface {
    private:
        LLRVDomainMask() : epsilon(0.0) {} 

    protected:
        const double epsilon; 

    public:
        LLRVDomainMask(double epsilon)
            : epsilon(epsilon)
        {}

        double mask(double d) const {
            if (d > 8.0*epsilon) {
                return 0.0; // we're safely outside
            }
            else if (d < -8.0*epsilon) {
                return 1.0; // inside
            }
            else {
                return 0.5 * (1.0 - tanh(3.0 * d / epsilon));
            }
        }

        double dmask(double d) const {
            if (fabs(d) > 8.0*epsilon) {
                return 0.0; // we're safely outside or inside
            }
            else {
                double tanh3d = tanh(3.0*d/epsilon);
                return 1.5*(tanh3d*tanh3d - 1.0) / epsilon;
            }
        }

        double surface(double d) const {
            double phi = mask(d);
            double phic = 1.0 - phi;
            return 36.0*phi*phi*phic*phic/epsilon;
        }

        double dsurface(double d) const {
            double phi = mask(d);
            double dphi = dmask(d);
            return 72.0*phi*(1.0-phi)*dphi*(1.0 - 2.0*phi)/epsilon;
        }

        virtual ~LLRVDomainMask() {}
    };

    class GaussianDomainMask : public DomainMaskInterface {
    private:
        GaussianDomainMask() : epsilon(0.0) {} 

    protected:
        const double epsilon; 

    public:
        GaussianDomainMask(double epsilon)
            : epsilon(epsilon)
        {}

        double mask(double d) const {
            if (d > 8.0*epsilon) {
                return 0.0; // we're safely outside
            }
            else if (d < -8.0*epsilon) {
                return 1.0; // inside
            }
            else {
                return (1.0 - erf(d / (sqrt(2.0) * epsilon))) * 0.5;
            }
        }

        double dmask(double d) const {
            if (fabs(d) > 8.0*epsilon) {
                return 0.0; // we're safely outside or inside
            }
            else {
                return -exp(-d*d*0.5/(epsilon*epsilon)) / (sqrt(2.0*constants::pi)
                    * epsilon);
            }
        }

        double surface(double d) const {
            return exp(-d*d*0.5/(epsilon*epsilon)) / (sqrt(2.0*constants::pi)
                * epsilon);
        }

        double dsurface(double d) const {
            return -exp(-d*d*0.5/(epsilon*epsilon)) * d / (sqrt(2.0*constants::pi)
                * epsilon*epsilon*epsilon);
        }

        virtual ~GaussianDomainMask() {}
    };

} // end of madness namespace

#endif // MADNESS_MRA_SDF_DOMAINMASK_H__INCLUDED