Page MenuHomec4science
Files F63702393

system_base.cc
No OneTemporary
Actions

File Metadata

Created
Tue, May 21, 23:08

system_base.cc
View Options

/**
* file system_base.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 09 May 2017
*
* @brief Implementation for system base class
*
* @section LICENCE
*
* Copyright (C) 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre 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 License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include <stdexcept>
#include <algorithm>
#include <utility>
#include "system/system_base.hh"
#include "iostream"
namespace muSpectre {
namespace Projection {
//----------------------------------------------------------------------------//
//! computes fft frequecies. see function fftfreq in numpy's helper.py
vector fft_freqs(Uint nb_samples) {
vector retval(1, nb_samples);
int N = (nb_samples-1)/2 + 1; // needs to be signed int for neg freqs
for (int i = 0; i < N; ++i) {
retval(0, i) = i;
}
for (int i = N; i < static_cast<int>(nb_samples); ++i) {
retval(0, i) = -(static_cast<int>(nb_samples)/2)+i-N;
}
return retval;
}
//----------------------------------------------------------------------------//
vector fft_freqs(Uint nb_samples, Real length) {
return fft_freqs(nb_samples)/length;
}
//----------------------------------------------------------------------------//
template <Dim_t Dim>
FreqStruc<Dim>::FreqStruc(const std::array<Real, Dim>& sizes,
const Index_t<Dim>& nb_pixels)
:nb_pixels(nb_pixels), sizes(sizes)
{
for (size_t i = 0; i < Dim; ++i) {
this->xi[i] = Projection::fft_freqs(this->nb_pixels[i]);
}
}
//----------------------------------------------------------------------------//
template <Dim_t Dim>
inline typename FreqStruc<Dim>::iterator FreqStruc<Dim>::begin() {
return iterator(*this, {0});
}
//----------------------------------------------------------------------------//
template <Dim_t Dim>
inline typename FreqStruc<Dim>::iterator FreqStruc<Dim>::end() {
return iterator(*this, this->nb_pixels);
}
//----------------------------------------------------------------------------//
template<Dim_t Dim>
FreqStruc<Dim>::iterator::iterator(const FreqStruc<Dim> & parent,
Index_t<Dim> index)
:parent(parent),
index(std::accumulate(index.begin(), index.end(), 1, std::multiplies<size_t>())){ }
//----------------------------------------------------------------------------//
template<Dim_t Dim>
inline Index_t<Dim> FreqStruc<Dim>::iterator::get_index() {
// implements row-major choice. might have to be templated at some point
static_assert(((Dim == twoD) || (Dim == threeD)), "only 2d or 3d");
auto & stride = this->parent.nb_pixels[1];
return Index_t<Dim>{static_cast<Dim_t>(this->index/stride),
static_cast<Dim_t>(this->index%stride)};
}
//----------------------------------------------------------------------------//
template<>
inline Index_t<threeD> FreqStruc<threeD>::iterator::get_index() {
// implements row-major choice. might have to be templated at some point
auto && kstride = this->parent.nb_pixels[2];
Dim_t && k = this->index%kstride;
auto && jstride = this->parent.nb_pixels[1];
Dim_t && j = (this->index/kstride)%jstride;
Dim_t && i = (this->index/kstride)/jstride;
return Index_t<threeD>{i, j, k};
}
//----------------------------------------------------------------------------//
template<Dim_t Dim>
inline Eigen::Matrix<Real, Dim, 1> FreqStruc<Dim>::iterator::operator*() {
Eigen::Matrix<Real, Dim, 1> retval;
auto && indices = this->get_index();
if (Dim == twoD) {
retval << parent.xi[0](indices[0]), parent.xi[1](indices[1]);
} else if (Dim == threeD) {
retval <<
parent.xi[0](indices[0]),
parent.xi[1](indices[1]),
parent.xi[2](indices[2]);
}
return retval;
}
//----------------------------------------------------------------------------//
template <Dim_t Dim>
inline typename FreqStruc<Dim>::iterator& FreqStruc<Dim>::iterator::operator++() {
this->index++;
return *this;
}
} // Projection
//----------------------------------------------------------------------------//
template<Dim_t DimS, Dim_t DimM>
SystemBase<DimS, DimM>::SystemBase(std::array<Real, DimS> sizes, Index nb_pixels,
Formulation form)
: sizes{sizes}, nb_pixels{nb_pixels},
nb_pixel(std::accumulate(nb_pixels.begin(), nb_pixels.end(), 1, std::multiplies<size_t>())),
form{form}
{
for (const auto & size: this->nb_pixels) {
if (size%2 != 1) {
throw std::runtime_error("Sorry, for now can only handle odd grid sizes");
}
}
if (form != Formulation::finite_strain) {
throw std::runtime_error("Sorry, only finite strain for the time being");
}
//resize and allocate
this->Ghats = T4Vector(DimM, DimM, DimM, DimM, this->size());
this->Ghats.setZero();
this->build_ghats();
}
//----------------------------------------------------------------------------//
template <Dim_t DimS, Dim_t DimM>
void SystemBase<DimS, DimM>::build_ghats() {
//! row-major storage for pixels
Projection::FreqStruc<DimS> freqs(this->sizes, this->nb_pixels);
auto it = freqs.begin();
//drop the first frequency, as it is zero
++it;
for (size_t pix_id = 1; pix_id < this->nb_pixel; ++pix_id, ++it) {
auto xi = *it;
xi/=xi.norm();
using T4shape = Eigen::Sizes<DimS, DimS, DimS, DimS>;
using T4 = Eigen::TensorFixedSize<Real, T4shape>;
Eigen::TensorMap<T4> Ghat_kk(&this->Ghats(0, 0, 0, 0, pix_id),
DimS, DimS, DimS, DimS);
for (Dim_t im = 0; im < DimS; ++im) {
for (Dim_t j = 0; j < DimS; ++j) {
for (Dim_t l = 0; l < DimS; ++l) {
Ghat_kk(im, j, l, im) = xi(j)*xi(l);
}
}
}
}
}
template class SystemBase<2, 2>;
template class SystemBase<2, 3>;
template class SystemBase<3, 3>;
} // muSpectre

Event Timeline

c4science · Help