Homec4science
Diffusion tamaas (features/declarative) Browse

tamaas/src/bem_fft_base.cpp
47f0e11cf893features/declarative

Last Change
Raw File

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 <omp.h>
/* -------------------------------------------------------------------------- */
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;
std::cerr << "Trying to set the number of threads to " << nthreads << std::endl;
omp_set_num_threads(nthreads);
std::cerr << "The number of threads have been set to " << omp_get_num_threads() << std::endl;
}
/* -------------------------------------------------------------------------- */
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 (unsigned int i = 1; i < n/2; ++i) {
for (unsigned int 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 (unsigned int 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 (unsigned int 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;
}
unsigned int size = n*n;
Real L = surface.getL();
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
surface_spectral_influence_disp(i) *= tmp_coeff/this->e_star*L;
}
STOPTIMER("computeSpectralInfluenceOverDisplacement");
}
/* -------------------------------------------------------------------------- */
// void BemFFTBase::computeSpectralDisplacements(){
// STARTTIMER("computeSpectralDisplacement");
// UInt n = surface.size();
// UInt size = n*n;
// surface_spectral_displacements = surface_spectral_tractions;
// #pragma omp parallel for
// for (UInt i = 0; i < size; ++i) {
// surface_spectral_displacements(i) *= surface_spectral_influence_disp(i);
// }
// STOPTIMER("computeSpectralDisplacement");
// }
/* -------------------------------------------------------------------------- */
// void BemFFTBase::computeSpectralTractions(){
// STARTTIMER("computeSpectralTractions");
// UInt n = surface.size();
// UInt size = n*n;
// surface_spectral_tractions = surface_spectral_displacements;
// #pragma omp parallel for
// for (UInt i = 0; i < size; ++i) {
// surface_spectral_tractions(i) /= surface_spectral_influence_disp(i);
// }
// STOPTIMER("computeSpectralTractions");
// }
/* -------------------------------------------------------------------------- */
// void BemFFTBase::computeSpectralTractionsFromReal(){
// STARTTIMER("computeSpectralTractionFromReal");
// surface_tractions.FFTTransform(surface_spectral_tractions,nthreads);
// STOPTIMER("computeSpectralTractionFromReal");
// }
/* -------------------------------------------------------------------------- */
// void BemFFTBase::computeSpectralDisplacementsFromReal(){
// STARTTIMER("computeSpectralDisplacementsFromReal");
// surface_displacements.FFTTransform(surface_spectral_displacements,nthreads);
// STOPTIMER("computeSpectralDisplacementsFromReal");
// }
/* -------------------------------------------------------------------------- */
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(const Surface<Real> & input, Surface<Real> & output){
STARTTIMER("applyInfluenceFunctions");
input.FFTTransform(surface_spectral_input,nthreads);
surface_spectral_output = surface_spectral_input;
surface_spectral_output *= surface_spectral_influence_disp;
surface_spectral_output.FFTITransform(output,nthreads);
STOPTIMER("applyInfluenceFunctions");
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::applyInverseInfluenceFunctions(const Surface<Real> & input, Surface<Real> & output){
STARTTIMER("applyInfluenceFunctions");
input.FFTTransform(surface_spectral_input,nthreads);
surface_spectral_output = surface_spectral_input;
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);
}
// //internal line (i or j == n/2)
// #pragma omp parallel for
// for (unsigned int i = 1; i < n/2; ++i) {
// surface_spectral_output(i,n/2) = 0.;
// surface_spectral_output(n-i,n/2) = 0.;
// surface_spectral_output(n/2,i) = 0.;
// surface_spectral_output(n/2,n-i) = 0.;
// }
// {
// //i == n/2 j == n/2
// surface_spectral_output(n/2,n/2) = 0.;
// //i == 0 j == n/2
// surface_spectral_output(0,n/2) = 0.;
// //i == n/2 j == 0
// surface_spectral_output(n/2,0) = 0.;
// }
surface_spectral_output.FFTITransform(output,nthreads);
STOPTIMER("applyInfluenceFunctions");
}
/* -------------------------------------------------------------------------- */
// void BemFFTBase::computeRealTractionsFromSpectral(){
// STARTTIMER("computeRealTractionsFromSpectral");
// surface_spectral_tractions.FFTITransform(surface_tractions,nthreads);
// STOPTIMER("computeRealDisplacementFromSpectral");
// }
/* -------------------------------------------------------------------------- */
// void BemFFTBase::computeRealDisplacementsFromSpectral(){
// STARTTIMER("computeRealDisplacementFromSpectral");
// surface_spectral_displacements.FFTITransform(surface_displacements,nthreads);
// STOPTIMER("computeRealDisplacementFromSpectral");
// }
/* -------------------------------------------------------------------------- */
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;
unsigned int n = surface.size();
unsigned int size = n*n;
Real shift = 1e300;
for (unsigned int i = 0; i < size; ++i) {
if (surface_displacements(i) - shift - surface(i) < 0.){
shift = surface_displacements(i) - surface(i);
}
}
for (unsigned int i = 0; i < size; ++i) {
true_displacements(i) = surface_displacements(i) - shift;
}
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::setFunctional(Functional *new_funct) {
delete this->functional;
this->functional = new_funct;
}
/* -------------------------------------------------------------------------- */
c4science · Help