diff --git a/src/bem/bem_fft_base.cpp b/src/bem/bem_fft_base.cpp
index d385590..de00540 100644
--- a/src/bem/bem_fft_base.cpp
+++ b/src/bem/bem_fft_base.cpp
@@ -1,193 +1,171 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 "bem_fft_base.hh"
/* -------------------------------------------------------------------------- */
#include "bem_meta_functional.hh"
#include "elastic_energy_functional.hh"
-#include "fft_plan_manager.hh"
-#include "surface_statistics.hh"
#include "legacy_types.hh"
+#include "surface_statistics.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
BemFFTBase::BemFFTBase(Surface& p)
: BemInterface(p), functional(new MetaFunctional(*this)),
surface_tractions(p.size(), p.getL()),
surface_displacements(p.size(), p.getL()),
surface_spectral_input(p.size(), p.getL()),
surface_spectral_output(p.size(), p.getL()),
surface_spectral_influence_disp(p.size(), p.getL()),
true_displacements(p.size(), p.getL()),
old_displacements(p.size(), p.getL()), gap(p.size(), p.getL()),
e_star(std::nan("")), tmp_coeff(1. / M_PI) {
// this->setNumberOfThreads(this->nthreads);
max_iterations = 5000;
dump_freq = 100;
this->functional->addFunctionalTerm(new ElasticEnergyFunctional(*this));
std::cerr << this << " is setting up the surfaces";
std::cerr << &p << " has size " << p.size() << std::endl;
}
/* -------------------------------------------------------------------------- */
-BemFFTBase::~BemFFTBase() { FFTPlanManager::get().clean(); }
-
-/* -------------------------------------------------------------------------- */
-
const Surface& BemFFTBase::getTractions() const {
return surface_tractions;
}
/* -------------------------------------------------------------------------- */
const Surface& BemFFTBase::getDisplacements() const {
return surface_displacements;
}
/* -------------------------------------------------------------------------- */
const Surface& BemFFTBase::getSpectralInfluenceOverDisplacement() {
return surface_spectral_influence_disp;
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::setEffectiveModulus(Real e_star) {
this->e_star = e_star;
this->computeSpectralInfluenceOverDisplacement();
}
/* -------------------------------------------------------------------------- */
Real BemFFTBase::getEffectiveModulus() { return this->e_star; }
/* -------------------------------------------------------------------------- */
void BemFFTBase::computeSpectralInfluenceOverDisplacement() {
// if (std::isnan(tmp_coeff)) SURFACE_FATAL("constant tmp_coeff is nan");
// if (std::isnan(e_star)) SURFACE_FATAL("constant e_star is nan");
const Real L = surface.getL();
- auto wavevectors = FFTransform::computeFrequencies(
+ auto wavevectors = FFTEngine::template computeFrequencies(
surface_spectral_influence_disp.sizes());
auto inf_it = surface_spectral_influence_disp.begin(),
inf_begin = surface_spectral_influence_disp.begin();
auto end = wavevectors.end(2);
//#pragma omp parallel for firstprivate(inf_it)
for (auto it = wavevectors.begin(2); it < end; ++it) {
inf_it = inf_begin;
inf_it += it - wavevectors.begin(2);
legacy::VectorProxy q_vec(&(*it), 2);
q_vec *= 2 * M_PI / L;
Real q = q_vec.l2norm();
// Influence function
*inf_it = 2. / (q * e_star);
}
surface_spectral_influence_disp(0) = 0.;
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::computeDisplacementsFromTractions() {
this->applyInfluenceFunctions(surface_tractions, surface_displacements);
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::computeTractionsFromDisplacements() {
this->applyInverseInfluenceFunctions(surface_displacements,
surface_tractions);
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::applyInfluenceFunctions(Surface& input,
Surface& output) {
- FFTransform& plan1 =
- FFTPlanManager::get().createPlan(input, surface_spectral_output);
-
- plan1.forwardTransform();
-
+ this->engine.forward(input, surface_spectral_output);
surface_spectral_output *= surface_spectral_influence_disp;
-
- FFTransform& plan2 =
- FFTPlanManager::get().createPlan(output, surface_spectral_output);
-
- plan2.backwardTransform();
+ this->engine.backward(output, surface_spectral_output);
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::applyInverseInfluenceFunctions(Surface& input,
Surface& output) {
-
- FFTransform& plan1 =
- FFTPlanManager::get().createPlan(input, surface_spectral_output);
-
- plan1.forwardTransform();
-
+ this->engine.forward(input, surface_spectral_output);
surface_spectral_output /= surface_spectral_influence_disp;
surface_spectral_output(0) = 0;
-
- FFTransform& plan2 =
- FFTPlanManager::get().createPlan(output, surface_spectral_output);
-
- plan2.backwardTransform();
+ this->engine.backward(output, surface_spectral_output);
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::computeGaps() {
UInt size = surface.size();
#pragma omp parallel for
for (UInt i = 0; i < size * size; i++) {
gap(i) = true_displacements(i) - surface(i);
}
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::computeTrueDisplacements() {
this->true_displacements = this->surface_displacements;
UInt n = surface.size();
UInt size = n * n;
Real shift = 1e300;
for (UInt i = 0; i < size; ++i) {
if (surface_displacements(i) - shift - surface(i) < 0.) {
shift = surface_displacements(i) - surface(i);
}
}
for (UInt i = 0; i < size; ++i) {
true_displacements(i) = surface_displacements(i) - shift;
}
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::addFunctional(Functional* new_funct) {
this->functional->addFunctionalTerm(new_funct);
}
/* -------------------------------------------------------------------------- */
} // namespace tamaas
diff --git a/src/bem/bem_fft_base.hh b/src/bem/bem_fft_base.hh
index 2ae7a40..f489bb1 100644
--- a/src/bem/bem_fft_base.hh
+++ b/src/bem/bem_fft_base.hh
@@ -1,150 +1,148 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 BEM_FFT_BASE_H
#define BEM_FFT_BASE_H
/* -------------------------------------------------------------------------- */
#include "bem_functional.hh"
#include "bem_interface.hh"
#include "bem_meta_functional.hh"
-#include "fftransform.hh"
#include "surface_statistics.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
class BemFFTBase : public BemInterface {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
BemFFTBase(Surface& p);
- virtual ~BemFFTBase();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
//! compute the displacement from traction input
virtual void computeDisplacementsFromTractions();
//! compute the tractions from displacement input
virtual void computeTractionsFromDisplacements();
//! apply influence funtions
virtual void applyInfluenceFunctions(Surface& input,
Surface& output);
//! apply inverse influence funtions
virtual void applyInverseInfluenceFunctions(Surface& input,
Surface& output);
//! compute the influcence coefficients of pressure over displacement
virtual void computeSpectralInfluenceOverDisplacement();
//! compute the displacements
virtual void computeTrueDisplacements();
//! compute the gap
virtual void computeGaps();
//! get the number of iterations
virtual UInt getNbIterations() const = 0;
//! get the convergence steps of iterations
virtual const std::vector& getConvergenceIterations() const = 0;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
//! return the influcence coefficients
virtual const Surface& getSpectralInfluenceOverDisplacement();
//! return the actual tractions
virtual const Surface& getTractions() const;
//! return the actual displacements
virtual const Surface& getDisplacements() const;
//! get the computed true displacement
const Surface& getTrueDisplacements() { return true_displacements; };
//! get the computed gap
const Surface& getGap() { return gap; };
//! set the elastic effective modulus
void setEffectiveModulus(Real e_star);
//! set the elastic effective modulus
Real getEffectiveModulus();
//! set the maximal number of iterations
void setMaxIterations(UInt max_iter) { this->max_iterations = max_iter; };
//! get the maximal number of iterations
UInt getMaxIterations() { return this->max_iterations; };
//! set dump freq
void setDumpFreq(UInt dump_freq) { this->dump_freq = dump_freq; };
//! get dump freq
UInt getDumpFreq() { return this->dump_freq; };
//! set functional
void addFunctional(::tamaas::Functional* new_funct);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
//! Functional for optimization
MetaFunctional* functional;
//! surface tractions
Surface surface_tractions;
//! surface displacement
Surface surface_displacements;
//! surface input in spectral
SurfaceComplex surface_spectral_input;
//! surface output in spectral
SurfaceComplex surface_spectral_output;
//! surface influence coefficient over displacement
Surface surface_spectral_influence_disp;
//! shifted displacement in order to get the true displacement
Surface true_displacements;
//! old displacements (for stopping criterion)
Surface old_displacements;
//! computing of the gap
Surface gap;
//! Effective young modulus
Real e_star;
//! maximum of iterations
UInt max_iterations;
//! number of iterations
UInt nb_iterations;
//! number of iterations
std::vector convergence_iterations;
//! frequency at which to dump iteration info
UInt dump_freq;
/* ------------------------------------------------------------------------ */
/* temporary trick about some obscure scaling factor */
/* ------------------------------------------------------------------------ */
/*
0 2 = 2*1
1 4 = 2*2
2 10 = 2*5
3 20 = 2*10
4 34 = 2*12
*/
//#define tmp_coeff (2.*n/M_PI)
//#define tmp_coeff (n/M_PI)
const Real tmp_coeff;
};
} // namespace tamaas
#endif /* BEM_FFT_BASE_H */
diff --git a/src/bem/bem_grid.cpp b/src/bem/bem_grid.cpp
index 6c3963c..8bcc155 100644
--- a/src/bem/bem_grid.cpp
+++ b/src/bem/bem_grid.cpp
@@ -1,442 +1,429 @@
/**
*
* @author Lucas Frérot
*
* @section LICENSE
*
* Copyright (©) 2016 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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "bem_grid.hh"
-#include "fft_plan_manager.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
BemGrid::BemGrid(Surface& surface)
: E(1.), nu(0.3), surface(), surface_normals(surface.sizes(), 3),
tractions(surface.sizes(), 3), displacements(surface.sizes(), 3),
gaps(surface.sizes(), 3), true_tractions(nullptr), true_displacements(nullptr),
dump_freq(100) {
this->surface.wrap(surface);
auto hermitian_sizes =
GridHermitian::hermitianDimensions(surface.sizes());
influence.setNbComponents(9);
spectral_tractions.setNbComponents(3);
spectral_displacements.setNbComponents(3);
influence.resize(hermitian_sizes);
spectral_tractions.resize(hermitian_sizes);
spectral_displacements.resize(hermitian_sizes);
// computeSurfaceNormals();
}
/* -------------------------------------------------------------------------- */
BemGrid::~BemGrid() {
- FFTPlanManager::get().clean();
delete true_tractions;
delete true_displacements;
}
/* -------------------------------------------------------------------------- */
const Grid& BemGrid::getTrueTractions() const {
if (true_tractions == nullptr)
TAMAAS_EXCEPTION("True tractions were not allocated");
return *true_tractions;
}
const Grid& BemGrid::getTrueDisplacements() const {
if (true_displacements == nullptr)
TAMAAS_EXCEPTION("True displacements were not allocated");
return *true_displacements;
}
/* -------------------------------------------------------------------------- */
/// This function computes surface slopes in fourier space and fills the normals
/// grid
void BemGrid::computeSurfaceNormals() {
SurfaceComplex spectral_surface(surface.size(), surface.getL());
GridHermitian spectral_slopes(spectral_surface.sizes(), 3);
- FFTransform& plan1 =
- FFTPlanManager::get().createPlan(surface, spectral_surface);
- FFTransform& plan2 =
- FFTPlanManager::get().createPlan(surface_normals, spectral_slopes);
-
- plan1.forwardTransform();
+ engine.forward(surface, spectral_surface);
// Loops should be with signed integers
std::array n;
std::copy(spectral_surface.sizes().begin(), spectral_surface.sizes().end(),
n.begin());
const Complex I(0, 1);
#pragma omp parallel
{
#pragma omp for collapse(2)
for (Int i = 0; i <= n[0] / 2; i++) {
for (Int j = 0; j < n[1]; j++) {
Real q1 = 2 * M_PI * j;
Real q2 = 2 * M_PI * i;
spectral_slopes(i, j, 0) = -q1 * I * spectral_surface(i, j);
spectral_slopes(i, j, 1) = -q2 * I * spectral_surface(i, j);
spectral_slopes(i, j, 2) = 0;
}
}
#pragma omp for collapse(2)
for (Int i = n[0] / 2 + 1; i < n[0]; i++) {
for (Int j = 0; j < n[1]; j++) {
Real q1 = 2 * M_PI * j;
Real q2 = 2 * M_PI * (i - n[0]);
spectral_slopes(i, j, 0) = -q1 * I * spectral_surface(i, j);
spectral_slopes(i, j, 1) = -q2 * I * spectral_surface(i, j);
spectral_slopes(i, j, 2) = 0;
}
}
} // end #pragma omp parallel
// Dirac
spectral_slopes(0, 0, 2) = n[0] * n[0];
- plan2.backwardTransform();
- FFTPlanManager::get().destroyPlan(surface, spectral_surface);
- FFTPlanManager::get().destroyPlan(surface_normals, spectral_slopes);
+ engine.backward(surface_normals, spectral_slopes);
// Make surface normal
legacy::VectorProxy sn(nullptr, 3);
const UInt N = surface_normals.getNbPoints();
#pragma omp parallel for firstprivate(sn)
for (UInt i = 0; i < N; i++) {
sn.setPointer(surface_normals.getInternalData() + i * sn.dataSize());
sn /= sn.l2norm();
}
}
/* -------------------------------------------------------------------------- */
void BemGrid::computeDisplacementsFromTractions() {
applyInfluenceFunctions(tractions, displacements);
}
void BemGrid::computeTractionsFromDisplacements() {
applyInverseInfluenceFunctions(tractions, displacements);
}
/* -------------------------------------------------------------------------- */
void BemGrid::computeGaps(Grid& disp) {
auto n = surface.sizes();
#pragma omp parallel for collapse(2)
for (UInt i = 0; i < n[0]; i++) {
for (UInt j = 0; j < n[1]; j++) {
gaps(i, j, 0) = disp(i, j, 0);
gaps(i, j, 1) = disp(i, j, 1);
gaps(i, j, 2) = disp(i, j, 2) - surface(i, j);
}
}
}
/* -------------------------------------------------------------------------- */
void BemGrid::applyInfluenceFunctions(Grid& input,
Grid& output) {
TAMAAS_ASSERT(input.dataSize() == output.dataSize(),
"input and output have different size");
- FFTransform& plan_1 =
- FFTPlanManager::get().createPlan(input, spectral_tractions);
- FFTransform& plan_2 =
- FFTPlanManager::get().createPlan(output, spectral_displacements);
- plan_1.forwardTransform();
+ engine.forward(input, spectral_tractions);
Complex *F = influence.getInternalData(),
*p = spectral_tractions.getInternalData(),
*u = spectral_displacements.getInternalData();
legacy::VectorProxy spectral_p(nullptr, 3), spectral_u(nullptr, 3);
legacy::MatrixProxy spectral_F(nullptr, 3, 3);
spectral_displacements = 0;
const UInt N = influence.getNbPoints();
// OpenMP not very cool with this: best to have a POD iteration variable
#pragma omp parallel for firstprivate(spectral_p, spectral_u, spectral_F)
for (UInt i = 0; i < N; ++i) {
spectral_p.setPointer(p + i * spectral_p.dataSize());
spectral_u.setPointer(u + i * spectral_u.dataSize());
spectral_F.setPointer(F + i * spectral_F.dataSize());
spectral_u.mul(spectral_F, spectral_p);
}
- plan_2.backwardTransform();
+ engine.backward(output, spectral_displacements);
}
/* -------------------------------------------------------------------------- */
void BemGrid::applyInverseInfluenceFunctions(
Grid& input[[gnu::unused]], Grid& output[[gnu::unused]]) {
TAMAAS_EXCEPTION("applyInverseInfluenceFunctions not yet implemented");
}
/* -------------------------------------------------------------------------- */
// warning: this function will suck your soul out
// loops must be written with signed integers
void BemGrid::computeInfluence() {
// Iterator over points of surface
legacy::MatrixProxy it(influence.getInternalData(), 3, 3);
// Imaginary unit
const Complex I(0, 1);
const Int size = this->surface.size();
const Real L_1 = 1.; // this->surface.lengths()[0];
const Real L_2 = 1.; // this->surface.lengths()[1];
this->influence = 0;
// Influence functions matrix, ask Lucas
#define INFLUENCE(q_1, q_2, q, E, nu) \
{ \
2. / (E * q * q * q) * (-nu * nu * q_1 * q_1 - nu * q_2 * q_2 + q * q), \
-q_1 *q_2 / (E * q * q * q) * (1 + nu) * 2 * nu, \
I *q_1 / (E * q * q) * (1 + nu) * (1 - 2 * nu), \
-q_1 *q_2 / (E * q * q * q) * (1 + nu) * 2 * nu, \
2. / (E * q * q * q) * \
(-nu * nu * q_2 * q_2 - nu * q_1 * q_1 + q * q), \
I *q_2 / (E * q * q) * (1 + nu) * (1 - 2 * nu), \
-I *q_1 / (E * q * q) * (1 + nu) * (1 - 2 * nu), \
-I *q_2 / (E * q * q) * (1 + nu) * (1 - 2 * nu), \
2. / (E * q) * (1 - nu * nu) \
}
// Fill influence for single point
#define FILL_INFLUENCE(k, l) \
do { \
Real q_1 = 2 * M_PI * (l) / L_1; \
Real q_2 = 2 * M_PI * (k) / L_2; \
Real q = std::sqrt(q_1 * q_1 + q_2 * q_2); \
Complex influence_1[] = INFLUENCE(q_1, q_2, q, E, nu); \
legacy::MatrixProxy inf_mat(influence_1, 3, 3); \
it.setPointer(&influence(i, j, 0)); \
it = inf_mat; \
} while (0)
// Fill influence points for "symetric" areas
#define FILL_SYM_INFLUENCE(k, l) \
do { \
Real q_1 = 2 * M_PI * (l) / L_1; \
Real q_2 = 2 * M_PI * (k) / L_2; \
Real q = std::sqrt(q_1 * q_1 + q_2 * q_2); \
Complex influence_1[] = INFLUENCE(q_1, q_2, q, E, nu); \
Complex influence_2[] = INFLUENCE(q_2, q_1, q, E, nu); \
legacy::MatrixProxy inf_mat(influence_1, 3, 3); \
it.setPointer(&influence(i, j, 0)); \
it = inf_mat; \
inf_mat.setPointer(influence_2); \
it.setPointer(&influence(j, i, 0)); \
it = inf_mat; \
} while (0)
#pragma omp parallel firstprivate( \
it) // all the following loops are independent
{
// Loops over parts of surface
/*
[ ][ ][ ][ ]
[ ][x][ ][ ]
[ ][ ][ ][ ]
[ ][ ][ ][o]
*/
#pragma omp for schedule(dynamic) // can't collapse triangluar loops
for (Int i = 1; i < size / 2; i++) {
for (Int j = i; j < size / 2; j++) {
FILL_SYM_INFLUENCE(i, j);
}
}
/*
[ ][ ][ ][ ]
[ ][ ][ ][o]
[ ][ ][ ][ ]
[ ][x][ ][ ]
*/
#pragma omp for collapse(2)
for (Int i = size / 2 + 1; i < size; i++) {
for (Int j = 1; j < size / 2; j++) {
FILL_INFLUENCE(i - size, j);
}
}
/*
[ ][x][x][o]
[x][ ][x][ ]
[x][x][x][ ]
[ ][ ][ ][ ]
*/
#pragma omp for
for (Int i = 1; i <= size / 2; i++) {
{
Int j = size / 2;
FILL_SYM_INFLUENCE(i, j);
}
{
Int j = 0;
FILL_SYM_INFLUENCE(i, j);
}
}
/*
[ ][ ][ ][ ]
[ ][ ][ ][ ]
[ ][ ][ ][o]
[x][ ][x][ ]
*/
#pragma omp for
for (Int i = size / 2 + 1; i < size; i++) {
{
Int j = size / 2;
FILL_INFLUENCE(i - size, j);
}
{
Int j = 0;
FILL_INFLUENCE(i - size, j);
}
}
} // End #pragma omp parallel
/* State of surface now:
[0][x][x][o]
[x][x][x][o]
[x][x][x][o]
[x][x][x][o]
Legend: x => explicitely filled
o => hermitian symetry
0 => actually 0
*/
computeInfluenceLipschitz();
#undef INFLUENCE
#undef FILL_INFLUENCE
#undef FILL_SYM_INFLUENCE
}
/* -------------------------------------------------------------------------- */
/// Computing L_2 norm of influence functions
void BemGrid::computeInfluenceLipschitz() {
const UInt N = influence.getNbPoints();
legacy::VectorProxy m(nullptr, 9);
Real _norm = 0.;
#pragma omp parallel for firstprivate(m) reduction(+ : _norm)
for (UInt i = 0; i < N; i++) {
m.setPointer(influence.getInternalData() + i * m.dataSize());
_norm += pow(norm(m.l2norm()), 2);
}
lipschitz = sqrt(_norm);
}
/* -------------------------------------------------------------------------- */
/* Friction related computations */
/* -------------------------------------------------------------------------- */
void BemGrid::projectOnFrictionCone(Grid& field, Real mu) {
legacy::VectorProxy p_proxy(nullptr, 3);
const UInt N = field.getNbPoints();
#pragma omp parallel for firstprivate(p_proxy)
for (UInt i = 0; i < N; i++) {
p_proxy.setPointer(field.getInternalData() + i * p_proxy.dataSize());
projectOnFrictionCone(p_proxy, mu);
}
}
/* -------------------------------------------------------------------------- */
void BemGrid::enforceMeanValue(Grid& field,
const Grid& target) {
Real red_0 = 0, red_1 = 0, red_2 = 0;
auto n = field.sizes();
#pragma omp parallel for collapse(2) reduction(+ : red_0, red_1, red_2)
for (UInt i = 0; i < n[0]; i++) {
for (UInt j = 0; j < n[1]; j++) {
red_0 += field(i, j, 0);
red_1 += field(i, j, 1);
red_2 += field(i, j, 2);
}
}
red_0 /= n[0] * n[1];
red_1 /= n[0] * n[1];
red_2 /= n[0] * n[1];
Real p1_data[3] = {red_0, red_1, red_2};
// Loop for pressure shifting
legacy::VectorProxy p1(p1_data, 3);
legacy::VectorProxy p0((Real*)target.getInternalData(), 3);
legacy::VectorProxy p(nullptr, 3);
const UInt N = field.getNbPoints();
#pragma omp parallel for firstprivate(p0, p1)
for (UInt i = 0; i < N; i++) {
p.setPointer(field.getInternalData() + i * p.dataSize());
p -= p1;
p += p0;
}
}
/* -------------------------------------------------------------------------- */
void BemGrid::computeMeanValueError(Grid& field,
const Grid& target,
Grid& error_vec) {
Real p0 = 0, p1 = 0, p2 = 0;
auto n = field.sizes();
#pragma omp parallel for collapse(2) reduction(+ : p0, p1, p2)
for (UInt i = 0; i < n[0]; i++) {
for (UInt j = 0; j < n[1]; j++) {
p0 += field(i, j, 0);
p1 += field(i, j, 1);
p2 += field(i, j, 2);
}
}
p0 /= n[0] * n[1];
p1 /= n[0] * n[1];
p2 /= n[0] * n[1];
Real error[3] = {p0 - target(0), p1 - target(1), p2 - target(2)};
legacy::VectorProxy e(error, 3);
error_vec = e;
}
/* -------------------------------------------------------------------------- */
Real BemGrid::computeMeanValueError(Grid& field,
const Grid& target) {
Real e[3] = {0.};
legacy::VectorProxy error(e, 3);
computeMeanValueError(field, target, error);
Real norm =
target(0) * target(0) + target(1) * target(1) + target(2) * target(2);
return error.l2norm() / std::sqrt(norm);
}
/* -------------------------------------------------------------------------- */
} // namespace tamaas
diff --git a/src/bem/bem_grid.hh b/src/bem/bem_grid.hh
index 844cd72..b8896e5 100644
--- a/src/bem/bem_grid.hh
+++ b/src/bem/bem_grid.hh
@@ -1,184 +1,187 @@
/**
*
* @author Lucas Frérot
*
* @section LICENSE
*
* Copyright (©) 2016 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 .
*
*/
/* -------------------------------------------------------------------------- */
-#ifndef __BEM_GRID_HH__
-#define __BEM_GRID_HH__
+#ifndef BEM_GRID_HH
+#define BEM_GRID_HH
/* -------------------------------------------------------------------------- */
+#include "fft_engine.hh"
#include "grid.hh"
#include "surface.hh"
#include "legacy_types.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/**
* @brief abstract class for BEMs that use multi component fields
*/
class BemGrid {
public:
/// Constructor
BemGrid(Surface& surface);
/// Destructor
virtual ~BemGrid();
public:
/* ------------------------------------------------------------------------ */
/* Actual computations */
/* ------------------------------------------------------------------------ */
/// Compute normal vectors to surface
void computeSurfaceNormals();
/// Compute influence coefficients
void computeInfluence();
/// Compute a Lipschitz constant of gradient
void computeInfluenceLipschitz();
/// Apply influence fonctions (p -> u)
void applyInfluenceFunctions(Grid& input, Grid& output);
/// Apply inverse influence fonctions (u -> p)
void applyInverseInfluenceFunctions(Grid& input,
Grid& output);
/// Compute true tractions
virtual void computeTrueTractions() = 0;
/// Compute true displacements
virtual void computeTrueDisplacements() = 0;
/// Compute equilibrium
virtual void computeEquilibrium(Real tol, const Grid& p_target) = 0;
/// Compute displacements
void computeDisplacementsFromTractions();
/// Compute tractions
void computeTractionsFromDisplacements();
/// Compute gaps
void computeGaps(Grid& disp);
/* ------------------------------------------------------------------------ */
/* Field related computations */
/* ------------------------------------------------------------------------ */
/// Project point friction cone
static inline void projectOnFrictionCone(legacy::VectorProxy& p, Real mu);
/// Project field friction cone
static void projectOnFrictionCone(Grid& field, Real mu);
/// Project on mean value space
static void enforceMeanValue(Grid& field,
const Grid& target);
/// Compute error relative to target mean value
static Real computeMeanValueError(Grid& field,
const Grid& target);
/// Compute error vector
static void computeMeanValueError(Grid& field,
const Grid& target,
Grid& error_vec);
/* ------------------------------------------------------------------------ */
/* Getters */
/* ------------------------------------------------------------------------ */
/// Get influence
const GridHermitian& getInfluence() const { return influence; }
/// Get normals
const Grid& getSurfaceNormals() const { return surface_normals; }
/// Get tractions
const Grid& getTractions() const { return tractions; }
/// Get displacements
const Grid& getDisplacements() const { return displacements; }
/// Get true tractions
const Grid& getTrueTractions() const;
/// Get true displacements
const Grid& getTrueDisplacements() const;
/// Get gaps
const Grid& getGaps() const { return gaps; }
/* ------------------------------------------------------------------------ */
/* Setters */
/* ------------------------------------------------------------------------ */
/// Elasticity parameters
void setElasticity(Real E_, Real nu_) {
E = E_;
nu = nu_;
}
/// Friction coefficient
void setMu(Real mu_) { mu = mu_; }
/// Max iterations
void setMaxIterations(UInt n) { max_iter = n; }
/// Dump frequency
void setDumpFrequency(UInt n) { dump_freq = n; }
protected:
/// Elasticity constants
Real E, nu;
/// Friction coefficient
Real mu;
/// Lipschitz constant of functional gradient
Real lipschitz;
/// Max iterations
UInt max_iter;
/// Surface
Surface surface;
/// Normal vectors
Grid surface_normals;
/// Influcence matrices
GridHermitian influence;
/// Tractions
Grid tractions;
/// Displacements
Grid displacements;
/// Gap
Grid gaps;
/// True tractions
Grid* true_tractions;
/// True displacements
Grid* true_displacements;
/// Spectral tractions
GridHermitian spectral_tractions;
/// Spectral displacements
GridHermitian spectral_displacements;
/// Dump frequency
UInt dump_freq;
+ /// FFT Engine
+ FFTEngine engine;
};
/* -------------------------------------------------------------------------- */
/* Inline definitions */
/* -------------------------------------------------------------------------- */
inline void BemGrid::projectOnFrictionCone(legacy::VectorProxy& p, Real mu) {
Real pt = std::sqrt(p(0) * p(0) + p(1) * p(1));
Real pn = p(2);
Real x = (pn + mu * pt) / (1 + mu * mu);
if (pt <= mu * pn && pn >= 0) // inside friction cone
return;
else if (mu * pt <= std::abs(pn) &&
pn < 0) // inside "project to origin" cone (undefined normal)
// else if (pn < 0) // project to origin all points with pn < 0
p = 0;
else { // projectable on cone
p(0) *= mu / pt * x;
p(1) *= mu / pt * x;
p(2) = x;
}
}
/* -------------------------------------------------------------------------- */
} // namespace tamaas
#endif // __BEM_GRID_HH__
diff --git a/src/bem/bem_interface.hh b/src/bem/bem_interface.hh
index 86f3877..2cc620b 100644
--- a/src/bem/bem_interface.hh
+++ b/src/bem/bem_interface.hh
@@ -1,100 +1,100 @@
/**
*
* @author Guillaume Anciaux
*
* @section LICENSE
*
* Copyright (©) 2016 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 .
*
*/
/* -------------------------------------------------------------------------- */
#ifndef BEM_INTERFACE_H
#define BEM_INTERFACE_H
-#include "fft_plan_manager.hh"
+#include "fft_engine.hh"
#include "fftransform.hh"
#include "surface.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
class BemInterface {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
BemInterface(Surface& surface)
: surface(), spectral_surface(surface.size(), surface.getL()) {
this->surface.wrap(surface);
this->computeSpectralSurface();
};
virtual ~BemInterface(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
protected:
//! compute the displacements from the tractions
virtual void computeDisplacementsFromTractions(){
SURFACE_FATAL("TO BE IMPLEMENTED")};
//! compute the tractions from the displacements
virtual void computeTractionsFromDisplacements(){
SURFACE_FATAL("TO BE IMPLEMENTED")};
public:
//! compute the spectral surface
virtual void computeSpectralSurface() {
- FFTransform& plan =
- FFTPlanManager::get().createPlan(surface, spectral_surface);
- plan.forwardTransform();
+ engine.forward(surface, spectral_surface);
};
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
//! return the actual tractions
virtual const Surface& getTractions() const = 0;
//! return the actual displacements
virtual const Surface& getDisplacements() const = 0;
//! return the manipulated surface
const Surface& getSurface() const { return this->surface; };
//! return the manipulated surface in fourier space
const SurfaceComplex& getSpectralSurface() const {
return this->spectral_surface;
};
//! return the manipulated surface
void setSurface(Surface& surface) { this->surface = surface; };
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
//! the surface profile we work with
Surface surface;
//! the surface profile in fourier space
SurfaceComplex spectral_surface;
+ //! fft engine
+ FFTEngine engine;
};
} // namespace tamaas
#endif /* BEM_INTERFACE_H */
diff --git a/src/core/surface_statistics.cpp b/src/core/surface_statistics.cpp
index 9cc2178..efdcf15 100644
--- a/src/core/surface_statistics.cpp
+++ b/src/core/surface_statistics.cpp
@@ -1,535 +1,530 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 "surface_statistics.hh"
-#include "fft_plan_manager.hh"
-#include "fftransform.hh"
+#include "fft_engine.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
+namespace {
+FFTEngine engine; ///< local to this translation unit
+}
+
/* -------------------------------------------------------------------------- */
std::vector SurfaceStatistics::computeSpectralDistribution(
const Surface& surface, unsigned int n_bins, UInt cutoff) {
if (n_bins == 0)
SURFACE_FATAL("badly set number of bins");
const int n = surface.size();
Real min_value = SurfaceStatistics::computeMinimum(surface);
Real max_value = SurfaceStatistics::computeMaximum(surface);
std::vector distrib_fft(n_bins);
const int N = n_bins;
Real L = max_value - min_value;
#pragma omp parallel for
for (int k = 0; k < N; ++k)
distrib_fft[k] = 0.;
#pragma omp parallel for
for (int k = 1; k < N / 2; ++k) {
for (int i = 0; i < n * n; ++i) {
Real value = surface(i);
Real angle = -2. * M_PI * k * (value - min_value) / L;
Complex phase = std::polar(1., angle);
distrib_fft[k] += phase;
distrib_fft[N - k] += phase;
}
}
for (int i = 0; i < n * n; ++i) {
Real value = surface(i);
Real angle = -2. * M_PI * N / 2 * (value - min_value) / L;
Complex phase = std::polar(1., angle);
distrib_fft[N / 2] += phase;
}
#pragma omp parallel for
for (int k = 1; k < N; ++k)
distrib_fft[k] /= n * n;
distrib_fft[0] = std::polar(1., 0.);
#pragma omp parallel for
for (int k = 1; k < N / 2; ++k) {
if (k >= (int)cutoff) {
distrib_fft[k] = 0.;
distrib_fft[N - k] = 0.;
}
}
- // proceed to inverse fourier transform
+ // proceed to inverse fourier transform
#pragma omp parallel for
for (int k = 0; k < N; ++k) {
distrib_fft[k] = conj(distrib_fft[k]);
}
auto* in = (fftw_complex*)&distrib_fft[0];
fftw_plan p = fftw_plan_dft_1d(N, in, in, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(p);
fftw_destroy_plan(p);
std::vector distrib;
distrib.resize(N);
#pragma omp parallel for
for (int k = 0; k < N; ++k) {
distrib[k] = conj(distrib_fft[k]);
}
return distrib;
}
/* -------------------------------------------------------------------------- */
Real SurfaceStatistics::computeAverageAmplitude(const Surface& surface) {
return surface.mean();
}
/* -------------------------------------------------------------------------- */
Complex SurfaceStatistics::computeAverageAmplitude(
const SurfaceComplex& surface) {
return surface.mean();
}
/* -------------------------------------------------------------------------- */
Real SurfaceStatistics::computeContactArea(const Surface& surface) {
Real L = surface.getL();
return SurfaceStatistics::computeContactAreaRatio(surface) * L * L;
}
/* -------------------------------------------------------------------------- */
Real SurfaceStatistics::computeContactAreaRatio(const Surface& surface) {
UInt area = 0;
UInt n = surface.size();
#pragma omp parallel for collapse(2), reduction(+ : area)
for (UInt i = 0; i < n; i++) {
for (UInt j = 0; j < n; j++) {
if (surface(i, j) > 1e-10)
area++;
}
}
return area / Real(n * n);
}
/* -------------------------------------------------------------------------- */
Surface
SurfaceStatistics::computeAutocorrelation(Surface& surface) {
SurfaceComplex correlation_heights = computePowerSpectrum(surface);
Surface acf(correlation_heights.size(), correlation_heights.getL());
- FFTransform& plan =
- FFTPlanManager::get().createPlan(acf, correlation_heights);
- plan.backwardTransform();
- FFTPlanManager::get().destroyPlan(acf, correlation_heights);
+ engine.backward(acf, correlation_heights);
return acf;
}
/* -------------------------------------------------------------------------- */
SurfaceComplex
SurfaceStatistics::computePowerSpectrum(Surface& surface) {
SurfaceComplex power_spectrum(surface.size(), surface.getL());
- FFTransform& plan =
- FFTPlanManager::get().createPlan(surface, power_spectrum);
- plan.forwardTransform();
-
+ engine.forward(surface, power_spectrum);
UInt n = power_spectrum.dataSize();
UInt factor = power_spectrum.size() * power_spectrum.size();
#pragma omp parallel for
for (UInt i = 0; i < n; ++i)
power_spectrum(i) /= factor;
power_spectrum.makeItRealBySquare();
- FFTPlanManager::get().destroyPlan(surface, power_spectrum);
return power_spectrum;
}
/* -------------------------------------------------------------------------- */
Real SurfaceStatistics::computeSkewness(const Surface& surface, Real avg,
Real std) {
Real skw = 0.0;
UInt n = surface.size();
#pragma omp parallel for reduction(+ : skw)
for (UInt i = 0; i < n * n; i++) {
Real tmp = surface(i) - avg;
skw += (tmp * tmp * tmp);
}
skw /= (n * std * std * std);
return skw;
}
-/* -------------------------------------------------------------------------- */
+/* -------------------------------------------------------------------------- */
Real SurfaceStatistics::computeKurtosis(const Surface& surface, Real avg,
Real std) {
-
Real kur = 0.0;
UInt n = surface.size();
#pragma omp parallel for reduction(+ : kur)
for (UInt i = 0; i < n * n; i++) {
Real tmp = surface.size() - avg;
kur += (tmp * tmp * tmp * tmp);
}
kur /= (n * std * std * std * std);
return kur - 3;
}
-/* -------------------------------------------------------------------------- */
+/* -------------------------------------------------------------------------- */
Real SurfaceStatistics::computeMaximum(const Surface& surface) {
Real _max = -1. * std::numeric_limits::max();
UInt n = surface.size();
#pragma omp parallel for reduction(max : _max)
for (UInt i = 0; i < n * n; i++) {
_max = std::max(surface(i), _max);
}
return _max;
}
/* -------------------------------------------------------------------------- */
Real SurfaceStatistics::computeMinimum(const Surface& surface) {
Real _min = std::numeric_limits::max();
UInt n = surface.size();
#pragma omp parallel for reduction(min : _min)
for (UInt i = 0; i < n * n; i++) {
_min = std::min(surface(i), _min);
}
return _min;
}
/* -------------------------------------------------------------------------- */
Real SurfaceStatistics::computeRMSSlope(const Surface& surface) {
UInt n = surface.size();
Surface grad(n, 1.);
for (UInt i = 0, m = 0; i < n; i++) {
for (UInt j = 0; j < n; j++) {
UInt k = i + n * j;
// xl = ((i==0)?(k+n-1):(k-1));
UInt xh = ((i == (n - 1)) ? (k - n + 1) : (k + 1));
// yl = ((j==0)?(k+n*n-n):(k-n));
UInt yh = ((j == (n - 1)) ? (k - n * n + n) : (k + n));
Real gx = surface(xh) - surface(k);
Real gy = surface(yh) - surface(k);
grad(m++) = gx * gx + gy * gy;
}
}
Real avg = SurfaceStatistics::computeAverage(grad);
Real slp = sqrt(avg); /* According to paper Hyun PRE 70:026117 (2004) */
Real scale_factor = surface.getL() / n;
return slp / scale_factor;
}
/* -------------------------------------------------------------------------- */
Real SurfaceStatistics::computeSpectralRMSSlope(Surface& surface) {
SurfaceComplex power = computePowerSpectrum(surface);
UInt n = surface.size();
Real rms_slopes = 0;
// /!\ commented code was causing race conditions with dummy variable in
// GridHermitian
#pragma omp parallel for collapse(2), reduction(+ : rms_slopes)
for (UInt i = 1; i < n / 2; ++i)
for (UInt j = 1; j < n / 2; ++j) {
rms_slopes += 2. * ((i * i) + (j * j)) * power(i, j).real();
// rms_slopes += 1.*((i*i)+(j*j))*power(n-i,n-j).real();
rms_slopes += 2. * ((i * i) + (j * j)) * power(n - i, j).real();
// rms_slopes += 1.*((i*i)+(j*j))*power(i,n-j).real();
}
- // external line (i or j == 0)
+ // external line (i or j == 0)
#pragma omp parallel for reduction(+ : rms_slopes)
for (UInt i = 1; i < n / 2; ++i) {
rms_slopes += 1. * (i * i) * power(i, 0).real();
rms_slopes += 1. * (i * i) * power(n - i, 0).real();
rms_slopes += 2. * (i * i) * power(0, i).real();
// rms_slopes += 1.*(i*i)*power(0,n-i).real();
}
- // internal line (i or j == n/2)
+ // internal line (i or j == n/2)
#pragma omp parallel for reduction(+ : rms_slopes)
for (UInt i = 1; i < n / 2; ++i) {
rms_slopes += 0.5 * (i * i + n * n / 4) * power(i, n / 2).real();
rms_slopes += 0.5 * (i * i + n * n / 4) * power(n - i, n / 2).real();
rms_slopes += 1. * (i * i + n * n / 4) * power(n / 2, i).real();
// rms_slopes += 0.5*(i*i+n*n/4)*power(n/2,n-i).real();
}
{
// i == n/2 j == n/2
rms_slopes += 0.25 * (n * n / 2) * power(n / 2, n / 2).real();
// i == 0 j == n/2
rms_slopes += 1. * (n * n / 4) * power(0, n / 2).real();
// i == n/2 j == 0
rms_slopes += 1. * (n * n / 4) * power(n / 2, 0).real();
}
rms_slopes = 2. * M_PI / n * sqrt(rms_slopes);
Real scale_factor = surface.getL() / n;
return rms_slopes / scale_factor;
}
/* -------------------------------------------------------------------------- */
Real SurfaceStatistics::computeStdev(const Grid& surface,
Real epsilon) {
UInt n = surface.dataSize();
UInt cpt = 0;
Real std = 0.0;
Real avg = 0.0;
#pragma omp parallel for reduction(+ : avg, std, cpt)
for (UInt i = 0; i < n; i++) {
Real val = surface(i);
if (val > epsilon) {
avg += val;
std += val * val;
++cpt;
}
}
avg /= cpt;
return sqrt(std / cpt - avg * avg);
}
/* -------------------------------------------------------------------------- */
Real SurfaceStatistics::computeSpectralStdev(Surface& surface) {
SurfaceComplex power = computePowerSpectrum(surface);
Real rms_heights = 0;
UInt n = surface.size();
#pragma omp parallel for collapse(2), reduction(+ : rms_heights)
for (UInt i = 0; i < n; ++i)
for (UInt j = 0; j < n; ++j) {
rms_heights += power(i, j).real();
}
rms_heights = sqrt(rms_heights);
return rms_heights;
}
/* -------------------------------------------------------------------------- */
Real SurfaceStatistics::computeSpectralMeanCurvature(Surface& surface) {
SurfaceComplex power = computePowerSpectrum(surface);
UInt n = surface.size();
Real rms_curvatures = 0;
#pragma omp parallel for collapse(2), reduction(+ : rms_curvatures)
for (UInt i = 1; i < n / 2; ++i)
for (UInt j = 1; j < n / 2; ++j) {
Real coeff = ((i * i) + (j * j));
coeff *= coeff;
rms_curvatures += 1. * coeff * power(i, j).real();
rms_curvatures += 1. * coeff * power(n - i, n - j).real();
rms_curvatures += 1. * coeff * power(n - i, j).real();
rms_curvatures += 1. * coeff * power(i, n - j).real();
}
- // external line (i or j == 0)
+ // external line (i or j == 0)
#pragma omp parallel for reduction(+ : rms_curvatures)
for (UInt i = 1; i < n / 2; ++i) {
Real coeff = i * i;
coeff *= coeff;
rms_curvatures += 4. / 3. * coeff * power(i, 0).real();
rms_curvatures += 4. / 3. * coeff * power(n - i, 0).real();
rms_curvatures += 4. / 3. * coeff * power(0, i).real();
rms_curvatures += 4. / 3. * coeff * power(0, n - i).real();
}
- // internal line (i or j == n/2)
+ // internal line (i or j == n/2)
#pragma omp parallel for reduction(+ : rms_curvatures)
for (UInt i = 1; i < n / 2; ++i) {
Real coeff = (i * i + n * n / 4);
coeff *= coeff;
rms_curvatures += 0.5 * coeff * power(i, n / 2).real();
rms_curvatures += 0.5 * coeff * power(n - i, n / 2).real();
rms_curvatures += 0.5 * coeff * power(n / 2, i).real();
rms_curvatures += 0.5 * coeff * power(n / 2, n - i).real();
}
{
// i == n/2 j == n/2
rms_curvatures +=
0.25 * (n * n / 2) * (n * n / 2) * power(n / 2, n / 2).real();
// i == 0 j == n/2
rms_curvatures +=
2. / 3. * (n * n / 4) * (n * n / 4) * power(0, n / 2).real();
// i == n/2 j == 0
rms_curvatures +=
2. / 3. * (n * n / 4) * (n * n / 4) * power(n / 2, 0).real();
}
rms_curvatures *= 9. / 4.;
rms_curvatures = M_PI / n / n * sqrt(rms_curvatures);
return rms_curvatures;
}
/* -------------------------------------------------------------------------- */
Real SurfaceStatistics::computePerimeter(const Surface& surface,
Real eps) {
int perimeter = 0;
UInt n = surface.size();
#pragma omp parallel for collapse(2), reduction(+ : perimeter)
for (UInt i = 0; i < n; ++i) {
for (UInt j = 0; j < n; ++j) {
if (j > 0 && ((surface(i, (j - 1)) < eps && surface(i, j) > eps) ||
(surface(i, (j - 1)) > eps && surface(i, j) < eps)))
perimeter++;
if (i > 0 && ((surface(j, (i - 1)) < eps && surface(j, i) > eps) ||
(surface(j, (i - 1)) > eps && surface(j, i) < eps)))
perimeter++;
}
}
return perimeter / Real(n);
}
/* -------------------------------------------------------------------------- */
Real SurfaceStatistics::computeSum(const Surface& surface) {
Real sum = 0.;
UInt n = surface.size();
#pragma omp parallel for reduction(+ : sum)
for (UInt i = 0; i < n * n; ++i)
sum += surface(i);
return sum;
}
/* -------------------------------------------------------------------------- */
std::vector SurfaceStatistics::computeMoments(Surface& surface) {
SurfaceComplex power = computePowerSpectrum(surface);
Real L = surface.getL();
UInt n = surface.size();
Real k_factor = 2. * M_PI / L;
std::vector> moment;
moment.resize(6);
for (int i = 0; i < 5; i++)
moment[i].resize(5);
for (UInt i = 0; i < n / 2; i++) {
Real i0 = pow(Real(i) * k_factor, 0);
Real i2 = pow(Real(i) * k_factor, 2);
Real i4 = pow(Real(i) * k_factor, 4);
Real im0 = pow(-Real(i) * k_factor, 0);
Real im2 = pow(-Real(i) * k_factor, 2);
Real im4 = pow(-Real(i) * k_factor, 4);
for (UInt j = 0; j < n / 2; j++) {
Real Phi = power(i, j).real();
Real Phi_10;
Real Phi_01;
Real Phi_11;
if (i == 0 && j == 0) {
Phi_10 = 0;
Phi_11 = 0;
Phi_01 = 0;
} else if (i == 0) {
Phi_10 = 0;
Phi_11 = 0;
Phi_01 = power(i, (n - j)).real();
} else if (j == 0) {
Phi_10 = power(n - i, j).real();
Phi_11 = 0;
Phi_01 = 0;
} else {
Phi_10 = power(n - i, j).real();
Phi_11 = power(n - i, (n - j)).real();
Phi_01 = power(i, (n - j)).real();
}
Real j0 = pow(Real(j) * k_factor, 0);
Real j2 = pow(Real(j) * k_factor, 2);
Real j4 = pow(Real(j) * k_factor, 4);
Real jm0 = pow(-Real(j) * k_factor, 0);
Real jm2 = pow(-Real(j) * k_factor, 2);
Real jm4 = pow(-Real(j) * k_factor, 4);
moment[0][0] += i0 * j0 * Phi + im0 * j0 * Phi_10 + i0 * jm0 * Phi_01 +
im0 * jm0 * Phi_11;
moment[0][2] += i0 * j2 * Phi + im0 * j2 * Phi_10 + i0 * jm2 * Phi_01 +
im0 * jm2 * Phi_11;
moment[2][0] += i2 * j0 * Phi + im2 * j0 * Phi_10 + i2 * jm0 * Phi_01 +
im2 * jm0 * Phi_11;
moment[2][2] += i2 * j2 * Phi + im2 * j2 * Phi_10 + i2 * jm2 * Phi_01 +
im2 * jm2 * Phi_11;
moment[0][4] += i0 * j4 * Phi + im0 * j4 * Phi_10 + i0 * jm4 * Phi_01 +
im0 * jm4 * Phi_11;
moment[4][0] += i4 * j0 * Phi + im4 * j0 * Phi_10 + i4 * jm0 * Phi_01 +
im4 * jm0 * Phi_11;
}
}
std::vector mean_moments;
// m2
Real mean_m2 = 0.5 * (moment[0][2] + moment[2][0]);
mean_moments.push_back(mean_m2);
// m4
Real mean_m4 = (3 * moment[2][2] + moment[0][4] + moment[4][0]) / 3.;
mean_moments.push_back(mean_m4);
return mean_moments;
}
/* -------------------------------------------------------------------------- */
Real SurfaceStatistics::computeAnalyticFractalMoment(UInt order, UInt k1,
UInt k2, Real hurst,
Real A, Real L) {
Real xi = k2 / k1;
Real _k1 = 2. * M_PI / L * k1;
Real Cq;
switch (order) {
case 0:
Cq = 2. * M_PI;
break;
case 2:
Cq = M_PI;
break;
case 4:
Cq = 3. / 4. * M_PI;
break;
default:
SURFACE_FATAL("invalid order");
}
Real exp = order - 2. * hurst;
Real res = A * Cq * pow(_k1, exp) / exp * (pow(xi, exp) - 1);
return res;
}
/* -------------------------------------------------------------------------- */
-}
+} // namespace tamaas
diff --git a/src/core/tamaas.cpp b/src/core/tamaas.cpp
index 04200e1..8ff2c25 100644
--- a/src/core/tamaas.cpp
+++ b/src/core/tamaas.cpp
@@ -1,56 +1,55 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 "tamaas.hh"
#include "fft_plan_manager.hh"
#include "loop.hh"
#include
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
void initialize(UInt num_threads) {
if (num_threads)
omp_set_num_threads(num_threads); // set user-defined number of threads
// fftw_import_wisdom_from_filename("fftransform_fftw_wisdom");
#ifndef USE_CUDA
if (!fftw_init_threads())
TAMAAS_EXCEPTION("FFTW could not initialize threads!");
fftw_plan_with_nthreads(omp_get_max_threads());
#endif
}
/* -------------------------------------------------------------------------- */
void finalize() {
- FFTPlanManager::get().clean();
// fftw_export_wisdom_to_filename("fftransform_fftw_wisdom");
fftw_cleanup_threads();
}
} // namespace tamaas
diff --git a/src/model/volume_potential.cpp b/src/model/volume_potential.cpp
index d30fb84..5cad2d7 100644
--- a/src/model/volume_potential.cpp
+++ b/src/model/volume_potential.cpp
@@ -1,115 +1,113 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 "volume_potential.hh"
#include "model.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
template
VolumePotential::VolumePotential(Model* model) : IntegralOperator(model) {
// Copy horizontal sizes
std::array sizes;
std::copy(model->getDiscretization().begin() + 1,
model->getDiscretization().end(), sizes.begin());
auto hermitian_sizes =
GridHermitian::hermitianDimensions(
sizes);
wavevectors =
- FFTransform::template computeFrequencies<
- true>(hermitian_sizes);
+ FFTEngine::template computeFrequencies(hermitian_sizes);
// Normalize wavevectors
auto system_size = model->getBoundarySystemSize();
VectorProxy boundary_domain{
system_size[0]};
wavevectors *= 2 * M_PI;
wavevectors /= boundary_domain;
wavevectors *= -1.; // < this is important for the convolution computation
}
/* -------------------------------------------------------------------------- */
template
model_type VolumePotential::getType() const {
return this->model->getType();
}
/* -------------------------------------------------------------------------- */
template
void VolumePotential::initialize(UInt source_components,
UInt out_components,
UInt out_buffer_size) {
auto hermitian_sizes =
GridHermitian::hermitianDimensions(
this->model->getBoundaryDiscretization());
// Initializing buffers
this->source_buffer.resize(this->model->getDiscretization()[0]);
this->out_buffer.resize(out_buffer_size);
// Resizing source buffer
std::for_each(this->source_buffer.begin(), this->source_buffer.end(),
[&](auto& buffer) {
buffer.setNbComponents(source_components);
buffer.resize(hermitian_sizes);
});
// Resizing output buffer
std::for_each(this->out_buffer.begin(), this->out_buffer.end(),
[&](auto& buffer) {
buffer.setNbComponents(out_components);
buffer.resize(hermitian_sizes);
});
}
/* -------------------------------------------------------------------------- */
template
void VolumePotential::transformSource(GridBase& in,
filter_t pred) const {
constexpr UInt dim = trait::dimension;
auto& i = dynamic_cast&>(in);
// Transforming source
for (UInt layer : Loop::range(i.sizes().front())) {
if (not pred(layer))
continue;
auto in_layer = make_view(i, layer);
- FFTPlanManager::get()
- .createPlan(in_layer, source_buffer[layer])
- .forwardTransform();
+ engine.forward(in_layer, source_buffer[layer]);
}
- for (auto && layer : out_buffer)
+ for (auto&& layer : out_buffer)
layer = 0;
}
/* -------------------------------------------------------------------------- */
/* Template instanciation */
/* -------------------------------------------------------------------------- */
template class VolumePotential;
template class VolumePotential;
} // namespace tamaas
diff --git a/src/model/volume_potential.hh b/src/model/volume_potential.hh
index f0d15f9..d2da1fd 100644
--- a/src/model/volume_potential.hh
+++ b/src/model/volume_potential.hh
@@ -1,142 +1,141 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 VOLUME_POTENTIAL_HH
#define VOLUME_POTENTIAL_HH
/* -------------------------------------------------------------------------- */
-#include "fft_plan_manager.hh"
+#include "fft_engine.hh"
#include "grid_hermitian.hh"
#include "grid_view.hh"
#include "integral_operator.hh"
#include "logger.hh"
#include "model_type.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/// Trait type for component management
template
struct derivative_traits;
template <>
struct derivative_traits<0> {
template
static constexpr auto source_components = model_type_traits::components;
template
static constexpr auto out_components = model_type_traits::components;
};
template <>
struct derivative_traits<1> {
template
static constexpr auto source_components = model_type_traits::voigt;
template
static constexpr auto out_components = model_type_traits::components;
};
template <>
struct derivative_traits<2> {
template
static constexpr auto source_components = model_type_traits::voigt;
template
static constexpr auto out_components = model_type_traits::voigt;
};
/**
* @brief Volume potential operator class. Applies the operators for computation
* of displacements and strains due to residual/eigen strains
*/
template
class VolumePotential : public IntegralOperator {
using trait = model_type_traits;
protected:
using filter_t = std::function;
public:
VolumePotential(Model* model);
/// Update from model (does nothing)
void updateFromModel() override {}
/// Kind
IntegralOperator::kind getKind() const override {
return IntegralOperator::neumann;
}
/// Type
model_type getType() const override;
/// Apply to all of the source layers
void apply(GridBase& input, GridBase& output) const override {
applyIf(input, output, [](UInt) { return true; });
}
protected:
/// Transform source layer-by-layer
void transformSource(GridBase& in, filter_t pred) const;
/// Transform all source
void transformSource(GridBase& in) const {
transformSource(in, [](auto) { return true; });
}
/// Transform output layer-by-layer
template
void transformOutput(Func func, GridBase& out) const;
/// Initialize fourier buffers
void initialize(UInt source_components, UInt out_components,
UInt source_buffer_size);
protected:
Grid wavevectors;
using BufferType = GridHermitian;
mutable std::vector source_buffer;
mutable std::vector out_buffer;
+ mutable FFTEngine engine;
};
/* -------------------------------------------------------------------------- */
/* Template implementation */
/* -------------------------------------------------------------------------- */
template
template
void VolumePotential::transformOutput(Func func,
GridBase& out) const {
constexpr UInt dim = trait::dimension;
auto& o = dynamic_cast&>(out);
// Transforming output
for (UInt layer : Loop::range(o.sizes().front())) {
Logger().get(LogLevel::debug)
<< "[VolumePotential] Computing layer " << layer << '\n';
auto out_layer = make_view(o, layer);
auto& fourier_out_layer = func(out_buffer, layer);
- FFTPlanManager::get()
- .createPlan(out_layer, fourier_out_layer)
- .backwardTransform();
+ engine.backward(out_layer, fourier_out_layer);
}
}
} // namespace tamaas
#endif // VOLUME_POTENTIAL_HH
diff --git a/src/solvers/contact_solver.cpp b/src/solvers/contact_solver.cpp
index 79c2317..4812cf8 100644
--- a/src/solvers/contact_solver.cpp
+++ b/src/solvers/contact_solver.cpp
@@ -1,71 +1,65 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 "contact_solver.hh"
#include "fft_plan_manager.hh"
#include "logger.hh"
#include
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
ContactSolver::ContactSolver(Model& model, const GridBase& surface,
Real tolerance)
: model(model), surface(), functional(), tolerance(tolerance),
surface_stddev(std::sqrt(surface.var())) {
auto discretization = model.getBoundaryDiscretization();
auto n_points = model.getTraction().getNbPoints();
if (n_points != surface.dataSize())
TAMAAS_EXCEPTION("Model size and surface size do not match!");
this->surface.wrap(surface);
_gap = allocateGrid(model.getType(), discretization);
model.registerField("gap", _gap);
}
/* -------------------------------------------------------------------------- */
-ContactSolver::~ContactSolver() {
- FFTPlanManager::get().clean(); /// TODO find a smarter way to clean up
-}
-
-/* -------------------------------------------------------------------------- */
-
void ContactSolver::printState(UInt iter, Real cost_f, Real error) {
// UInt precision = std::cout.precision();
if (iter % dump_frequency)
return;
Logger().get(LogLevel::info)
<< std::setw(5) << iter << " " << std::setw(15) << std::scientific
<< cost_f << " " << std::setw(15) << error << std::endl
<< std::fixed;
}
/* -------------------------------------------------------------------------- */
void ContactSolver::addFunctionalTerm(functional::Functional* func) {
functional.addFunctionalTerm(func, false);
}
} // namespace tamaas
/* -------------------------------------------------------------------------- */
diff --git a/src/solvers/contact_solver.hh b/src/solvers/contact_solver.hh
index 45c2fa5..51b28f8 100644
--- a/src/solvers/contact_solver.hh
+++ b/src/solvers/contact_solver.hh
@@ -1,77 +1,77 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 __CONTACT_SOLVER_HH__
#define __CONTACT_SOLVER_HH__
/* -------------------------------------------------------------------------- */
#include "model.hh"
#include "tamaas.hh"
#include "meta_functional.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
class ContactSolver {
public:
/// Constructor
ContactSolver(Model& model, const GridBase& surface, Real tolerance);
/// Destructor
- virtual ~ContactSolver();
+ virtual ~ContactSolver() = default;
public:
/// Print state of solve
virtual void printState(UInt iter, Real cost_f, Real error);
/// Set maximum number of iterations
void setMaxIterations(UInt n) { max_iterations = n; }
/// Set dump_frequency
void setDumpFrequency(UInt n) { dump_frequency = n; }
/// Add term to functional
void addFunctionalTerm(functional::Functional* func);
// Virtual functions
public:
virtual Real solve(std::vector /*load*/) {
TAMAAS_EXCEPTION("solve() not implemented");
}
virtual Real solve(Real load) {
return solve(std::vector{load});
}
public:
GridBase& getSurface() { return surface; }
Model& getModel() { return model; }
protected:
Model& model;
GridBase surface;
std::shared_ptr> _gap = nullptr;
functional::MetaFunctional functional;
Real tolerance;
UInt max_iterations = 1000;
Real surface_stddev;
UInt dump_frequency = 100;
};
} // namespace tamaas
/* -------------------------------------------------------------------------- */
#endif // __CONTACT_SOLVER_HH__
diff --git a/src/surface/filter.hh b/src/surface/filter.hh
index e4f9af1..aedff77 100644
--- a/src/surface/filter.hh
+++ b/src/surface/filter.hh
@@ -1,65 +1,65 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 FILTER_H
#define FILTER_H
/* -------------------------------------------------------------------------- */
-#include "fftransform.hh"
+#include "fft_engine.hh"
#include "grid.hh"
#include "grid_hermitian.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
template
class Filter {
public:
/// Default constructor
Filter() = default;
/// Destructor
virtual ~Filter() = default;
public:
/// Compute Filter coefficients
virtual void
computeFilter(GridHermitian& filter_coefficients) const = 0;
protected:
/// Compute filter coefficients using lambda
template
void computeFilter(T&& f, GridHermitian& filter) const;
};
template
template
void Filter::computeFilter(T&& f, GridHermitian& filter) const {
auto wavevectors =
- FFTransform::template computeFrequencies(filter.sizes());
+ FFTEngine::template computeFrequencies(filter.sizes());
Loop::loop(std::forward(f), filter,
range>(wavevectors));
}
} // namespace tamaas
#endif // FILTER_H
diff --git a/tests/test_integration.cpp b/tests/test_integration.cpp
index 253df49..ff4726e 100644
--- a/tests/test_integration.cpp
+++ b/tests/test_integration.cpp
@@ -1,164 +1,164 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 "fftransform.hh"
+#include "fft_engine.hh"
#include "integration/accumulator.hh"
#include "test.hh"
#include
#include
/* -------------------------------------------------------------------------- */
using namespace tamaas;
namespace ex = expolit;
static Real tol = 1e-13;
class AccumulatorTest : public ::testing::Test {
protected:
void SetUp() override {
nodal_values.resize(n);
for (auto&& grid : nodal_values) {
grid.setNbComponents(6);
grid.resize({nh, nh});
grid = 1;
}
accumulator.makeUniformMesh(nodal_values.size(), size);
- wavevectors = FFTransform::computeFrequencies({nh, nh});
+ wavevectors = FFTEngine::computeFrequencies({nh, nh});
}
UInt n = 10, nh = 10;
Real size = 3;
std::vector> nodal_values;
Grid wavevectors;
Accumulator> accumulator;
};
TEST_F(AccumulatorTest, uniformMesh) {
auto& pos = this->accumulator.nodePositions();
std::vector sol(pos.size());
std::iota(sol.begin(), sol.end(), 0);
std::transform(sol.begin(), sol.end(), sol.begin(),
[&](auto& x) { return size * x / (sol.size() - 1); });
ASSERT_TRUE(compare(sol, pos, AreFloatEqual())) << "Uniform mesh fail";
}
// Litteral type
using Q = ex::Litteral;
TEST_F(AccumulatorTest, forward) {
// Setup for layer checking
std::vector layers_seen, layers_to_see(n);
std::iota(layers_to_see.begin(), layers_to_see.end(), 0);
// Setup for integral checking
constexpr Q q;
constexpr auto z = ex::Polynomial({0, 1});
constexpr auto g0 = ex::exp(q * z);
constexpr auto g1 = q * z * ex::exp(q * z);
for (auto&& tuple :
accumulator.forward(this->nodal_values, this->wavevectors)) {
auto&& l = std::get<0>(tuple);
auto&& xl_ = std::get<1>(tuple);
auto&& acc_g0 = std::get<2>(tuple);
auto&& acc_g1 = std::get<3>(tuple);
layers_seen.push_back(l);
// fundamental mode
auto testing = acc_g0(0, 0, 0);
ASSERT_NEAR(testing.real(), xl_, tol) << "Fundamental mode G0 fail";
testing = acc_g1(0, 0, 0);
ASSERT_NEAR(testing.real(), 0, tol) << "Fundamental mode G1 fail";
// other modes
testing = acc_g0(1, 0, 0);
ASSERT_NEAR(testing.real(), acc_g0(0, 1, 0).real(), tol) << "Symmetry fail";
ASSERT_NEAR(testing.imag(), acc_g0(0, 1, 0).imag(), tol) << "Symmetry fail";
testing = acc_g0(1, 1, 0);
Real sqrt_two = std::sqrt(Real(2));
auto ref_value = ex::definite_integral(
std::pair(0, xl_),
ex::substitute(ex::Constant({sqrt_two}), g0));
ASSERT_NEAR(testing.real(), ref_value, tol)
<< "Integration of exp(qy)*ϕ(y) fail";
testing = acc_g1(1, 1, 0);
ref_value = ex::definite_integral(
std::pair(0, xl_),
ex::substitute(ex::Constant