Homec4science
Diffusion LibMultiScale (features/pybind) Browse

libmultiscale/src/math/dispersion_relation.hh
b05d5f2aba78features/pybind

Last Change
Raw File

dispersion_relation.hh
View Options

/**
* @file dispersion_relation.hh
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
* @author Srinivasa Babu Ramisetti <srinivasa.ramisetti@epfl.ch>
*
* @date Mon Oct 28 19:23:14 2013
*
* @brief This function implements the dispersion relation of various crytals
*
* @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_DISPERSION_RELATION_HH__
#define __LIBMULTISCALE_DISPERSION_RELATION_HH__
/* -------------------------------------------------------------------------- */
#include "function_interface.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
/** Class DispersionRelation
*
* This class implements access to the dispersion relation
*
*/
class DispersionRelationInterface : public FunctionInterface
{
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
DispersionRelationInterface(){
this->elastic_stiffness = 1.0;
this->mass = 1.0;
this->r0 = 1.0;
};
DispersionRelationInterface(Real stiffness,
const Quantity<Mass> & mass,
const Quantity<Length> & r0){
this->elastic_stiffness = stiffness;
DUMP(" STIFFNESS " << this->elastic_stiffness, DBG_MESSAGE);
this->mass = mass;
DUMP(" ATOM MASS " << this->mass, DBG_MESSAGE);
this->r0 = r0;
DUMP(" R0 " << this->r0, DBG_MESSAGE);
};
virtual ~DispersionRelationInterface(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
// inline Real compute(Real * k, UInt branch);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
Real elastic_stiffness;
Quantity<Length> r0;
Quantity<Mass> mass;
};
/* -------------------------------------------------------------------------- */
class DispersionRelation : public DispersionRelationInterface
{
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
DispersionRelation(){};
DispersionRelation(Real stiffness,
const Quantity<Mass> & mass,
const Quantity<Length> & r0)
: DispersionRelationInterface(stiffness,mass,r0){};
virtual ~DispersionRelation(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
template <LatticeType lattice_type>
inline Real compute(Real * k, UInt branch);
inline std::vector<Real> compute(std::vector<Real> & wave_numbers,
const LatticeType & type, UInt branch);
template <LatticeType lattice_type>
inline std::vector<Real> compute(std::vector<Real> & wave_numbers, UInt branch);
template <LatticeType lattice_type>
inline void computePolarization(Real * k, Real * pol, UInt branch);
inline std::vector<Real> computePolarization(std::vector<Real> & wave_numbers,
const LatticeType & type,
UInt branch);
template <LatticeType lattice_type>
inline std::vector<Real> computePolarization(std::vector<Real> & wave_numbers,
UInt branch);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
private:
};
/* -------------------------------------------------------------------------- */
template<LatticeType lattice_type>
inline Real DispersionRelation::compute(Real * k, UInt branch){
LM_TOIMPLEMENT;
}
/* -------------------------------------------------------------------------- */
template<LatticeType lattice_type>
inline void DispersionRelation::computePolarization(Real * k, Real * pol, UInt branch){
static const UInt Dim = LatticeDim<lattice_type>::value;
pol[0] = 1.;
for (UInt i = 1; i < Dim; ++i) pol[i] = 0.;
}
/* -------------------------------------------------------------------------- */
template<>
inline void DispersionRelation::computePolarization<hex_2d>(Real * k, Real * pol, UInt branch){
Real A = cos(k[0]*r0/2)*cos(k[1]*sqrt(3)*r0/2);
Real B = sin(k[0]*r0/2)*sin(sqrt(3)*k[1]*r0/2);
Real C = cos(k[0]*r0);
Real D = C-A;
Real D2 = D*D;
Real B2 = B*B;
Real sign = pow(-1.,branch);
DUMP("A = " << A,DBG_DETAIL);
DUMP("B = " << B,DBG_DETAIL);
DUMP("C = " << C,DBG_DETAIL);
DUMP("D = " << D,DBG_DETAIL);
DUMP("sign = " << sign,DBG_DETAIL);
if (B != 0){
Real Delta = D2+3.*B2;
DUMP("Delta = " << Delta,DBG_DETAIL);
pol[0] = A - C - sqrt(Delta);
pol[1] = sqrt(3.)* B;
Real n = sqrt(pol[0]*pol[0] + pol[1]*pol[1]);
DUMP("n = " << n,DBG_DETAIL);
pol[0] /= n;
pol[1] /= n;
}
else {
if (sign < 0){
pol[0] = 0;
pol[1] = 1;
}
else{
pol[0] = 1;
pol[1] = 0;
}
}
DUMP("k = " << k[0] << " " << k[1]
<< " polarization = " << pol[0] << " " << pol[1],DBG_MESSAGE);
}
/* -------------------------------------------------------------------------- */
template <>
inline Real DispersionRelation::compute<mono_atom_1d>(Real * k, UInt branch){
Real omega=0.0;
omega = 2.0
* sqrt(elastic_stiffness*code_unit_system.e_dd2m_tt/mass)
*(std::abs(sin(k[0]/2.0*r0)));
DUMP("k = " << k[0] << " omega = " << omega,DBG_MESSAGE);
return omega;
}
/* -------------------------------------------------------------------------- */
template <>
inline Real DispersionRelation::compute<hex_2d>(Real * k,UInt branch){
Real A = cos(k[0]*r0/2)*cos(k[1]*sqrt(3)*r0/2);
Real B = sin(k[0]*r0/2)*sin(sqrt(3)*k[1]*r0/2);
Real C = cos(k[0]*r0);
Real D = C-A;
Real D2 = D*D;
Real B2 = B*B;
Real Delta = D2+3.*B2;
Real sign = pow(-1.,branch);
Real omega =
sqrt(
(elastic_stiffness*code_unit_system.e_dd2m_tt/mass)
*
(3.0 - 2.0*A - C + sign*sqrt(Delta))
);
// Real omega1D = compute<mono_atom_1d>(k,branch);
DUMP("k = " << k[0] << " " << k[1] << " omega = " << omega,DBG_MESSAGE);
// DUMP("k = " << k[0] << " " << k[1] << " omega1D = " << omega1D,DBG_MESSAGE);
return omega;
}
/* -------------------------------------------------------------------------- */
template <>
inline Real DispersionRelation::compute<hex_2d_radial>(Real * k, UInt branch){
Real temp0=(cos(2*k[0])-cos(k[0])*cos(k[1]));
Real temp1=3.0*sin(k[0])*sin(k[0])*sin(k[1])*sin(k[1]);
Real omega =
sqrt(
(elastic_stiffness/mass)
*
(3.0-cos(2*k[0])-2.0*cos(k[0])*cos(k[1])-sqrt(temp0*temp0+temp1))
);
return omega;
}
/* -------------------------------------------------------------------------- */
inline std::vector<Real> DispersionRelation::compute(std::vector<Real> & wave_numbers,
const LatticeType & type,
UInt branch){
switch (type){
case mono_atom_1d :return compute<mono_atom_1d>(wave_numbers,branch); break;
case hex_2d :return compute<hex_2d>(wave_numbers,branch); break;
case hex_2d_radial :return compute<hex_2d_radial>(wave_numbers,branch); break;
case fcc_3d :return compute<fcc_3d>(wave_numbers,branch); break;
default: LM_FATAL("internal error: unknown lattice");
}
}
/* -------------------------------------------------------------------------- */
inline std::vector<Real> DispersionRelation::computePolarization(std::vector<Real> & wave_numbers,
const LatticeType & type,
UInt branch){
switch (type){
case mono_atom_1d :return computePolarization<mono_atom_1d>(wave_numbers,branch); break;
case hex_2d :return computePolarization<hex_2d>(wave_numbers,branch); break;
case hex_2d_radial :return computePolarization<hex_2d_radial>(wave_numbers,branch); break;
case fcc_3d :return computePolarization<fcc_3d>(wave_numbers,branch); break;
default: LM_FATAL("internal error: unknown lattice");
}
}
/* -------------------------------------------------------------------------- */
template <LatticeType type>
inline std::vector<Real> DispersionRelation::compute(std::vector<Real> & wave_numbers,
UInt branch){
static const UInt Dim = LatticeDim<type>::value;
UInt total_nb_modes = wave_numbers.size()/Dim;
std::vector<Real> res;
for (UInt nb = 0; nb < total_nb_modes; ++nb) {
Real * k = &wave_numbers[nb*Dim];
res.push_back(this->compute<type>(k,branch));
}
return res;
}
/* -------------------------------------------------------------------------- */
template <LatticeType type>
inline std::vector<Real> DispersionRelation::computePolarization(std::vector<Real> & wave_numbers, UInt branch){
static const UInt Dim = LatticeDim<type>::value;
Real pol[Dim];
UInt total_nb_modes = wave_numbers.size()/Dim;
std::vector<Real> res;
for (UInt nb = 0; nb < total_nb_modes; ++nb) {
Real * k = &wave_numbers[nb*Dim];
this->computePolarization<type>(k,pol,branch);
for (UInt i = 0; i < Dim; ++i) res.push_back(pol[i]);
}
return res;
}
/* -------------------------------------------------------------------------- */
__END_LIBMULTISCALE__
#endif /* __LIBMULTISCALE_DISPERSION_RELATION_HH__ */
c4science · Help