Page MenuHomec4science
Files F64418656

bem_fft_base.cpp
No OneTemporary
Actions

File Metadata

Created
Sun, May 26, 18:16

bem_fft_base.cpp
View Options

/**
*
* @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 "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("")),
nthreads(1),
tmp_coeff(1./M_PI){
// this->setNumberOfThreads(this->nthreads);
max_iterations = 100000;
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() {
delete functional;
}
/* -------------------------------------------------------------------------- */
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::setNumberOfThreads(UInt nthreads){
this->nthreads = nthreads;
omp_set_num_threads(nthreads);
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::setEffectiveModulus(Real e_star){
this->e_star = e_star;
this->computeSpectralInfluenceOverDisplacement();
}
/* -------------------------------------------------------------------------- */
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");
UInt n = surface.size();
surface_spectral_influence_disp(0) = 0.;
// everything but the external lines and central lines
#pragma omp parallel for
for (UInt i = 1; i < n/2; ++i) {
for (UInt j = 1; j < n/2; ++j) {
Real d = sqrt(i*i+j*j);
surface_spectral_influence_disp(i,j) = 1./d;
surface_spectral_influence_disp(n-i,j) = 1./d;
surface_spectral_influence_disp(i,n-j) = 1./d;
surface_spectral_influence_disp(n-i,n-j) = 1./d;
}
}
// external line (i or j == 0)
#pragma omp parallel for
for (UInt i = 1; i < n/2; ++i) {
Real d = i;
surface_spectral_influence_disp(i,0) = 1./d;
surface_spectral_influence_disp(n-i,0) = 1./d;
surface_spectral_influence_disp(0,i) = 1./d;
surface_spectral_influence_disp(0,n-i) = 1./d;
}
//internal line (i or j == n/2)
#pragma omp parallel for
for (UInt i = 1; i < n/2; ++i) {
Real d = sqrt(i*i+n/2*n/2);
surface_spectral_influence_disp(i,n/2) = 1./d;
surface_spectral_influence_disp(n-i,n/2) = 1./d;
surface_spectral_influence_disp(n/2,i) = 1./d;
surface_spectral_influence_disp(n/2,n-i) = 1./d;
}
{
//i == n/2 j == n/2
Real d = sqrt(2*n/2*n/2);
surface_spectral_influence_disp(n/2,n/2) = 1./d;
//i == 0 j == n/2
d = n/2;
surface_spectral_influence_disp(0,n/2) = 1./d;
//i == n/2 j == 0
surface_spectral_influence_disp(n/2,0) = 1./d;
}
UInt size = n*n;
Real L = surface.getL();
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
surface_spectral_influence_disp(i) *= tmp_coeff/this->e_star*L;
}
STOPTIMER("computeSpectralInfluenceOverDisplacement");
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::computeDisplacementsFromTractions(){
STARTTIMER("computeDisplacementsFromTractions");
this->applyInfluenceFunctions(surface_tractions,surface_displacements);
STOPTIMER("computeDisplacementFromTractions");
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::computeTractionsFromDisplacements(){
STARTTIMER("computeDisplacementsFromTractions");
this->applyInverseInfluenceFunctions(true_displacements,surface_tractions);
STOPTIMER("computeDisplacementFromTractions");
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::applyInfluenceFunctions(Surface<Real> & input, Surface<Real> & output){
STARTTIMER("applyInfluenceFunctions");
FFTransform<Real> & plan1 = FFTPlanManager::get().createPlan(input,surface_spectral_output);
plan1.forwardTransform();
surface_spectral_output *= surface_spectral_influence_disp;
FFTransform<Real> & 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> & plan1 = FFTPlanManager::get().createPlan(input,surface_spectral_output);
plan1.forwardTransform();
UInt n = surface_spectral_input.size();
UInt size = n*n;
surface_spectral_output(0) = 0;
#pragma omp parallel for
for (UInt i = 1; i < size; ++i) {
surface_spectral_output(i) /= surface_spectral_influence_disp(i);
}
FFTransform<Real> & 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);
gap(i) = surface(i) - true_displacements(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__

Event Timeline

c4science · Help