Homec4science
Diffusion LibMultiScale (features/coupler_component) Browse

libmultiscale/src/md/spatial_filter_atomic.hh
1e97b3b33259features/coupler_compon…

Last Change
Raw File

spatial_filter_atomic.hh
View Options

/**
* @file spatial_filter_atomic.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch>
*
* @date Fri Oct 11 13:00:16 2013
*
* @brief This class allows to perform spatial filtering of atomic fields
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale is free software: you can redistribute it and/or modify it
* under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* LibMultiScale 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 Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
#ifndef __LIBMULTISCALE_SPATIAL_FILTER_ATOMIC_HH__
#define __LIBMULTISCALE_SPATIAL_FILTER_ATOMIC_HH__
/* -------------------------------------------------------------------------- */
#include "container_neighbor_atoms.hh"
#include "container_neighbor_atoms_interface.hh"
#include "cubic_spline.hh"
#include "hamming_window.hh"
#include "hanning_window.hh"
#include "math.h"
#include "ref_atom.hh"
#include "spatial_grid_libmultiscale.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
class SpatialFilterAtomic {
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
public:
enum WindowFunction { Hanning, Hamming, Rectangular };
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
SpatialFilterAtomic() {
cSpline = new CubicSpline;
kernel_length = 10;
kernel_distance = 10;
kernel_distance_sq = 10;
neighbors = NULL;
density = 1.0;
nReplica[0] = 0;
nReplica[1] = 0;
nReplica[2] = 0;
};
virtual ~SpatialFilterAtomic(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
//! set the density per unit length/area/volume
void setDensity(Real dens) { density = dens; };
//! set the number of replicas in each direction
void setNumberofReplica(UInt *replica) {
for (UInt i = 0; i < 3; ++i)
nReplica[i] = replica[i];
};
template <UInt Dim, typename Container>
void buildFilteredField(std::vector<Real> &filtered, FieldType ftype,
bool *pbc_flag, Container &container);
//! read memory kernel from a file
void readKernelFromFile(const std::string &filename);
//! read memory kernel from a file
template <UInt Dim, typename F>
void buildKernelFromFunction(F &foo, Real xmin, Real xmax, UInt n_points,
WindowFunction window = Rectangular);
//! build the Cubic Spline from a set of points
void buildCSpline(std::vector<double> &x, std::vector<double> &memory_kernel);
//! build the neighbor list information
template <UInt Dim, typename Container>
void buildNeighbors(Real rcut, bool *pbc, Container &cont);
//! function to know if everything is ready for filtering
bool isInitialized() { return (neighbors != NULL); };
private:
//! build a spatially filtered version of the field provided
template <UInt Dim, bool pbc_x, bool pbc_y, bool pbc_z, typename Container>
void buildFilteredField(std::vector<Real> &filtered, FieldType ftype,
Container &container);
//! build a spatially filtered version of the field provided
template <UInt Dim, FieldType ftype, bool pbc_x, bool pbc_y, bool pbc_z,
typename Container>
void buildFilteredField(std::vector<Real> &filtered, Container &container);
template <UInt Dim, typename FuncWindow> void applyWindow();
template <UInt Dim>
void getReplicatedPosition(Real *replica, Real *pos, int *n);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
//! get the memory kernel
std::vector<Real> &getMemoryKernel();
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
//! the memory kernel content
std::vector<Real> memory_kernel;
// kernel length
UInt kernel_length;
// length of finite impulse response array in terms of lattice^2
double kernel_distance, kernel_distance_sq;
// double kernel_distance_first_neigh;
//! verbose flag
// UInt verbose;
//! frequency to write to file
// UInt freq;
//! cubic spline object
CubicSpline *cSpline;
//! domain bounding box
Cube domain_bounding_box;
//! pointer to neighbors information
ContainerNeighborAtomsInterface *neighbors;
//! the integral density
Real density;
//! cutoff radius to take into account the
Real cutoff2;
//! Number of periodic replica in each direction
int nReplica[3]; // dont change to UInt
};
/* -------------------------------------------------------------------------- */
/* inline functions */
/* -------------------------------------------------------------------------- */
inline std::vector<Real> &SpatialFilterAtomic::getMemoryKernel() {
return memory_kernel;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim, FieldType ftype, bool pbc_x, bool pbc_y, bool pbc_z,
typename Container>
inline void SpatialFilterAtomic::buildFilteredField(std::vector<Real> &filtered,
Container &container) {
typedef typename Container::Ref Ref;
ContainerNeighborAtoms<Ref, Dim> *neigh;
neighbors->castPointer(neigh);
filtered.resize(container.size() * Dim);
UInt index_atom = 0, index = 0;
Real wt = 0.0, d2 = 0.0, dist = 0.0;
typename Container::iterator it = container.begin();
for (Ref at = it.getFirst(); !it.end(); at = it.getNext()) {
index = index_atom * Dim;
wt = density * memory_kernel[0];
for (UInt i = 0; i < Dim; ++i) {
filtered[index + i] = at.template field<ftype>(i) * wt;
}
Real pos[Dim];
at.getPositions0(pos);
// replicas of itself
for (int nx = -nReplica[0]; nx < nReplica[0] + 1; ++nx) {
for (int ny = -nReplica[1]; ny < nReplica[1] + 1; ++ny) {
for (int nz = -nReplica[2]; nz < nReplica[2] + 1; ++nz) {
int n[3] = {nx, ny, nz};
Real replica[Dim];
getReplicatedPosition<Dim>(replica, pos, n);
d2 = at.template computeInitialDistance<pbc_x, pbc_y, pbc_z>(
replica, domain_bounding_box);
if (d2 > kernel_distance_sq || d2 == 0.0) {
continue;
} // rejected
// if (d2 < kernel_distance_sq && d2 != 0.0)
dist = sqrt(d2);
// spline interpolation
cSpline->fastEvaluateSpline(dist, wt);
for (UInt i = 0; i < Dim; ++i) {
filtered[index + i] += at.template field<ftype>(i) * wt;
}
}
}
}
// replicas of neighbor atoms
typename ContainerNeighborAtoms<Ref, Dim>::iterator it = neigh->begin();
Ref at_neighbor;
for (at_neighbor = it.getFirst0(at); !it.end();
at_neighbor = it.getNext()) {
Real pos_neighbor[Dim];
at_neighbor.getPositions0(pos_neighbor);
for (int nx = -nReplica[0]; nx < nReplica[0] + 1; ++nx) {
for (int ny = -nReplica[1]; ny < nReplica[1] + 1; ++ny) {
for (int nz = -nReplica[2]; nz < nReplica[2] + 1; ++nz) {
int n[3] = {nx, ny, nz};
Real replica[Dim];
getReplicatedPosition<Dim>(replica, pos_neighbor, n);
d2 = at.template computeInitialDistance<pbc_x, pbc_y, pbc_z>(
replica, domain_bounding_box);
if (d2 > kernel_distance_sq || d2 == 0.0) {
continue;
} // rejected
dist = sqrt(d2);
// spline interpolation
cSpline->fastEvaluateSpline(dist, wt);
wt *= density;
for (UInt i = 0; i < Dim; ++i) {
filtered[index + i] += at_neighbor.template field<ftype>(i) * wt;
}
}
}
}
}
++index_atom;
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim, typename Container>
void SpatialFilterAtomic::buildNeighbors(Real rcut, bool *pbc,
Container &cont) {
domain_bounding_box = cont.getBoundingBox();
if (Dim != domain_bounding_box.getDim())
LM_FATAL("mismatched dimensions " << Dim << " "
<< domain_bounding_box.getDim());
if (neighbors)
delete neighbors;
typedef typename Container::Ref Ref;
ContainerNeighborAtoms<Ref, Dim> *tmp;
if (nReplica[0] != 0 || nReplica[1] != 0 || nReplica[0] != 0) {
Real *dsize = domain_bounding_box.getSize();
UInt ncells[3] = {1, 1, 1};
for (UInt i = 0; i < Dim; ++i) {
int temp = (int)(dsize[i] / rcut);
if (temp > 0)
ncells[i] = temp;
if (nReplica[i] && temp > 0)
DUMP("Direction " << i << " domain size: " << dsize[i]
<< " rcut: " << rcut
<< "Use of replica is not advised in this case."
<< "This may lead to wrong results!",
DBG_WARNING);
// force pbc to false
pbc[i] = false;
}
tmp = new ContainerNeighborAtoms<Ref, Dim>(domain_bounding_box, ncells);
tmp->setPBCFlag(pbc);
tmp->fillGridWithInitialPositions(cont);
neighbors = tmp;
} else {
tmp = new ContainerNeighborAtoms<Ref, Dim>(domain_bounding_box, rcut);
tmp->setPBCFlag(pbc);
// tmp->fillGrid(cont);
tmp->fillGridWithInitialPositions(cont);
neighbors = tmp;
// force replica to zero (may change in future)
for (UInt i = 0; i < Dim; ++i)
nReplica[i] = 0;
}
cutoff2 = rcut * rcut;
}
/* -------------------------------------------------------------------------- */
template <UInt Dim, typename Container>
void SpatialFilterAtomic::buildFilteredField(std::vector<Real> &filtered,
FieldType ftype, bool *pbc_flag,
Container &container) {
UInt pbc[3] = {0, 0, 0};
for (UInt i = 0; i < Dim; ++i)
pbc[i] = pbc_flag[i];
if (pbc[0] == 0 && pbc[1] == 0 && pbc[2] == 0)
buildFilteredField<Dim, false, false, false>(filtered, ftype, container);
if (pbc[0] == 1 && pbc[1] == 0 && pbc[2] == 0)
buildFilteredField<Dim, true, false, false>(filtered, ftype, container);
if (pbc[0] == 1 && pbc[1] == 1 && pbc[2] == 0)
buildFilteredField<Dim, true, true, false>(filtered, ftype, container);
if (pbc[0] == 1 && pbc[1] == 1 && pbc[2] == 1)
buildFilteredField<Dim, true, true, true>(filtered, ftype, container);
if (pbc[0] == 1 && pbc[1] == 0 && pbc[2] == 1)
buildFilteredField<Dim, true, false, true>(filtered, ftype, container);
if (pbc[0] == 0 && pbc[1] == 1 && pbc[2] == 0)
buildFilteredField<Dim, false, true, false>(filtered, ftype, container);
if (pbc[0] == 0 && pbc[1] == 1 && pbc[2] == 1)
buildFilteredField<Dim, false, true, true>(filtered, ftype, container);
if (pbc[0] == 0 && pbc[1] == 0 && pbc[2] == 1)
buildFilteredField<Dim, false, false, true>(filtered, ftype, container);
}
/* -------------------------------------------------------------------------- */
//! build a spatially filtered version of the field provided
template <UInt Dim, bool pbc_x, bool pbc_y, bool pbc_z, typename Container>
void SpatialFilterAtomic::buildFilteredField(std::vector<Real> &filtered,
FieldType ftype,
Container &container) {
switch (ftype) {
case _displacement:
buildFilteredField<Dim, _displacement, pbc_x, pbc_y, pbc_z>(filtered,
container);
break;
case _velocity:
buildFilteredField<Dim, _velocity, pbc_x, pbc_y, pbc_z>(filtered,
container);
break;
case _force:
buildFilteredField<Dim, _force, pbc_x, pbc_y, pbc_z>(filtered, container);
break;
case _stress:
buildFilteredField<Dim, _stress, pbc_x, pbc_y, pbc_z>(filtered, container);
break;
default:
LM_FATAL("this field cannot be filtered in space");
}
}
/* -------------------------------------------------------------------------- */
template <UInt Dim, typename F>
void SpatialFilterAtomic::buildKernelFromFunction(F &foo, Real xmin, Real xmax,
UInt n_points,
WindowFunction window) {
std::vector<Real> x;
Real pos = xmin;
Real step = (xmax - xmin) / n_points;
for (UInt i = 0; i < n_points; ++i, pos += step) {
x.push_back(pos);
Real res = foo.template compute<0>(pos);
memory_kernel.push_back(res);
}
switch (window) {
case Hanning:
applyWindow<Dim, HanningWindow>();
break;
case Hamming:
applyWindow<Dim, HammingWindow>();
break;
case Rectangular:
break;
}
buildCSpline(x, memory_kernel);
}
/* -------------------------------------------------------------------------- */
template <UInt Dim, typename FuncWindow>
void SpatialFilterAtomic::applyWindow() {
UInt size = memory_kernel.size();
FuncWindow f(size - 1);
for (UInt i = 0; i < size; ++i) {
for (UInt j = 0; j < Dim; ++j)
memory_kernel[i] *= f.template compute<0>(i);
}
}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template <UInt Dim>
void SpatialFilterAtomic::getReplicatedPosition(Real *replica, Real *pos,
int *n) {
for (UInt i = 0; i < Dim; ++i) {
replica[i] = pos[i] + n[i] * domain_bounding_box.getSize(i);
}
}
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_SPATIAL_FILTER_ATOMIC_HH__ */
c4science · Help