Homec4science
Diffusion tamaas (features/declarative) Browse

tamaas/src/model/integration/accumulator.hh
383b32a7e7befeatures/declarative

Last Change
Raw File

accumulator.hh
View Options

/**
* @file
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2018 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Tamaas 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.
*
* Tamaas 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 Tamaas. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef ACCUMULATOR_HH
#define ACCUMULATOR_HH
/* -------------------------------------------------------------------------- */
#include "grid_hermitian.hh"
#include "integrator.hh"
#include "model_type.hh"
#include "static_types.hh"
#include <array>
#include <vector>
/* -------------------------------------------------------------------------- */
namespace tamaas {
/// Accumulation integral manager
template <model_type type, typename Source,
typename = std::enable_if_t<is_proxy_v<Source>>>
class Accumulator {
using trait = model_type_traits<type>;
using BufferType = GridHermitian<Real, trait::boundary_dimension>;
public:
/// Constructor
Accumulator() {
// Resizing accumulators (with source)
for (auto&& acc : accumulator)
acc.setNbComponents(Source::size);
}
/// Initialize uniform mesh
void makeUniformMesh(UInt N, Real domain_size) {
Real dx = domain_size / (N - 1);
Real x = 0;
node_positions.resize(N);
for (UInt i = 0; i < N; ++i, x += dx) {
node_positions[i] = x;
}
}
const std::vector<Real>& nodePositions() const {
return node_positions;
}
/// Direction flag
enum class direction { forward, backward };
/// direction <=> upper
static constexpr bool is_upper(direction dir) {
return dir == direction::backward;
}
/// Forward/Backward iterator for integration passes
template <direction dir>
struct iterator {
using integ = Integrator<1>;
/// Constructor
iterator(Accumulator& acc, Int k) : acc(acc), k(k), l(k) {
l = (dir != direction::forward) ? k + 1 : k;
}
/// Copy constructor
iterator(const iterator& o)
: acc(o.acc), k(o.k), l(o.l), r(o.r), xc(o.xc) {}
/// Compare
bool operator!=(const iterator& o) const { return k != o.k; }
/// Increment
iterator& operator++() {
// update current layer and element. skip last operation
if (not setup()) {
next();
return *this;
}
// part of the domain to be integrated
constexpr bool upper = is_upper(dir);
Loop::loop(
[&](auto qv, auto wk0, auto wk1, auto acc_g0, auto acc_g1) {
Real q = qv.l2norm();
acc_g0 += wk0 * integ::template G0<upper, 0>(q, r, xc);
acc_g0 += wk1 * integ::template G0<upper, 1>(q, r, xc);
acc_g1 += wk0 * integ::template G1<upper, 0>(q, r, xc);
acc_g1 += wk1 * integ::template G1<upper, 1>(q, r, xc);
},
range<VectorProxy<const Real, trait::dimension - 1>>(acc.waveVectors()),
range<Source>(acc.nodeValues()[k]),
range<Source>(acc.nodeValues()[k + 1]),
range<Source>(acc.accumulator.front()),
range<Source>(acc.accumulator.back()));
next();
return *this;
}
/// Dereference iterator
std::tuple<int, Real, BufferType&, BufferType&> operator*() {
return {l, acc.node_positions[l], acc.accumulator.front(),
acc.accumulator.back()};
}
protected:
/// Update index layer and element info
bool setup() {
if (k + 1 >= static_cast<int>(acc.node_positions.size()) or k < 0)
return false;
r = 0.5 * (acc.node_positions[k + 1] - acc.node_positions[k]);
xc = 0.5 * (acc.node_positions[k + 1] + acc.node_positions[k]);
return true;
}
void next() {
l = (dir == direction::forward) ? k + 1 : k;
k += (dir == direction::forward) ? 1 : -1;
}
protected:
Accumulator& acc; /// accumulator holder
Int k = 0; ///< element index
Int l = 0; ///< layer index
Real r = 0; ///< element radius
Real xc = 0; ///< element center
};
/// Forward/Backward iterator for integration of fundamental mode
template <direction dir>
struct fundamental_iterator : iterator<dir> {
using iterator<dir>::iterator;
using integ = typename iterator<dir>::integ;
fundamental_iterator& operator++() {
this->next();
//TODO fix
Source wk0(&this->acc.nodeValues()[this->k](0)),
wk1(&this->acc.nodeValues()[this->k + 1](0));
acc_tensor += wk0 * integ::template F<0>(this->r);
acc_tensor += wk1 * integ::template F<1>(this->r);
return *this;
}
std::tuple<int, Real, Source> operator*() {
return {this->l, this->acc.node_positions[this->l], Source(acc_tensor)};
}
private:
typename Source::stack_type acc_tensor;
};
/// Range for convenience
template <typename Iterator>
struct acc_range {
Iterator _begin, _end;
Iterator begin() const { return _begin; }
Iterator end() const { return _end; }
};
/// Prepare forward loop
acc_range<iterator<direction::forward>>
forward(std::vector<BufferType>& nvalues,
const Grid<Real, trait::boundary_dimension>& wvectors) {
using it = iterator<direction::forward>;
reset(nvalues, wvectors);
return acc_range<it>{it(*this, 0), it(*this, node_positions.size())};
}
/// Prepare backward loop
acc_range<iterator<direction::backward>>
backward(std::vector<BufferType>& nvalues,
const Grid<Real, trait::boundary_dimension>& wvectors) {
using it = iterator<direction::backward>;
reset(nvalues, wvectors);
// don't ask why the range has those values
// if it works, don't change it
return acc_range<it>{it(*this, node_positions.size() - 2), it(*this, -2)};
}
/// Forward loop for fundamental mode
acc_range<fundamental_iterator<direction::forward>>
fundamentalForward(std::vector<BufferType>& nvalues) {
using it = fundamental_iterator<direction::forward>;
reset(nvalues, *this->wavevectors);
return acc_range<it>{it(*this, 0), it(*this, node_positions.size())};
}
/// Backward loop for fundamental mode
acc_range<fundamental_iterator<direction::backward>>
fundamentalBackward(std::vector<BufferType>& nvalues) {
using it = fundamental_iterator<direction::backward>;
reset(nvalues, *this->wavevectors);
return acc_range<it>{it(*this, node_positions.size() - 1), it(*this, -1)};
}
void reset(std::vector<BufferType>& nvalues,
const Grid<Real, trait::boundary_dimension>& wvectors) {
node_values = &nvalues;
wavevectors = &wvectors;
for (auto&& acc : accumulator) {
acc.resize(wvectors.sizes());
acc = 0;
}
}
std::vector<BufferType>& nodeValues() {
TAMAAS_ASSERT(node_values, "Node values is invalid");
return *node_values;
}
const Grid<Real, trait::boundary_dimension>& waveVectors() {
TAMAAS_ASSERT(wavevectors, "Wavevectors is invalid");
return *wavevectors;
}
private:
std::array<BufferType, 2> accumulator;
std::vector<Real> node_positions;
std::vector<BufferType>* node_values;
const Grid<Real, trait::boundary_dimension>* wavevectors;
};
} // namespace tamaas
#endif // ACCUMULATOR_HH
c4science · Help