FChebMapping.hpp 5 KB
Newer Older
BRAMAS Berenger's avatar
BRAMAS Berenger committed
1
// ===================================================================================
2 3 4 5 6 7 8 9 10 11 12 13 14
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.  
// 
// 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 and CeCILL-C Licenses for more details.
// "http://www.cecill.info". 
// "http://www.gnu.org/licenses".
BRAMAS Berenger's avatar
BRAMAS Berenger committed
15
// ===================================================================================
16 17 18 19 20
#ifndef FCHEBMAPPING_HPP
#define FCHEBMAPPING_HPP

#include <limits>

21
#include "../../Utils/FNoCopyable.hpp"
22
#include "../../Utils/FPoint.hpp"
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

/**
 * @author Matthias Messner (matthias.matthias@inria.fr)
 * Please read the license
 */

/**
 * @class FChebMapping
 *
 * The class @p FChebMapping is the base class for the affine mapping
 * \f$\Phi:[-1,1]\rightarrow[a,b] \f$ and the inverse affine mapping
 * \f$\Phi^{-1}:[a,b]\rightarrow[-1,1]\f$.
 */
class FChebMapping : FNoCopyable
{
protected:
39 40
  FPoint a;
  FPoint b;
41
 
42
  explicit FChebMapping(const FPoint& center,
43 44 45 46 47 48 49 50
											 const FReal width)
		: a(center.getX() - width / FReal(2.),
				center.getY() - width / FReal(2.),
				center.getZ() - width / FReal(2.)),
			b(center.getX() + width / FReal(2.),
				center.getY() + width / FReal(2.),
				center.getZ() + width / FReal(2.)) {}

51
  virtual void operator()(const FPoint&, FPoint&) const = 0;
52 53 54 55 56 57 58 59 60

public:
	/**
	 * Checks wheter @p position is within cluster, ie within @p a and @p b, or
	 * not.
	 *
	 * @param[in] position position (eg of a particle)
	 * @return @p true if position is in cluster else @p false
	 */
61
	bool is_in_cluster(const FPoint& position) const
62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99
	{
		// Set numerical limit
		const FReal epsilon = FReal(10.) * std::numeric_limits<FReal>::epsilon();
		// Return false if x-coordinate is not in cluster
		if (a.getX()-position.getX()>epsilon ||	position.getX()-b.getX()>epsilon) {
			std::cout << a.getX()-position.getX() << "\t"
								<< position.getX()-b.getX() << "\t"	<< epsilon << std::endl;
			return false;
		}
		// Return false if x-coordinate is not in cluster
		if (a.getY()-position.getY()>epsilon ||	position.getY()-b.getY()>epsilon) {
			std::cout << a.getY()-position.getY() << "\t"
								<< position.getY()-b.getY() << "\t"	<< epsilon << std::endl;
			return false;
		}
		// Return false if x-coordinate is not in cluster
		if (a.getZ()-position.getZ()>epsilon ||	position.getZ()-b.getZ()>epsilon) {
			std::cout << a.getZ()-position.getZ() << "\t"
								<< position.getZ()-b.getZ() << "\t"	<< epsilon << std::endl;
			return false;
		}
		// Position is in cluster, return true
		return true;
	}

};


/**
 * @class map_glob_loc
 *
 * This class defines the inverse affine mapping
 * \f$\Phi^{-1}:[a,b]\rightarrow[-1,1]\f$. It maps from global coordinates to
 * local ones.
 */
class map_glob_loc : public FChebMapping
{
public:
100
  explicit map_glob_loc(const FPoint& center, const FReal width)
101 102 103 104 105 106 107 108 109
		: FChebMapping(center, width) {}
  
	/**
	 * Maps from a global position to its local position: \f$\Phi^{-1}(x) =
	 * \frac{2x-b-a}{b-a}\f$.
	 *
	 * @param[in] globPos global position
	 * @param[out] loclPos local position
	 */
110
  void operator()(const FPoint& globPos, FPoint& loclPos) const
111
  {
112 113 114
		loclPos.setX((FReal(2.)*globPos.getX()-b.getX()-a.getX()) / (b.getX()-a.getX())); // 5 flops
		loclPos.setY((FReal(2.)*globPos.getY()-b.getY()-a.getY()) / (b.getY()-a.getY())); // 5 flops
		loclPos.setZ((FReal(2.)*globPos.getZ()-b.getZ()-a.getZ()) / (b.getZ()-a.getZ())); // 5 flops
115 116 117
	}

	// jacobian = 2 / (b - a);
118
	void computeJacobian(FPoint& jacobian) const
119
  {
120 121 122
		jacobian.setX(FReal(2.) / (b.getX() - a.getX())); // 2 flops
		jacobian.setY(FReal(2.) / (b.getY() - a.getY()));	// 2 flops
		jacobian.setZ(FReal(2.) / (b.getZ() - a.getZ())); // 2 flops
123 124 125 126 127 128 129 130 131 132 133 134 135
  }
};


/**
 * @class map_loc_glob
 *
 * This class defines the affine mapping \f$\Phi:[-1,1]\rightarrow[a,b]\f$. It
 * maps from local coordinates to global ones.
 */
class map_loc_glob : public FChebMapping
{
public:
136
  explicit map_loc_glob(const FPoint& center, const FReal width)
137 138 139 140 141 142 143 144 145 146
		: FChebMapping(center, width) {}
  
	// globPos = (a + b) / 2 + (b - a) * loclPos / 2;
	/**
	 * Maps from a local position to its global position: \f$\Phi(\xi) =
	 * \frac{1}{2}(a+b)+ \frac{1}{2}(b-a)\xi\f$.
	 *
	 * @param[in] loclPos local position
	 * @param[out] globPos global position
	 */
147
  void operator()(const FPoint& loclPos, FPoint& globPos) const
148 149 150 151 152 153 154 155 156 157 158 159 160
  {
		globPos.setX((a.getX()+b.getX())/FReal(2.)+
								 (b.getX()-a.getX())*loclPos.getX()/FReal(2.));
		globPos.setY((a.getY()+b.getY())/FReal(2.)+
								 (b.getY()-a.getY())*loclPos.getY()/FReal(2.));
		globPos.setZ((a.getZ()+b.getZ())/FReal(2.)+
								 (b.getZ()-a.getZ())*loclPos.getZ()/FReal(2.));
	}
};



#endif