FChebInterpolator.hpp 14 KB
Newer Older
1 2 3 4 5 6 7 8
#ifndef FCHEBINTERPOLATOR_HPP
#define FCHEBINTERPOLATOR_HPP


#include "./FChebMapping.hpp"
#include "./FChebTensor.hpp"
#include "./FChebRoots.hpp"

9
#include "../../Utils/FBlas.hpp"
10

11 12 13 14 15 16 17 18 19


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

/**
 * @class FChebInterpolator
20
 *
21
 * The class @p FChebInterpolator defines the anterpolation (M2M) and
22
 * interpolation (L2L) concerning operations.
23 24 25 26 27 28 29 30 31 32 33
 */
template <int ORDER>
class FChebInterpolator : FNoCopyable
{
  // compile time constants and types
  enum {nnodes = TensorTraits<ORDER>::nnodes};
  typedef FChebRoots< ORDER>  BasisType;
  typedef FChebTensor<ORDER> TensorType;

  FReal T_of_roots[ORDER][ORDER];
	unsigned int node_ids[nnodes][3];
34 35 36 37 38 39 40
	FReal* ChildParentInterpolator[8];

	/**
	 * Initialize the child - parent - interpolator, it is basically the matrix
	 * S which is precomputed and reused for all M2M and L2L operations, ie for
	 * all non leaf inter/anterpolations.
	 */
41
	void initM2MandL2L()
42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
	{
		F3DPosition ParentRoots[nnodes], ChildRoots[nnodes];
		const FReal ParentWidth(2.);
		const F3DPosition ParentCenter(0., 0., 0.);
		FChebTensor<ORDER>::setRoots(ParentCenter, ParentWidth, ParentRoots);

		F3DPosition ChildCenter;
		const FReal ChildWidth(1.);
		
		// loop: child cells
		for (unsigned int child=0; child<8; ++child) {

			// allocate memory
			ChildParentInterpolator[child] = new FReal [nnodes * nnodes];

			// set child info
			FChebTensor<ORDER>::setRelativeChildCenter(child, ChildCenter);
			FChebTensor<ORDER>::setRoots(ChildCenter, ChildWidth, ChildRoots);

			// assemble child - parent - interpolator
			assembleInterpolator(nnodes, ChildRoots, ChildParentInterpolator[child]);
		}
	}

66 67 68 69


public:
	/**
70
	 * Constructor: Initialize the Chebyshev polynomials at the Chebyshev
71 72 73 74 75 76 77
	 * roots/interpolation point
	 */
	explicit FChebInterpolator()
	{
		// initialize chebyshev polynomials of root nodes: T_o(x_j)
    for (unsigned int o=1; o<ORDER; ++o)
      for (unsigned int j=0; j<ORDER; ++j)
messner's avatar
messner committed
78
        T_of_roots[o][j] = FReal(BasisType::T(o, FReal(BasisType::roots[j])));
79 80 81

		// initialize root node ids
		TensorType::setNodeIds(node_ids);
82 83 84

		// initialize interpolation operator for non M2M and L2L (non leaf
		// operations)
85
		this -> initM2MandL2L();
86 87 88 89 90 91 92 93 94 95
	}

	
	/**
	 * Destructor: Delete dynamically allocated memory for M2M and L2L operator
	 */
	~FChebInterpolator()
	{
		for (unsigned int child=0; child<8; ++child)
			delete [] ChildParentInterpolator[child];
96 97 98
	}


99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139
	/**
	 * Assembles the interpolator \f$S_\ell\f$ of size \f$N\times
	 * \ell^3\f$. Here local points is meant as points whose global coordinates
	 * have already been mapped to the reference interval [-1,1].
	 *
	 * @param[in] NumberOfLocalPoints
	 * @param[in] LocalPoints
	 * @param[out] Interpolator
	 */
	void assembleInterpolator(const unsigned int NumberOfLocalPoints,
														const F3DPosition *const LocalPoints,
														FReal *const Interpolator) const
	{
		// values of chebyshev polynomials of source particle: T_o(x_i)
		FReal T_of_x[ORDER][3];

		// loop: local points (mapped in [-1,1])
		for (unsigned int m=0; m<NumberOfLocalPoints; ++m) {
			// evaluate chebyshev polynomials at local points
			for (unsigned int o=1; o<ORDER; ++o) {
				T_of_x[o][0] = BasisType::T(o, LocalPoints[m].getX());
				T_of_x[o][1] = BasisType::T(o, LocalPoints[m].getY());
				T_of_x[o][2] = BasisType::T(o, LocalPoints[m].getZ());
			}
			
			// assemble interpolator
			for (unsigned int n=0; n<nnodes; ++n) {
				Interpolator[n*nnodes + m] = FReal(1.);
				for (unsigned int d=0; d<3; ++d) {
					const unsigned int j = node_ids[n][d];
					FReal S_d = FReal(1.) / ORDER;
					for (unsigned int o=1; o<ORDER; ++o)
						S_d += FReal(2.) / ORDER * T_of_x[o][d] * T_of_roots[o][j];
					Interpolator[n*nnodes + m] *= S_d;
				}
			}
			
		}
		
	}

140 141
	
	/**
142 143
	 * Particle to moment: application of \f$S_\ell(y,\bar y_n)\f$
	 * (anterpolation, it is the transposed interpolation)
144 145
	 */
	template <class ContainerClass>
146 147 148
	void applyP2M(const F3DPosition& center,
								const FReal width,
								FReal *const multipoleExpansion,
149
								const ContainerClass *const sourceParticles) const;
150 151 152 153


	
	/**
154
	 * Local to particle operation: application of \f$S_\ell(x,\bar x_m)\f$ (interpolation)
155 156
	 */
	template <class ContainerClass>
157 158 159
	void applyL2P(const F3DPosition& center,
								const FReal width,
								const FReal *const localExpansion,
160
								ContainerClass *const localParticles) const;
161

162 163 164 165 166 167 168 169

	/**
	 * Local to particle operation: application of \f$\nabla_x S_\ell(x,\bar x_m)\f$ (interpolation)
	 */
	template <class ContainerClass>
	void applyL2PGradient(const F3DPosition& center,
												const FReal width,
												const FReal *const localExpansion,
170
												ContainerClass *const localParticles) const;
171

172 173 174 175 176 177 178
	
	void applyM2M(const unsigned int ChildIndex,
								const FReal *const ChildExpansion,
								FReal *const ParentExpansion) const
	{
		FBlas::gemtva(nnodes, nnodes, FReal(1.),
									ChildParentInterpolator[ChildIndex],
messner's avatar
messner committed
179
									const_cast<FReal*>(ChildExpansion), ParentExpansion);
180
	}
181

182 183 184 185 186 187
	void applyL2L(const unsigned int ChildIndex,
								const FReal *const ParentExpansion,
								FReal *const ChildExpansion) const
	{
		FBlas::gemva(nnodes, nnodes, FReal(1.),
								 ChildParentInterpolator[ChildIndex],
messner's avatar
messner committed
188
								 const_cast<FReal*>(ParentExpansion), ChildExpansion);
189 190 191
	}
	
};
192 193 194



195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393




/**
 * Particle to moment: application of \f$S_\ell(y,\bar y_n)\f$
 * (anterpolation, it is the transposed interpolation)
 */
template <int ORDER>
template <class ContainerClass>
inline void FChebInterpolator<ORDER>::applyP2M(const F3DPosition& center,
																							 const FReal width,
																							 FReal *const multipoleExpansion,
																							 const ContainerClass *const sourceParticles) const
{
	// set all multipole expansions to zero
	FBlas::setzero(nnodes, multipoleExpansion);

	// allocate stuff
	const map_glob_loc map(center, width);
	F3DPosition localPosition;
	FReal T_of_x[ORDER][3];
	FReal S[3];

	// loop over source particles
	typename ContainerClass::ConstBasicIterator iter(*sourceParticles);
	while(iter.hasNotFinished()){

		// map global position to [-1,1]
		map(iter.data().getPosition(), localPosition);

		// evaluate chebyshev polynomials of source particle: T_o(x_i)
		T_of_x[0][0] = FReal(1.);	T_of_x[1][0] = localPosition.getX();
		T_of_x[0][1] = FReal(1.);	T_of_x[1][1] = localPosition.getY();
		T_of_x[0][2] = FReal(1.);	T_of_x[1][2] = localPosition.getZ();
		for (unsigned int o=2; o<ORDER; ++o) {
			T_of_x[o][0] = FReal(2.) * localPosition.getX() * T_of_x[o-1][0] - T_of_x[o-2][0];
			T_of_x[o][1] = FReal(2.) * localPosition.getY() * T_of_x[o-1][1] - T_of_x[o-2][1];
			T_of_x[o][2] = FReal(2.) * localPosition.getZ() * T_of_x[o-1][2] - T_of_x[o-2][2];
		}
		
		// anterpolate
		const FReal sourceValue = iter.data().getPhysicalValue();
		for (unsigned int n=0; n<nnodes; ++n) {
			const unsigned int j[3] = {node_ids[n][0], node_ids[n][1], node_ids[n][2]};
			S[0] = T_of_x[1][0] * T_of_roots[1][j[0]];
			S[1] = T_of_x[1][1] * T_of_roots[1][j[1]];
			S[2] = T_of_x[1][2] * T_of_roots[1][j[2]];
			for (unsigned int o=2; o<ORDER; ++o) {
				S[0] += T_of_x[o][0] * T_of_roots[o][j[0]];
				S[1] += T_of_x[o][1] * T_of_roots[o][j[1]];
				S[2] += T_of_x[o][2] * T_of_roots[o][j[2]];
			}
			// gather contributions
			S[0] *= FReal(2.); S[0] += FReal(1.);
			S[1] *= FReal(2.); S[1] += FReal(1.);
			S[2] *= FReal(2.); S[2] += FReal(1.);
			multipoleExpansion[n]	+= FReal(1.) / nnodes *	S[0] * S[1] * S[2] *	sourceValue;
		}
		// increment source iterator
		iter.gotoNext();
	}
}


/**
 * Local to particle operation: application of \f$S_\ell(x,\bar x_m)\f$ (interpolation)
 */
template <int ORDER>
template <class ContainerClass>
inline void FChebInterpolator<ORDER>::applyL2P(const F3DPosition& center,
																							 const FReal width,
																							 const FReal *const localExpansion,
																							 ContainerClass *const localParticles) const
{
	// allocate stuff
	const map_glob_loc map(center, width);
	F3DPosition localPosition;
	FReal T_of_x[ORDER][3];
	FReal S[3];
		
	typename ContainerClass::BasicIterator iter(*localParticles);
	while(iter.hasNotFinished()){
			
		// map global position to [-1,1]
		map(iter.data().getPosition(), localPosition);

		// evaluate chebyshev polynomials of source particle: T_o(x_i)
		T_of_x[0][0] = FReal(1.);	T_of_x[1][0] = localPosition.getX();
		T_of_x[0][1] = FReal(1.);	T_of_x[1][1] = localPosition.getY();
		T_of_x[0][2] = FReal(1.);	T_of_x[1][2] = localPosition.getZ();
		for (unsigned int o=2; o<ORDER; ++o) {
			T_of_x[o][0] = FReal(2.) * localPosition.getX() * T_of_x[o-1][0] - T_of_x[o-2][0];
			T_of_x[o][1] = FReal(2.) * localPosition.getY() * T_of_x[o-1][1] - T_of_x[o-2][1];
			T_of_x[o][2] = FReal(2.) * localPosition.getZ() * T_of_x[o-1][2] - T_of_x[o-2][2];
		}

		// interpolate and increment target value
		FReal targetValue = iter.data().getPotential();
		for (unsigned int n=0; n<nnodes; ++n) {
			const unsigned int j[3] = {node_ids[n][0], node_ids[n][1], node_ids[n][2]};
			S[0] = T_of_x[1][0] * T_of_roots[1][j[0]];
			S[1] = T_of_x[1][1] * T_of_roots[1][j[1]];
			S[2] = T_of_x[1][2] * T_of_roots[1][j[2]];
			for (unsigned int o=2; o<ORDER; ++o) {
				S[0] += T_of_x[o][0] * T_of_roots[o][j[0]];
				S[1] += T_of_x[o][1] * T_of_roots[o][j[1]];
				S[2] += T_of_x[o][2] * T_of_roots[o][j[2]];
			}
			// gather contributions
			S[0] *= FReal(2.); S[0] += FReal(1.);
			S[1] *= FReal(2.); S[1] += FReal(1.);
			S[2] *= FReal(2.); S[2] += FReal(1.);
			targetValue	+= FReal(1.) / nnodes *	S[0] * S[1] * S[2] * localExpansion[n];
		}
		// set potential
		iter.data().setPotential(targetValue);
		// increment target iterator
		iter.gotoNext();
	}
}






/**
 * Local to particle operation: application of \f$\nabla_x S_\ell(x,\bar x_m)\f$ (interpolation)
 */
template <int ORDER>
template <class ContainerClass>
inline void FChebInterpolator<ORDER>::applyL2PGradient(const F3DPosition& center,
																											 const FReal width,
																											 const FReal *const localExpansion,
																											 ContainerClass *const localParticles) const
{
	// setup local to global mapping
	const map_glob_loc map(center, width);
	F3DPosition Jacobian;
	map.computeJacobian(Jacobian);
	const FReal jacobian[3] = {Jacobian.getX(), Jacobian.getY(), Jacobian.getZ()}; 
	F3DPosition localPosition;
	FReal T_of_x[ORDER][3];
	FReal U_of_x[ORDER][3];
	FReal P[3];

	typename ContainerClass::BasicIterator iter(*localParticles);
	while(iter.hasNotFinished()){
			
		// map global position to [-1,1]
		map(iter.data().getPosition(), localPosition);
			
		// evaluate chebyshev polynomials of source particle
		// T_0(x_i) and T_1(x_i)
		T_of_x[0][0] = FReal(1.);	T_of_x[1][0] = localPosition.getX();
		T_of_x[0][1] = FReal(1.);	T_of_x[1][1] = localPosition.getY();
		T_of_x[0][2] = FReal(1.);	T_of_x[1][2] = localPosition.getZ();
		// U_0(x_i) and U_1(x_i)
		U_of_x[0][0] = FReal(1.);	U_of_x[1][0] = localPosition.getX() * FReal(2.);
		U_of_x[0][1] = FReal(1.);	U_of_x[1][1] = localPosition.getY() * FReal(2.);
		U_of_x[0][2] = FReal(1.);	U_of_x[1][2] = localPosition.getZ() * FReal(2.);
		for (unsigned int o=2; o<ORDER; ++o) {
			// T_o(x_i)
			T_of_x[o][0] = FReal(2.)*localPosition.getX()*T_of_x[o-1][0] - T_of_x[o-2][0];
			T_of_x[o][1] = FReal(2.)*localPosition.getY()*T_of_x[o-1][1] - T_of_x[o-2][1];
			T_of_x[o][2] = FReal(2.)*localPosition.getZ()*T_of_x[o-1][2] - T_of_x[o-2][2];
			// U_o(x_i)
			U_of_x[o][0] = FReal(2.)*localPosition.getX()*U_of_x[o-1][0] - U_of_x[o-2][0];
			U_of_x[o][1] = FReal(2.)*localPosition.getY()*U_of_x[o-1][1] - U_of_x[o-2][1];
			U_of_x[o][2] = FReal(2.)*localPosition.getZ()*U_of_x[o-1][2] - U_of_x[o-2][2];
		}

		// scale, because dT_o/dx = oU_{o-1}
		for (unsigned int o=2; o<ORDER; ++o) {
			U_of_x[o-1][0] *= FReal(o);
			U_of_x[o-1][1] *= FReal(o);
			U_of_x[o-1][2] *= FReal(o);
		}

		// apply P and increment forces
		FReal forces[3] = {FReal(0.), FReal(0.), FReal(0.)};
		for (unsigned int n=0; n<nnodes; ++n) {
			
			// tensor indices of chebyshev nodes
			const unsigned int j[3] = {node_ids[n][0], node_ids[n][1], node_ids[n][2]};

			// f0 component //////////////////////////////////////
			P[0] = U_of_x[0][0] * T_of_roots[1][j[0]];
			P[1] = T_of_x[1][1] * T_of_roots[1][j[1]];
			P[2] = T_of_x[1][2] * T_of_roots[1][j[2]];
			for (unsigned int o=2; o<ORDER; ++o) {
				P[0] += U_of_x[o-1][0] * T_of_roots[o][j[0]];
				P[1] += T_of_x[o  ][1] * T_of_roots[o][j[1]];
				P[2] += T_of_x[o  ][2] * T_of_roots[o][j[2]];
			}
			P[0] *= FReal(2.);
			P[1] *= FReal(2.); P[1] += FReal(1.);
			P[2] *= FReal(2.); P[2] += FReal(1.);
394
			forces[0]	+= P[0] * P[1] * P[2] * localExpansion[n];
395 396 397 398 399 400 401 402 403 404 405 406 407

			// f1 component //////////////////////////////////////
			P[0] = T_of_x[1][0] * T_of_roots[1][j[0]];
			P[1] = U_of_x[0][1] * T_of_roots[1][j[1]];
			P[2] = T_of_x[1][2] * T_of_roots[1][j[2]];
			for (unsigned int o=2; o<ORDER; ++o) {
				P[0] += T_of_x[o  ][0] * T_of_roots[o][j[0]];
				P[1] += U_of_x[o-1][1] * T_of_roots[o][j[1]];
				P[2] += T_of_x[o  ][2] * T_of_roots[o][j[2]];
			}
			P[0] *= FReal(2.); P[0] += FReal(1.);
			P[1] *= FReal(2.); 
			P[2] *= FReal(2.); P[2] += FReal(1.);
408
			forces[1]	+= P[0] * P[1] * P[2] * localExpansion[n];
409 410 411 412 413 414 415 416 417 418 419 420 421

			// f2 component //////////////////////////////////////
			P[0] = T_of_x[1][0] * T_of_roots[1][j[0]];
			P[1] = T_of_x[1][1] * T_of_roots[1][j[1]];
			P[2] = U_of_x[0][2] * T_of_roots[1][j[2]];
			for (unsigned int o=2; o<ORDER; ++o) {
				P[0] += T_of_x[o  ][0] * T_of_roots[o][j[0]];
				P[1] += T_of_x[o  ][1] * T_of_roots[o][j[1]];
				P[2] += U_of_x[o-1][2] * T_of_roots[o][j[2]];
			}
			P[0] *= FReal(2.); P[0] += FReal(1.);
			P[1] *= FReal(2.); P[1] += FReal(1.);
			P[2] *= FReal(2.);
422
			forces[2]	+= P[0] * P[1] * P[2] * localExpansion[n];
423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440
		}

		// scale forces
		forces[0] *= jacobian[0] / nnodes;
		forces[1] *= jacobian[1] / nnodes;
		forces[2] *= jacobian[2] / nnodes;

		// set computed forces
		iter.data().incForces(forces[0] * iter.data().getPhysicalValue(),
													forces[1] * iter.data().getPhysicalValue(),
													forces[2] * iter.data().getPhysicalValue());

		// increment iterator
		iter.gotoNext();
	}
}


441
#endif