Homec4science
Diffusion tamaas (maxwell_viscoelastic-branch) Browse

tamaas/src/bem/bem_fft_base.cpp
cb09ff877dabmaxwell_viscoelastic…

Last Change
Raw File

bem_fft_base.cpp
View Options

/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-19 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 <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#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 <omp.h>
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
BemFFTBase::BemFFTBase(Surface<Real>& 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<Real>& BemFFTBase::getTractions() const {
return surface_tractions;
}
/* -------------------------------------------------------------------------- */
const Surface<Real>& BemFFTBase::getDisplacements() const {
return surface_displacements;
}
/* -------------------------------------------------------------------------- */
const Surface<Real>& 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<Real, 2>::computeFrequencies<true>(
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<Real> 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<Real>& input,
Surface<Real>& output) {
FFTransform<Real, 2>& plan1 =
FFTPlanManager::get().createPlan(input, surface_spectral_output);
plan1.forwardTransform();
surface_spectral_output *= surface_spectral_influence_disp;
FFTransform<Real, 2>& plan2 =
FFTPlanManager::get().createPlan(output, surface_spectral_output);
plan2.backwardTransform();
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::applyInverseInfluenceFunctions(Surface<Real>& input,
Surface<Real>& output) {
FFTransform<Real, 2>& plan1 =
FFTPlanManager::get().createPlan(input, surface_spectral_output);
plan1.forwardTransform();
surface_spectral_output /= surface_spectral_influence_disp;
surface_spectral_output(0) = 0;
FFTransform<Real, 2>& plan2 =
FFTPlanManager::get().createPlan(output, surface_spectral_output);
plan2.backwardTransform();
}
/* -------------------------------------------------------------------------- */
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
c4science · Help