Homec4science
Diffusion tamaas (archive/legacy_mpi) Browse

tamaas/src/bem/bem_fft_base.cpp
a757e2311dc9archive/legacy_mpi

Last Change
Raw File

bem_fft_base.cpp
View Options

/**
* @file
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "bem_fft_base.hh"
/* -------------------------------------------------------------------------- */
#define TIMER
#include "surface_timer.hh"
#include "surface_statistics.hh"
#include "bem_meta_functional.hh"
#include "elastic_energy_functional.hh"
#include "fft_plan_manager.hh"
/* -------------------------------------------------------------------------- */
#include <omp.h>
/* -------------------------------------------------------------------------- */
__BEGIN_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(){
STARTTIMER("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();
#pragma omp parallel for firstprivate(inf_it)
for (auto it = wavevectors.begin(2) ;
it < wavevectors.end(2) ;
++it) {
inf_it = inf_begin;
inf_it += it - wavevectors.begin(2);
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.;
STOPTIMER("computeSpectralInfluenceOverDisplacement");
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::computeDisplacementsFromTractions(){
STARTTIMER("computeDisplacementsFromTractions");
this->applyInfluenceFunctions(surface_tractions,surface_displacements);
STOPTIMER("computeDisplacementFromTractions");
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::computeTractionsFromDisplacements(){
STARTTIMER("computeDisplacementsFromTractions");
this->applyInverseInfluenceFunctions(surface_displacements,surface_tractions);
STOPTIMER("computeDisplacementFromTractions");
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::applyInfluenceFunctions(Surface<Real> & input, Surface<Real> & output){
STARTTIMER("applyInfluenceFunctions");
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();
STOPTIMER("applyInfluenceFunctions");
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::applyInverseInfluenceFunctions(Surface<Real> & input, Surface<Real> & output){
STARTTIMER("applyInfluenceFunctions");
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();
STOPTIMER("applyInfluenceFunctions");
}
/* -------------------------------------------------------------------------- */
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);
}
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
c4science · Help