Commit 46d9ec36 authored by COULAUD Olivier's avatar COULAUD Olivier

Add data/methods to use non cuboid geometry

parent a036c5dd
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
// Copyright ScalFmm 2011 INRIA
// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
// This software is a computer program whose purpose is to compute the FMM.
//
......
......@@ -28,8 +28,9 @@ enum KERNEL_FUNCTION_IDENTIFIER {ONE_OVER_R,
ONE_OVER_R_SQUARED,
LENNARD_JONES_POTENTIAL,
R_IJ,
R_IJK};
R_IJK,
ONE_OVER_RH
};
// probably not extendable :)
enum KERNEL_FUNCTION_TYPE {HOMOGENEOUS, NON_HOMOGENEOUS};
......@@ -50,7 +51,7 @@ enum KERNEL_FUNCTION_TYPE {HOMOGENEOUS, NON_HOMOGENEOUS};
* Let there be a kernel \f$K\f$ such that \f$X_i=K_{ij}Y_j\f$
* with \f$X\f$ the lhs of size NLHS and \f$Y\f$ the rhs of size NRHS.
* The table applyTab provides the indices in the reduced storage table
* corresponding to the application scheme depicted ealier.
* corresponding to the application scheme depicted earlier.
*
* PB: BEWARE! Homogeneous matrix kernels do not support cell width extension
* yet. Is it possible to find a reference width and a scale factor such that
......@@ -113,10 +114,18 @@ struct FInterpMatrixKernelR : FInterpAbstractMatrixKernel
/// One over r when the box size is rescaled to 1
struct FInterpMatrixKernelRH :FInterpMatrixKernelR{
static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
static const KERNEL_FUNCTION_IDENTIFIER Identifier = ONE_OVER_RH;
static const unsigned int NCMP = 1; //< number of components
static const unsigned int NPV = 1; //< dim of physical values
static const unsigned int NPOT = 1; //< dim of potentials
static const unsigned int NRHS = 1; //< dim of mult exp
static const unsigned int NLHS = 1; //< dim of loc exp
FReal LX,LY,LZ ;
FInterpMatrixKernelRH(const FReal = 0.0, const unsigned int = 0) : FInterpMatrixKernelR(),
LX(1.0),LY(1.0),LZ(1.0)
{}
FInterpMatrixKernelRH(const FReal = 0.0, const unsigned int = 0) : LX(1.0),LY(1.0),LZ(1.0)
{ }
FReal evaluate(const FPoint& x, const FPoint& y) const
{
const FPoint xy(x-y);
......@@ -126,6 +135,26 @@ struct FInterpMatrixKernelRH :FInterpMatrixKernelR{
}
void setCoeff(const FReal& a, const FReal& b, const FReal& c)
{LX= a ; LY = b ; LZ = c ;}
// returns position in reduced storage
int getPosition(const unsigned int) const
{return 0;}
void evaluateBlock(const FPoint& x, const FPoint& y, FReal* block) const
{
block[0]=this->evaluate(x,y);
}
FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
{
const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
return getScaleFactor(CellWidth);
}
FReal getScaleFactor(const FReal CellWidth) const
{
return FReal(2.) / CellWidth;
}
};
......
......@@ -29,6 +29,25 @@ struct DirectInteractionComputer<ONE_OVER_R, 1>
}
};
//! Specialization for Laplace potential on Non uniform domain
template <>
struct DirectInteractionComputer<ONE_OVER_RH, 1>
{
template <typename ContainerClass>
static void P2P( ContainerClass* const FRestrict TargetParticles,
ContainerClass* const NeighborSourceParticles[27]){
FP2P::FullMutual(TargetParticles,NeighborSourceParticles,14);
}
//
template <typename ContainerClass>
static void P2PRemote( ContainerClass* const FRestrict inTargets,
ContainerClass* const inNeighbors[27],
const int inSize){
FP2P::FullRemote(inTargets,inNeighbors,inSize);
}
};
/*! Specialization for Lennard-Jones potential */
template <>
......
......@@ -41,7 +41,7 @@ private:
public:
/** Default constructor (sets position to 0/0/0) */
FPoint(){
data[0] = data[1] = data[2] = FReal(0);
data[0] = data[1] = data[2] = FReal(0.0);
}
/** Constructor from an array */
......@@ -243,7 +243,7 @@ public:
/**
* Affect to all dim the other position
* @param other the value to afect
* @param other the value to affect
* @return the current object after being affected
*/
FPoint& operator+=(const FPoint& other){
......@@ -252,10 +252,20 @@ public:
this->data[2] += other.data[2];
return *this;
}
/**
* Divide to all dim the other position
* @param other the value to affect
* @return the current object after being affected
*/
FPoint& operator/=(const FPoint& other){
this->data[0] /= other.data[0];
this->data[1] /= other.data[1];
this->data[2] /= other.data[2];
return *this;
}
/**
* Affect to all dim the other position
* @param other the value to afect
* @param other the value to affect
* @return the current object after being affected
*/
FPoint& operator*=(const FReal value){
......@@ -289,7 +299,7 @@ public:
/**
* Operator F3Position minus F3Position
* This substract one from anther
* This subtract one from anther
* @param inPosition the position to reduce
* @param inOther the position to decrease/substract inPosition
* @return the resulting position
......@@ -300,7 +310,7 @@ public:
/**
* Operator F3Position plus F3Position
* This substract one from anther
* This subtract one from anther
* @param inPosition the position to reduce
* @param inOther the position to increase inPosition
* @return the resulting position
......@@ -309,6 +319,7 @@ public:
return FPoint(inPosition.data[0] + inOther.data[0], inPosition.data[1] + inOther.data[1], inPosition.data[2] + inOther.data[2]);
}
/**
* Operator stream FPoint to std::ostream
* This can be used to simpldata[1] write out a position
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment