Homec4science
Diffusion tamaas (master) Browse

tamaas/src/bem_fft_base.cpp
fbfcad1b6cb9master

Last Change
Raw File

bem_fft_base.cpp
View Options

#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;
}
/* -------------------------------------------------------------------------- */
Real BemFFTBase::getF(){
return computed_F;
}
/* -------------------------------------------------------------------------- */
Real BemFFTBase::getOldF(){
return computed_oldF;
}
/* -------------------------------------------------------------------------- */
const Surface<Real> & BemFFTBase::getGradF() const{
return gradient_F;
}
/* -------------------------------------------------------------------------- */
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 (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::computeGradF(){
STARTTIMER("computeGradient");
unsigned int n = surface.size();
unsigned int size = n*n;
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
gradient_F(i) =
(surface_displacements(i) - surface(i));
}
STOPTIMER("computeGradient");
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::computeF(){
STARTTIMER("computeF");
unsigned int n = surface.size();
unsigned int size = n*n;
Real res = 0;
Real d_max = SurfaceStatistics::computeMaximum(surface_displacements);
Real t_sum = std::abs(SurfaceStatistics::computeSum(surface_tractions));
Real s_max = SurfaceStatistics::computeMaximum(surface);
#pragma omp parallel for reduction(+:res)
for (unsigned int i = 0; i < size; ++i) {
res += std::abs(surface_tractions(i)
*(surface_displacements(i)-d_max-surface(i)+s_max)
);
// res +=
// this->surface_tractions[i].real()
// *(surface_displacements[i].real() - surface[i].real());
}
STOPTIMER("computeF");
// this->computed_F = -1.0*res;
this->computed_F = res/std::abs(t_sum);
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::computeOldF(){
STARTTIMER("computeF");
unsigned int n = surface.size();
unsigned int size = n*n;
Real res = 0;
// Real d_max = SurfaceStatistics::computeMaximum(surface_displacements);
//Real t_sum = std::abs(SurfaceStatistics::computeSum(surface_tractions));
//Real s_max = SurfaceStatistics::computeMaximum(surface);
#pragma omp parallel for reduction(+:res)
for (unsigned int i = 0; i < size; ++i) {
res +=
this->surface_tractions(i)
*(surface_displacements(i) - surface(i));
}
STOPTIMER("computeF");
this->computed_oldF = -1.0*res;
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::computeGaps() {
int size = surface.size();
for (int 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;
// choose a point in contact
//UInt index = 0;
//Real max_pressure = 0;
unsigned int n = surface.size();
unsigned int size = n*n;
// for (unsigned int i = 0; i < size; ++i) {
// if (surface_tractions(i).real() > max_pressure) {
// index = i;
// max_pressure = surface_tractions(i).real();
// }
// }
// Real shift = surface_displacements[index].real() - surface[index].real();
Real shift = 0;
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;
}
}
/* -------------------------------------------------------------------------- */
c4science · Help