sdf_shape_3D.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$
*/

#ifndef MADNESS_MRA_SDF_SHAPE_3D_H__INCLUDED
#define MADNESS_MRA_SDF_SHAPE_3D_H__INCLUDED

#include <madness/mra/sdf_domainmask.h>

namespace madness {

    class SDFPlane : public SignedDFInterface<3> {
    protected:
        const coord_3d normal; 
        const coord_3d point; 

    public:
        
        SDFPlane(const coord_3d& normal, const coord_3d& point)
            : normal(normal*(1.0/sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2])))
            , point(point) 
        {}

        double sdf(const coord_3d& pt) const {
            return (pt[0]-point[0])*normal[0] + (pt[1]-point[1])*normal[1] + (pt[2]-point[2])*normal[2];
        }

        coord_3d grad_sdf(const coord_3d& pt) const {
            MADNESS_EXCEPTION("gradient method is not yet implemented for this shape",0);
        }
    };

    class SDFSphere : public SignedDFInterface<3> {
    protected:
        const double radius; 
        const coord_3d center; 

    public:
        SDFSphere(const double radius, const coord_3d &center) 
            : radius(radius)
            , center(center) 
        {}

        double sdf(const coord_3d& pt) const {
            double temp, r;
            int i;
            
            r = 0.0;
            for(i = 0; i < 3; ++i) {
                temp = pt[i] - center[i];
                r += temp * temp;
            }
            
            return sqrt(r) - radius;
        }

        coord_3d grad_sdf(const coord_3d& pt) const {
            double x = pt[0] - center[0];
            double y = pt[1] - center[1];
            double z = pt[2] - center[2];
            double r = sqrt(x*x + y*y + z*z);
            coord_3d g;
            g[0] = x/r;
            g[1] = y/r;
            g[2] = z/r;
            return g;
        }
    };

    class SDFCone : public SignedDFInterface<3> {
    protected:
        const coord_3d apex; 
        const double c; 
        const coord_3d dir; 

    public:
        SDFCone(const double c, const coord_3d &apex, const coord_3d &direc) 
            : apex(apex)
            , c(c)
            , dir(direc*(1.0/sqrt(direc[0]*direc[0] + direc[1]*direc[1] + direc[2]*direc[2])))
        {}

        double sdf(const coord_3d& pt) const {
            coord_3d diff;
            unsigned int i;
            double dotp;
            
            for(i = 0; i < 3; ++i)
                diff[i] = pt[i] - apex[i];
            
            dotp = diff[0]*dir[0] + diff[1]*dir[1] + diff[2]*dir[2];
            
            for(i = 0; i < 3; ++i)
                diff[i] -= dotp * dir[i];
            
            return sqrt(diff[0]*diff[0] + diff[1]*diff[1] + diff[2]*diff[2])
                - c * dotp;
        }

        coord_3d grad_sdf(const coord_3d& pt) const {
            MADNESS_EXCEPTION("gradient method is not yet implemented for this shape",0);
        }
    };

    class SDFParaboloid : public SignedDFInterface<3> {
    protected:
        const coord_3d apex; 
        const double c; 
        const coord_3d dir;
        const long double zero; 
        const long double rootzero; 

    public:
        SDFParaboloid(const double c, const coord_3d &apex, const coord_3d &direc) 
            : apex(apex)
            , c(c)
            , dir(direc*(1.0/sqrt(direc[0]*direc[0] + direc[1]*direc[1] + direc[2]*direc[2])))
            , zero(1.0e-10L)
            , rootzero(1.0e-14L)
        {}

        double sdf(const coord_3d& pt) const {
            // work in extended precision since this calculation can have
            // large cancellations
            Vector<long double, 3> diff;
            unsigned int i, n;
            long double dotp, d;
            long double lambda, dist, temp[2], upper;
            std::vector<long double> roots;

            // silence warnings
            lambda = 0.0;
            dist = 0.0;

            // avoid lots of casting
            long double c = (long double)this->c;

            for(i = 0; i < 3; ++i)
                diff[i] = pt[i] - apex[i];

            dotp = diff[0]*dir[0] + diff[1]*dir[1] + diff[2]*dir[2];
            
            for(i = 0; i < 3; ++i)
                diff[i] -= dotp * dir[i];

            d = diff[0]*diff[0] + diff[1]*diff[1] + diff[2]*diff[2]
                - c * dotp;
            //printf("   %.16Le %.16Le %.16Le\n", c, d, dotp);

            // get the real roots (Lagrange multipliers) and find the one
            // that minimizes the distance.
            // note that real roots must be no bigger than (4dotp+c)/(2c)
            // from the analytical form
            roots = get_roots(c, d, dotp);
            upper = (4.0L*dotp + c) * 0.5L / c;
            n = 0;
            for(std::vector<long double>::iterator iter = roots.begin();
                iter != roots.end(); ++iter) {

                temp[0] = *iter;

                // are we in range?
                if(temp[0] <= upper) {

                    // calculate the distance
                    temp[1] = c*(dotp - temp[0]*c*0.5L + c*0.25L);
                    if(temp[1] < 0.0L) {
                        // temp[1] is analytically non-negative
                        // if we're here, it's noise
                        temp[1] = 0.0L;
                    }
                    else {
                        temp[1] = fabs(temp[0])*sqrt(temp[1]);
                    }

                    if(n == 0) {
                        lambda = temp[0];
                        dist = temp[1];
                    }
                    else {
                        if(temp[1] < dist) {
                            lambda = temp[0];
                            dist = temp[1];
                        }
                    }
                    ++n;
                }
            }
            //printf("   L = %.16Le D = %.16Le\n", lambda, dist);

            MADNESS_ASSERT(n > 0); // we should have at least one real root...

            if(d > 0.0L)
                return (double)dist;
            else
                return (double)(-dist);
        }

        coord_3d grad_sdf(const coord_3d& pt) const {
            MADNESS_EXCEPTION("gradient method is not yet implemented for this shape",0);
        }

        protected:
        std::vector<long double> get_roots(const long double c,
                             const long double d, const long double z) const {

            long double lower, upper, bound;
            std::vector<long double> roots(0);
            long double derivroot1, derivroot2, temp;
            bool hasregion2, keepgoing, foundroot;

            // calculate the region boundaries
            lower = 2.0L*(c+z);
            upper = fabs(c - 2.0L*z);
            hasregion2 = (upper >= 1.0e-10L);

            derivroot1 = (lower - upper) / (3.0L*c);
            derivroot2 = (lower + upper) / (3.0L*c);

            // first region -------------------------------------------------
            foundroot = false;
            upper = derivroot1;

            // if cubic(upper) < 0.0, there's no root in this region
            bound = eval_cubic(upper, c, d, z);
            if(fabs(bound) < zero) {
                // this is the root!
                roots.push_back(upper);
                hasregion2 = false;
            }
            if(bound > 0.0L) {
                // there's a root in region 1: find it!

                // need to find a lower bound
                temp = -1.0;
                do {
                    lower = upper + temp;
                    bound = eval_cubic(lower, c, d, z);
                    if(fabs(bound) < zero) {
                        // found the root!
                        roots.push_back(lower);
                        foundroot = true;
                        break;
                    }

                    keepgoing = (bound > 0.0L);
                    // if we didn't cross the root, we have a better upperbound
                    if(keepgoing)
                        upper = lower;
                } while(keepgoing);

                // with a lower and upper bound, find the root!
                if(!foundroot) {
                    temp = find_root(lower, upper, c, d, z, true);
                    roots.push_back(temp);
                }
            }
            else {
                // no root in region 1; hence no root in region 2
                hasregion2 = false;
            }

            // second region ------------------------------------------------
            if(hasregion2) {
                lower = derivroot1;
                upper = derivroot2;
                // the lower range has to be positive the upper range has to be
                // negative
                if(lower > 0.0L && upper < 0.0L) {
                    temp = find_root(lower, upper, c, d, z, false);
                    roots.push_back(temp);
                }
            }

            // third region -------------------------------------------------
            foundroot = false;
            lower = derivroot2;

            // if cubic(upper) > 0.0, there's no root in this region
            bound = eval_cubic(lower, c, d, z);
            if(fabs(bound) < zero) {
                // this is the root!
                roots.push_back(lower);
            }
            if(bound < 0.0L) {
                // there's a root in region 3: find it!

                // need to find an upper bound
                temp = 1.0;
                do {
                    upper = lower + temp;
                    bound = eval_cubic(upper, c, d, z);
                    if(fabs(bound) < zero) {
                        // found the root!
                        roots.push_back(upper);
                        foundroot = true;
                        break;
                    }

                    keepgoing = (bound < 0.0L);
                    // if we didn't cross the root, we have a better lowerbound
                    if(keepgoing)
                        lower = upper;
                } while(keepgoing);

                // with a lower and upper bound, find the root!
                if(!foundroot) {
                    temp = find_root(lower, upper, c, d, z, true);
                    roots.push_back(temp);
                }
            }

            return roots;
        }

        long double find_root(long double lower, long double upper,
                              const long double c, const long double d,
                              const long double z, bool dir) const {

            long double value, middle;
            do {
                middle = 0.5L*(upper + lower);
                value = eval_cubic(middle, c, d, z);

                if(fabs(value) < zero)
                    return middle;

                if(value < 0.0L) {
                    if(dir)
                        lower = middle;
                    else
                        upper = middle;
                }
                else {
                    if(dir)
                        upper = middle; 
                    else
                        lower = middle;
                }
            } while(fabs(upper-lower) > rootzero);

            return middle;
        }

        long double eval_cubic(const long double x, const long double c,
                               const long double d, const long double z)
                    const {
            return ((c*0.5L*x - (c + z))*x + (2.0L*z + c*0.5L))*c*x + d;
        }

        long double eval_cubic_deriv(const long double x, const long double c,
                               const long double d, const long double z)
                    const {
            return c*((1.5L*c*x - 2.0L*(c+z))*x + 2.0L*z + 0.5L*c);
        }
    };

    class SDFBox : public SignedDFInterface<3> {
    protected:
        const coord_3d lengths;  
        const coord_3d center; 

    public:
        SDFBox(const coord_3d& length, const coord_3d& center) 
            : lengths(length*0.5), center(center) 
        {}

        double sdf(const coord_3d& pt) const {
            double diff, max;
            int i;
            
            max = fabs(pt[0] - center[0]) - lengths[0];
            for(i = 1; i < 3; ++i) {
                diff = fabs(pt[i] - center[i]) - lengths[i];
                if(diff > max)
                    max = diff;
            }
            
            return max;
        }

        coord_3d grad_sdf(const coord_3d& pt) const {
            MADNESS_EXCEPTION("gradient method is not yet implemented for this shape",0);
        }
    };

    class SDFCube : public SDFBox {
    public:
        SDFCube(const double length, const coord_3d& center) 
            : SDFBox(coord_3d(length), center)
        {}
    };

    class SDFEllipsoid : public SignedDFInterface<3> {
    protected:
        coord_3d radii; 
        coord_3d center; 
        
    public:
        SDFEllipsoid(const coord_3d& radii, const coord_3d& center) 
            : radii(radii)
            , center(center) 
        {}
        
        double sdf(const coord_3d& pt) const {
            double quot, sum;
            int i;
            
            sum = 0.0;
            for(i = 0; i < 3; ++i) {
                quot = (pt[i] - center[i]) / radii[i];
                sum += quot * quot;
            }
            
            return sum - 1.0;
        }

        coord_3d grad_sdf(const coord_3d& pt) const {
            MADNESS_EXCEPTION("gradient method is not yet implemented for this shape",0);
        }
    };
    
    class SDFCylinder : public SignedDFInterface<3> {
    protected:
        double radius; 
        double a; 
        coord_3d center; 
        coord_3d axis; 

    public:
        SDFCylinder(const double radius, const double length, const coord_3d& axpt, const coord_3d& axis) 
            : radius(radius)
            , a(length / 2.0)
            , center(axpt)
            , axis(axis*(1.0/sqrt(axis[0]*axis[0] + axis[1]*axis[1] + axis[2]*axis[2]))) 
        {}
        
        double sdf(const coord_3d& pt) const {
            double dist;
            coord_3d rel, radial;
            int i;
            
            // axial distance
            dist = 0.0;
            for(i = 0; i < 3; ++i) {
                rel[i] = pt[i] - center[i];
                dist += rel[i] * axis[i];
            }
            
            // get the radial component
            for(i = 0; i < 3; ++i)
                radial[i] = rel[i] - dist * axis[i];
            
            return std::max(fabs(dist) - a, sqrt(radial[0]*radial[0] + radial[1]*radial[1]
                                                 + radial[2]*radial[2]) - radius);
        }

        coord_3d grad_sdf(const coord_3d& pt) const {
            MADNESS_EXCEPTION("gradient method is not yet implemented for this shape",0);
        }
    };

} // end of madness namespace

#endif // MADNESS_MRA_SDF_SHAPE_3D_H__INCLUDED