diff --git a/src/core/grid.hh b/src/core/grid.hh
index 4229715..84f1e83 100644
--- a/src/core/grid.hh
+++ b/src/core/grid.hh
@@ -1,234 +1,235 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2022 Lucas Frérot
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#ifndef GRID_HH
#define GRID_HH
/* -------------------------------------------------------------------------- */
#include "grid_base.hh"
#include "tamaas.hh"
#include
#include
#include
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/**
* @brief Multi-dimensional & multi-component array class
*
* This class is a container for multi-component data stored on a multi-
* dimensional grid.
*
* The access function is the parenthesis operator. For a grid of dimension d,
* the operator takes d+1 arguments: the first d arguments are the position on
* the grid and the last one is the component of the stored data.
*
* It is also possible to use the access operator with only one argument, it is
* then considering the grid as a flat array, accessing the given cell of the
* array.
*/
template
class Grid : public GridBase {
template
using is_valid_index = fold_trait;
public:
/* ------------------------------------------------------------------------ */
/* Types */
/* ------------------------------------------------------------------------ */
using value_type = T;
using reference = value_type&;
static constexpr UInt dimension = dim;
/* ------------------------------------------------------------------------ */
/* Constructors */
/* ------------------------------------------------------------------------ */
public:
/// Constructor by default (empty array)
Grid();
+
/// Constructor with shape from iterators
template
Grid(RandomAccessIterator begin, RandomAccessIterator end,
UInt nb_components);
/// Construct with shape from iterators on data view
template
Grid(RandomAccessIterator begin, RandomAccessIterator end, UInt nb_components,
span data);
- /// Constructor
- Grid(const std::array& n, UInt nb_components)
+ /// Constructor with container as shape
+ template
+ Grid(Container&& n, UInt nb_components)
: Grid(std::begin(n), std::end(n), nb_components) {}
- /// Wrap on given data
- Grid(const std::array& n, UInt nb_components,
- span data)
+ /// Constructor with shape and wrapped data
+ template
+ Grid(Container&& n, UInt nb_components, span data)
: Grid(std::begin(n), std::end(n), nb_components, data) {}
- /// Constructor with vectors
- Grid(const std::vector& n, UInt nb_components)
- : Grid(std::begin(n), std::end(n), nb_components) {}
/// Constructor with initializer list
Grid(const std::initializer_list& n, UInt nb_components)
: Grid(std::begin(n), std::end(n), nb_components) {}
+
/// Copy constructor
Grid(const Grid& o) : GridBase(), n(o.n), strides(o.strides) {}
/// Move constructor (transfers data ownership)
Grid(Grid&& o) noexcept
: GridBase(std::move(o)), n(std::move(o.n)),
strides(std::move(o.strides)) {}
+
/// Destructor
~Grid() override = default;
private:
/// Init from standard container
template
void init(const Container& n, UInt nb_components);
public:
/* ------------------------------------------------------------------------ */
/* Common operations */
/* ------------------------------------------------------------------------ */
/// Resize array
void resize(const std::array& n);
/// Resize array (from std::vector)
void resize(const std::vector& n);
/// Resize array (from initializer list)
void resize(std::initializer_list n);
/// Resize array (from iterators)
template
void resize(ForwardIt begin, ForwardIt end);
/// Compute size
inline UInt computeSize() const;
/// Get grid dimension
inline UInt getDimension() const override { return dim; }
/// Compute strides
void computeStrides();
/// Print
virtual void printself(std::ostream& str) const;
/// Get sizes
const std::array& sizes() const { return n; }
/// Get strides
const std::array& getStrides() const { return this->strides; }
/* ------------------------------------------------------------------------ */
/* Access operators (these cannot be moved outside the class) */
/* ------------------------------------------------------------------------ */
/// Variadic access operator (non-const)
template
inline std::enable_if_t::value, T&>
operator()(T1... args) {
/// Checking access integrity
constexpr UInt nargs = sizeof...(T1);
static_assert(nargs == dim + 1 || nargs == 1 || nargs == dim,
"number of arguments in operator() does not match dimension");
constexpr UInt start = (nargs == 1) ? dim : 0;
UInt offset = unpackOffset(0, start, args...);
return this->data[offset];
}
/// Variadic access operator
template
inline std::enable_if_t::value, const T&>
operator()(T1... args) const {
/// Checking access integrity
constexpr UInt nargs = sizeof...(T1);
static_assert(nargs == dim + 1 || nargs == 1 || nargs == dim,
"number of arguments in operator() does not match dimension");
constexpr UInt start = (nargs == 1) ? dim : 0;
UInt offset = unpackOffset(0, start, args...);
return this->data[offset];
}
/// Tuple index access operator
template
inline T& operator()(std::array tuple);
template
inline const T& operator()(std::array tuple) const;
private:
/// Unpacking the arguments of variadic ()
template
inline UInt unpackOffset(UInt offset, UInt index_pos, UInt index,
T1... rest) const;
/// End case for recursion
template
inline UInt unpackOffset(UInt offset, UInt index_pos, UInt index) const;
/// Computing offset for a tuple index
template
inline UInt computeOffset(std::array tuple) const;
/* ------------------------------------------------------------------------ */
/* Move/Copy operators */
/* ------------------------------------------------------------------------ */
public:
using GridBase::operator=;
// = operator
Grid& operator=(const Grid& other);
// = operator (move)
Grid& operator=(Grid&& other) noexcept;
// Copy data from another grid
template
void copy(const Grid& other);
// Move data from another grid
template
void move(Grid&& other) noexcept;
template
void wrap(GridBase& other, Container&& n) {
GridBase::wrap(other);
std::copy(n.begin(), n.end(), this->n.begin());
this->computeStrides();
}
template
void wrap(const GridBase& other, Container&& n) {
GridBase::wrap(other);
std::copy(n.begin(), n.end(), this->n.begin());
this->computeStrides();
}
// Wrap memory (non-const)
void wrap(Grid& other) { wrap(other, other.n); }
// Wrap memory
void wrap(const Grid& other) { wrap(other, other.n); }
/* ------------------------------------------------------------------------ */
/* Member variables */
/* ------------------------------------------------------------------------ */
protected:
std::array n; ///< shape of grid: size per dimension
std::array strides; ///< strides for access
};
} // namespace tamaas
/* -------------------------------------------------------------------------- */
/* Inline/template function definitions */
/* -------------------------------------------------------------------------- */
#include "grid_tmpl.hh"
/* -------------------------------------------------------------------------- */
#endif // GRID_HH
diff --git a/src/model/kelvin.cpp b/src/model/kelvin.cpp
index 2e8c2e7..96cb8a6 100644
--- a/src/model/kelvin.cpp
+++ b/src/model/kelvin.cpp
@@ -1,160 +1,159 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2022 Lucas Frérot
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#include "kelvin.hh"
#include "logger.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
/* Constructors */
/* -------------------------------------------------------------------------- */
template
Kelvin::Kelvin(Model* model) : VolumePotential(model) {
setIntegrationMethod(integration_method::linear, 1e-12);
}
template
void Kelvin::setIntegrationMethod(integration_method method,
Real cutoff) {
this->method = method;
this->cutoff = cutoff;
Logger logger;
if (this->method == integration_method::linear) {
logger.get(LogLevel::debug)
<< TAMAAS_DEBUG_MSG("Setting linear integration method");
this->initialize(dtrait::template source_components,
dtrait::template out_components,
this->model->getDiscretization()[0]);
}
else {
logger.get(LogLevel::debug) << TAMAAS_DEBUG_MSG(
"Setting cutoff integration method (cutoff " << this->cutoff << ')');
this->initialize(dtrait::template source_components,
dtrait::template out_components, 1);
}
auto max_q = Loop::reduce(
[] CUDA_LAMBDA(VectorProxy qv) {
return qv.l2norm();
},
range>(
this->wavevectors));
if (this->method == integration_method::linear and
not std::isfinite(std::exp(max_q * this->model->getSystemSize()[0])))
logger.get(LogLevel::warning)
<< "Probable overflow of integral computation (consider "
"changing integration method to integration_method::cutoff or "
"compiling with real_type='long double')\n";
}
/* -------------------------------------------------------------------------- */
/* Operator implementation */
/* -------------------------------------------------------------------------- */
template
void Kelvin::applyIf(GridBase& source,
GridBase& out,
filter_t pred) const {
KelvinInfluence kelvin(this->model->getShearModulus(),
this->model->getPoissonRatio());
parent::transformSource(source, pred);
// Reset buffer values
for (auto&& layer : this->out_buffer)
layer = 0;
if (method == integration_method::linear) {
linearIntegral(out, kelvin);
} else {
cutoffIntegral(out, kelvin);
}
}
template
void Kelvin::linearIntegral(GridBase& out,
KelvinInfluence& kelvin) const {
detail::KelvinHelper helper;
helper.applyIntegral(this->source_buffer, this->out_buffer, this->wavevectors,
this->model->getSystemSize().front(), kelvin);
helper.applyFreeTerm(this->source_buffer, this->out_buffer, kelvin);
helper.makeFundamentalGreatAgain(this->out_buffer);
parent::transformOutput(
[](auto&& out_buffer, auto layer) ->
typename parent::BufferType& { return out_buffer[layer]; },
out);
}
template
void Kelvin::cutoffIntegral(GridBase& out,
KelvinInfluence& kelvin) const {
detail::KelvinHelper helper;
auto func = [&](auto&& out_buffer,
auto layer) -> typename parent::BufferType& {
auto&& out = out_buffer.front();
out = 0;
helper.applyIntegral(this->source_buffer, out, layer, this->wavevectors,
this->model->getSystemSize().front(), cutoff, kelvin);
helper.applyFreeTerm(this->source_buffer[layer], out, kelvin);
helper.makeFundamentalGreatAgain(out);
return out;
};
parent::transformOutput(func, out);
}
/* -------------------------------------------------------------------------- */
template
std::pair Kelvin::matvecShape() const {
const auto N = (*this->model)["displacement"].getGlobalNbPoints();
return std::make_pair(dtrait::template out_components * N,
dtrait::template source_components * N);
}
/* -------------------------------------------------------------------------- */
template
GridBase Kelvin::matvec(GridBase& X) const {
- const auto dim = trait::dimension;
+ constexpr auto dim = trait::dimension;
const auto shape = this->model->getDiscretization();
- Grid x{shape.begin(), shape.end(),
- dtrait::template source_components, X.view()},
- y{shape, dtrait::template out_components};
+ Grid x{shape, dtrait::template source_components, X.view()};
+ Grid y{shape, dtrait::template out_components};
this->apply(x, y);
- return y;
+ return std::move(y);
}
/* -------------------------------------------------------------------------- */
/* Template instanciation */
/* -------------------------------------------------------------------------- */
template class Kelvin;
template class Kelvin;
template class Kelvin;
/* -------------------------------------------------------------------------- */
} // namespace tamaas
diff --git a/src/model/westergaard.cpp b/src/model/westergaard.cpp
index 30ced2d..6d33600 100644
--- a/src/model/westergaard.cpp
+++ b/src/model/westergaard.cpp
@@ -1,306 +1,305 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2022 Lucas Frérot
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#include "grid_view.hh"
#include "influence.hh"
#include "loop.hh"
#include "model.hh"
#include "static_types.hh"
#include "westergaard.hh"
#include
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
template
Westergaard::Westergaard(Model* model)
: IntegralOperator(model), influence(), buffer(),
engine(FFTEngine::makeEngine()) {
// Copy sizes
auto bdisc = model->getBoundaryDiscretization();
auto boundary_hermitian_sizes =
GridHermitian::hermitianDimensions(
bdisc);
constexpr UInt nb_components = trait::components;
buffer.setNbComponents(nb_components);
buffer.resize(boundary_hermitian_sizes);
influence.setNbComponents(nb_components * nb_components);
influence.resize(boundary_hermitian_sizes);
initInfluence();
}
/* -------------------------------------------------------------------------- */
namespace detail {
template
struct boundary_fft_helper {
template
static void backwardTransform(FFTEngine& e, Buffer&& buffer, Out&& out) {
auto view = make_view(out, 0);
e.backward(view, buffer);
}
};
template
struct boundary_fft_helper {
template
static void backwardTransform(FFTEngine& e, Buffer&& buffer, Out&& out) {
e.backward(out, buffer);
}
};
} // namespace detail
template
template
void Westergaard::fourierApply(Functor func, GridBase& in,
GridBase& out) const try {
auto& i = dynamic_cast&>(in);
auto& full_out = dynamic_cast&>(out);
engine->forward(i, buffer);
/// Applying influence
func(buffer, influence);
/// Backward fourier transform on boundary only
detail::boundary_fft_helper::backwardTransform(*engine, buffer,
full_out);
} catch (const std::bad_cast& e) {
TAMAAS_EXCEPTION("Neumann and dirichlet types do not match model type");
}
/* -------------------------------------------------------------------------- */
template
template
void Westergaard::initFromFunctor(Functor func) {
// Compute wavevectors for influence
auto wavevectors = FFTEngine::template computeFrequencies(
influence.sizes());
// Get boundary physical size
auto system_size = this->model->getBoundarySystemSize();
VectorProxy domain(system_size[0]);
// Normalize wavevectors
wavevectors *= 2 * M_PI;
wavevectors /= domain;
// Apply functor
Loop::loop(func, range>(wavevectors),
range>(influence));
if (mpi::rank() == 0) {
// Set fundamental mode to zero
MatrixProxy mat(influence(0));
mat = 0;
}
}
/* -------------------------------------------------------------------------- */
/// \cond DO_NOT_DOCUMENT
#define NEUMANN_BASIC(type) \
template <> \
void Westergaard::initInfluence() { \
auto E_star = model->getHertzModulus(); \
auto basic = [E_star] CUDA_LAMBDA(VectorProxy q, \
MatrixProxy k) { \
k(0, 0) = 2. / (E_star * q.l2norm()); \
}; \
initFromFunctor(basic); \
}
#define DIRICHLET_BASIC(type) \
template <> \
void Westergaard::initInfluence() { \
auto E_star = model->getHertzModulus(); \
auto basic = [E_star] CUDA_LAMBDA(VectorProxy q, \
MatrixProxy k) { \
k(0, 0) = E_star * q.l2norm() / 2; \
}; \
initFromFunctor(basic); \
}
NEUMANN_BASIC(model_type::basic_1d)
NEUMANN_BASIC(model_type::basic_2d)
DIRICHLET_BASIC(model_type::basic_1d)
DIRICHLET_BASIC(model_type::basic_2d)
#undef NEUMANN_BASIC
#undef DIRICHLET_BASIC
/* -------------------------------------------------------------------------- */
template <>
void Westergaard::initInfluence() {
auto E = model->getYoungModulus();
auto nu = model->getPoissonRatio();
const Complex I(0, 1);
auto surface = [E, nu, I] CUDA_LAMBDA(VectorProxy q_,
MatrixProxy F) {
Real q = q_(0);
q *= 1. / q_.l2norm();
F(0, 0) = 2 * (1 + nu) * (1 - nu * q * q);
F(1, 1) = 2 * (1 - nu * nu);
F(0, 1) = I * q * (1 + nu) * (1 - 2 * nu);
F(1, 0) = -F(0, 1);
F *= 1. / (E * q_.l2norm());
};
initFromFunctor(surface);
}
/* -------------------------------------------------------------------------- */
template <>
void Westergaard::initInfluence() {
auto E = model->getYoungModulus();
auto nu = model->getPoissonRatio();
const Complex I(0, 1);
auto surface = [E, nu, I] CUDA_LAMBDA(VectorProxy q_,
MatrixProxy F) {
Vector q(q_);
q *= 1 / q_.l2norm();
F(0, 0) = 2 * (1 + nu) * (1 - nu * q(0) * q(0));
F(1, 1) = 2 * (1 + nu) * (1 - nu * q(1) * q(1));
F(2, 2) = 2 * (1 - nu * nu);
F(0, 1) = F(1, 0) = -q(0) * q(1) * 2 * nu * (1 + nu);
F(0, 2) = I * q(0) * (1 + nu) * (1. - 2. * nu);
F(1, 2) = I * q(1) * (1 + nu) * (1. - 2. * nu);
F(2, 0) = -F(0, 2);
F(2, 1) = -F(1, 2);
F *= 1. / (E * q_.l2norm());
};
initFromFunctor(surface);
}
/* -------------------------------------------------------------------------- */
template
void Westergaard::initInfluence() {}
/// \endcond
/* ------------------------------------------------------------------------ */
template
void Westergaard::apply(GridBase& input,
GridBase& output) const {
auto apply = [](decltype(buffer)& buffer,
const decltype(influence)& influence) {
Loop::loop(
[] CUDA_LAMBDA(VectorProxy i,
MatrixProxy F) { i = F * i; },
range>(buffer),
range>(influence));
};
fourierApply(apply, input, output);
}
/* -------------------------------------------------------------------------- */
template
Real Westergaard::getOperatorNorm() {
constexpr UInt comp = model_type_traits::components;
Real _norm = 0.0;
_norm = Loop::reduce(
[] CUDA_LAMBDA(MatrixProxy m) {
Real n = thrust::norm(m.l2squared());
return std::isnan(n) ? 0 : n;
},
range>(influence));
auto size = model->getSystemSize();
auto disc = model->getDiscretization();
auto dim = model_type_traits::dimension;
/// TODO: why?
switch (dim) {
case 1:
_norm /= (size[0] / disc[0]) * (size[0] / disc[0]);
break;
default:
for (UInt i = 0; i < dim; i++) {
_norm /= size[i] / disc[i];
}
}
return std::sqrt(_norm);
}
/* -------------------------------------------------------------------------- */
template
std::pair Westergaard::matvecShape() const {
auto n = model->getBoundaryDiscretization();
auto N = std::accumulate(n.begin(), n.end(), comp, std::multiplies<>());
return std::make_pair(N, N);
}
/* -------------------------------------------------------------------------- */
template
GridBase Westergaard::matvec(GridBase& X) const {
auto n = model->getBoundaryDiscretization();
// Creating correct input
- Grid x(n, comp);
- thrust::copy(X.begin(), X.end(), x.begin());
+ Grid x{n, comp, X.view()};
// Creating correct output
if (bdim != dim)
n.insert(n.begin(), 1);
Grid y(n, comp);
apply(x, y);
return std::move(y);
}
/* -------------------------------------------------------------------------- */
/* Template instanciation */
/* -------------------------------------------------------------------------- */
template class Westergaard;
template class Westergaard;
template class Westergaard;
template class Westergaard;
template class Westergaard;
template class Westergaard;
template class Westergaard;
template class Westergaard;
template class Westergaard;
template class Westergaard;
template class Westergaard;
template class Westergaard;
/* -------------------------------------------------------------------------- */
} // namespace tamaas
/* -------------------------------------------------------------------------- */