Page MenuHomec4science

bem_fft_base.cpp
No OneTemporary

File Metadata

Created
Mon, May 27, 04:13

bem_fft_base.cpp

/**
*
* @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;
}
/* -------------------------------------------------------------------------- */
const SurfaceComplex<Real> & BemFFTBase::getSpectralDisplacements()const{
return surface_spectral_displacements;
}
/* -------------------------------------------------------------------------- */
const SurfaceComplex<Real> & BemFFTBase::getSpectralTractions()const{
return surface_spectral_tractions;
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::setNumberOfThreads(UInt nthreads){
this->nthreads = nthreads;
std::cerr << "Trying to set the number of threads to " << nthreads << std::endl;
#pragma omp parallel
{
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;
// std::cerr << "000000000000 " << surface_spectral_tractions(n*n-1) << std::endl;
// std::cerr << "000000000000 " << surface_spectral_displacements(n*n-1) << std::endl;
surface_spectral_displacements = surface_spectral_tractions;
// std::cerr << "AAAAAAAAAAAA " << surface_spectral_tractions(n*n-1) << std::endl;
// std::cerr << "AAAAAAAAAAAA " << surface_spectral_displacements(n*n-1) << std::endl;
// std::cerr << "AAAAAAAAAAAA " << &surface_spectral_displacements(n*n-1) << std::endl;
// std::cerr << "AAAAAAAAAAAA " << &surface_spectral_displacements(0) << std::endl;
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
surface_spectral_displacements(i) *= surface_spectral_influence_disp(i);
}
// std::cerr << "BBBBBBBBBBBB " << surface_spectral_displacements(n*n-1) << std::endl;
// std::cerr << "BBBBBBBBBBBB " << surface_spectral_influence_disp(n*n-1) << std::endl;
STOPTIMER("computeSpectralDisplacement");
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::computeSpectralTractions(){
STARTTIMER("computeSpectralTraction");
surface_tractions.FFTTransform(surface_spectral_tractions,nthreads);
STOPTIMER("computeSpectralTraction");
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::computeDisplacements(){
STARTTIMER("computeDisplacement");
this->computeSpectralTractions();
this->computeSpectralDisplacements();
this->computeRealDisplacementsFromSpectral();
STOPTIMER("computeDisplacement");
}
/* -------------------------------------------------------------------------- */
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();
for (UInt i = 0 ; i < size*size ; i++) {
gap(i) = true_displacements(i) - surface(i);
// if (std::norm(surface_tractions(i)) > zero_threshold*zero_threshold)
// gap(i) = 0.;
}
}
/* -------------------------------------------------------------------------- */
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;
}
/* -------------------------------------------------------------------------- */

Event Timeline