diff --git a/src/bem_fft_base.cpp b/src/bem_fft_base.cpp
index 2e2893e0..77506ada 100644
--- a/src/bem_fft_base.cpp
+++ b/src/bem_fft_base.cpp
@@ -1,301 +1,299 @@
/**
*
* @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 "elastic_energy_functional.hh"
/* -------------------------------------------------------------------------- */
#include <omp.h>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
BemFFTBase::BemFFTBase(Surface<Real> & p):
BemInterface(p),
functional(new ElasticEnergyFunctional(*this)),
surface_tractions(p.size(),p.getL()),
surface_displacements(p.size(),p.getL()),
surface_spectral_tractions(p.size(),p.getL()),
surface_spectral_displacements(p.size(),p.getL()),
pressure_transform(p.size(), surface_tractions, surface_spectral_tractions),
displacement_transform(p.size(), surface_displacements, surface_spectral_displacements),
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;
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;
}
/* -------------------------------------------------------------------------- */
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;
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 (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::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");
pressure_transform.forwardTransform();
STOPTIMER("computeSpectralTraction");
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::computeDisplacements(){
STARTTIMER("computeDisplacement");
this->computeSpectralTractions();
this->computeSpectralDisplacements();
this->computeRealDisplacementsFromSpectral();
STOPTIMER("computeDisplacement");
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::computeRealTractionsFromSpectral(){
STARTTIMER("computeRealTractionsFromSpectral");
pressure_transform.backwardTransform();
STOPTIMER("computeRealDisplacementFromSpectral");
}
/* -------------------------------------------------------------------------- */
void BemFFTBase::computeRealDisplacementsFromSpectral(){
STARTTIMER("computeRealDisplacementFromSpectral");
displacement_transform.backwardTransform();
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);
}
}
/* -------------------------------------------------------------------------- */
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::setFunctional(Functional *new_funct) {
delete this->functional;
this->functional = new_funct;
}
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
diff --git a/src/bem_polonski.cpp b/src/bem_polonski.cpp
index 22ab9c88..d7b6c343 100644
--- a/src/bem_polonski.cpp
+++ b/src/bem_polonski.cpp
@@ -1,309 +1,309 @@
/**
*
* @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 <vector>
#include "surface.hh"
#include "bem_polonski.hh"
#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <cmath>
/* -------------------------------------------------------------------------- */
#define TIMER
#include "surface_timer.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
BemPolonski::BemPolonski(Surface<Real> & p) :
BemFFTBase(p),
surface_t(p.size(),p.getL()),
surface_r(p.size(),p.getL()),
surface_spectral_r(p.size(),p.getL()),
surface_r_transform(p.size(), surface_r, surface_spectral_r),
pold(p.size(),p.getL()){
e_star = 1.;
max_iterations = 10000;
}
BemPolonski::~BemPolonski() {}
Real BemPolonski::computeEquilibrium(Real epsilon,
Real pressure){
this->computeSpectralInfluenceOverDisplacement();
this->surface_t = 0.;
this->surface_r = 0.;
this->surface_spectral_r = 0.;
this->pold = 0.;
this->surface_displacements = 0.;
Real delta = 0.;
Real Gold = 1.;
Real f = 1e300;
Real fPrevious = 1e300;
Real current_pressure = SurfaceStatistics::computeAverage(surface_tractions);
if (current_pressure <= 0.) surface_tractions = pressure;
this->enforcePressureBalance(pressure);
convergence_iterations.clear();
nb_iterations = 0;
while(f > epsilon && nb_iterations < max_iterations) {
fPrevious = f;
this->computeDisplacements();
this->computeGaps();
Real gbar = this->computeMeanGapsInContact();
this->gap -= gbar;
Real G = this->computeG();
this->updateT(G,Gold,delta);
Real tau = this->computeTau();
pold = this->surface_tractions;
Gold = G;
delta = this->updateTractions(tau);
this->enforcePressureBalance(pressure);
f = this->computeF();
if (nb_iterations % dump_freq == 0){
Real A = SurfaceStatistics::computeContactAreaRatio(surface_tractions);
std::cout << std::scientific << std::setprecision(10)
<< nb_iterations << " "
<< f << " " << f-fPrevious << " " << A
<< std::endl;
}
convergence_iterations.push_back(f);
++nb_iterations;
}
return f;
}
/* -------------------------------------------------------------------------- */
void BemPolonski::computeGaps() {
UInt size = surface.size();
#pragma omp parallel for
for (UInt i=0; i < size*size; i++) {
gap(i) = surface_displacements(i) - surface(i);
}
}
/* -------------------------------------------------------------------------- */
Real BemPolonski::computeMeanGapsInContact(){
STARTTIMER("computeMeanGapsInContact");
UInt n = surface.size();
UInt size = n*n;
Real res = 0.;
UInt nb_contact = 0;
#pragma omp parallel for reduction(+: nb_contact,res)
for (UInt i = 0; i < size; ++i) {
if (this->surface_tractions(i) <= 0) continue;
++nb_contact;
res += this->gap(i);
}
res /= nb_contact;
STOPTIMER("computeMeanGapsInContact");
return res;
}
/* -------------------------------------------------------------------------- */
Real BemPolonski::computeG(){
STARTTIMER("computeG");
UInt n = surface.size();
UInt size = n*n;
Real res = 0.;
#pragma omp parallel for reduction(+: res)
for (UInt i = 0; i < size; ++i) {
if (this->surface_tractions(i) <= 0) continue;
Real val = this->gap(i);
res += val*val;
}
STOPTIMER("computeG");
return res;
}
/* -------------------------------------------------------------------------- */
void BemPolonski::updateT(Real G, Real Gold, Real delta){
STARTTIMER("updateT");
UInt n = surface.size();
UInt size = n*n;
Real factor = delta*G/Gold;
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
if (this->surface_tractions(i) <= 0) this->surface_t(i) = 0.;
else{
this->surface_t(i) *= factor;
this->surface_t(i) += this->gap(i);
}
}
STOPTIMER("updateT");
}
/* -------------------------------------------------------------------------- */
Real BemPolonski::computeTau(){
STARTTIMER("computeTau");
- surface_spectral_r = surface_t;
+ surface_r = surface_t;
surface_r_transform.forwardTransform();
UInt n = surface.size();
UInt size = n*n;
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
surface_spectral_r(i) *= this->surface_spectral_influence_disp(i);
}
surface_r_transform.backwardTransform();
Real rbar = 0;
UInt nb_contact = 0;
#pragma omp parallel for reduction(+: nb_contact,rbar)
for (UInt i = 0; i < size; ++i) {
if (this->surface_tractions(i) <= 0) continue;
++nb_contact;
rbar += surface_r(i);
}
rbar /= nb_contact;
surface_r -= rbar;
Real tau_sum1 = 0.,tau_sum2 = 0.;
#pragma omp parallel for reduction(+: tau_sum1, tau_sum2)
for (UInt i = 0; i < size; ++i) {
if (this->surface_tractions(i) <= 0) continue;
tau_sum1 += this->gap(i)*surface_t(i);
tau_sum2 += surface_r(i)*surface_t(i);
}
Real tau = tau_sum1/tau_sum2;
STOPTIMER("computeTau");
return tau;
}
/* -------------------------------------------------------------------------- */
Real BemPolonski::updateTractions(Real tau){
STARTTIMER("updateTractions");
UInt n = surface.size();
UInt size = n*n;
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
if (this->surface_tractions(i) <= 0) continue;
this->surface_tractions(i) -= tau*this->surface_t(i);
if (this->surface_tractions(i) < 0.) this->surface_tractions(i) = 0;
}
//compute number of interpenetration without contact
UInt nb_iol = 0;
#pragma omp parallel for reduction(+: nb_iol)
for (UInt i = 0; i < size; ++i) {
if (this->surface_tractions(i) <= 0. &&
this->gap(i) < 0.){
this->surface_tractions(i) -= tau*this->gap(i);
++nb_iol;
}
}
Real delta = 0;
if (nb_iol > 0) delta = 0.;
else delta = 1.;
STOPTIMER("updateTractions");
return delta;
}
/* -------------------------------------------------------------------------- */
void BemPolonski::enforcePressureBalance(Real applied_pressure){
STARTTIMER("enforcePressureBalance");
UInt n = surface.size();
UInt size = n*n;
Real pressure = 0.;
#pragma omp parallel for reduction(+: pressure)
for (UInt i = 0; i < size; ++i) {
pressure += this->surface_tractions(i);
}
pressure *= 1./size;
for (UInt i = 0; i < size; ++i) {
this->surface_tractions(i) *= applied_pressure/pressure;
}
STOPTIMER("enforcePressureBalance");
}
/* -------------------------------------------------------------------------- */
Real BemPolonski::computeF() {
UInt n = surface.size();
UInt 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 (UInt 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());
}
return res/std::abs(t_sum);
}
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
diff --git a/src/fftransform.cpp b/src/fftransform.cpp
index 553aec9b..0ee0c942 100644
--- a/src/fftransform.cpp
+++ b/src/fftransform.cpp
@@ -1,130 +1,154 @@
/**
*
* @author Lucas Frérot <lucas.frerot@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 "fftransform.hh"
__BEGIN_TAMAAS__
template<typename T>
FFTransform<T>::FFTransform(UInt size, Surface<T> & real, SurfaceComplex<T> & spectral):
real(real),
spectral(spectral),
forward_plan(NULL),
backward_plan(NULL),
size(size),
reduced_size(size / 2 + 1)
{
spectral_buffer = new fftw_complex[size * reduced_size];
Real * in = const_cast<Real *>(real.getInternalData());
forward_plan = fftw_plan_dft_r2c_2d(size, size,
in,
spectral_buffer,
- FFTW_MEASURE);
+ FFTW_ESTIMATE);
backward_plan = fftw_plan_dft_c2r_2d(size, size,
spectral_buffer,
in,
- FFTW_MEASURE);
+ FFTW_ESTIMATE);
}
/* -------------------------------------------------------------------------- */
template<typename T>
FFTransform<T>::~FFTransform() {
fftw_destroy_plan(forward_plan);
fftw_destroy_plan(backward_plan);
delete[] spectral_buffer;
}
/* -------------------------------------------------------------------------- */
template<typename T>
void FFTransform<T>::forwardTransform() {
fftw_execute(forward_plan);
writeBuffer();
}
/* -------------------------------------------------------------------------- */
template<typename T>
void FFTransform<T>::backwardTransform() {
loadBuffer();
fftw_execute(backward_plan);
normalize();
}
/* -------------------------------------------------------------------------- */
template<typename T>
void FFTransform<T>::normalize() {
#pragma omp parallel for
for (UInt i = 0 ; i < size * size ; ++i)
real(i) /= size * size;
}
/* -------------------------------------------------------------------------- */
template <typename T>
void FFTransform<T>::writeBuffer() {
-#pragma omp parallel for
- for (UInt i = 0 ; i < size ; i++) {
+#pragma omp parallel
+{
+#pragma omp for
+ // First line of spectral surface is different
+ for (UInt j = 1 ; j < reduced_size-1 ; j++) {
+ std::complex<T> z(spectral_buffer[j][0],
+ spectral_buffer[j][1]);
+ spectral(0, j) = z;
+ spectral(0, size-j) = std::conj(z);
+ }
+
+#pragma omp for
+ // First column just needs copying
+ for (UInt i = 0 ; i < size ; ++i) {
+ UInt index = i * reduced_size;
+ std::complex<T> z(spectral_buffer[index][0],
+ spectral_buffer[index][1]);
+ spectral(i, 0) = z;
+ }
+
+#pragma omp for
+ // Copy reduced_size-1 column
+ for (UInt i = 0 ; i < size ; ++i) {
+ UInt index = i * reduced_size + reduced_size-1;
+ std::complex<T> z(spectral_buffer[index][0],
+ spectral_buffer[index][1]);
+ spectral(i, reduced_size-1) = z;
+ }
+
+#pragma omp for collapse(2)
+ // Rest of spectral surface is entertwined
+ for (UInt i = 1 ; i < size ; i++) {
for (UInt j = 1 ; j < reduced_size-1 ; j++) {
UInt index = i * reduced_size + j;
std::complex<T> z(spectral_buffer[index][0],
spectral_buffer[index][1]);
spectral(i, j) = z;
- spectral(i, size-j) = std::conj(z);
+ spectral(size-i, size-j) = std::conj(z);
}
}
-#pragma omp parallel for
- for (UInt i = 0 ; i < size ; i++) {
- UInt index = i * reduced_size;
- std::complex<T> z(spectral_buffer[index][0],
- spectral_buffer[index][1]);
- spectral(i, 0) = z;
- }
+} // end #pramga omp parallel
}
/* -------------------------------------------------------------------------- */
template <typename T>
void FFTransform<T>::loadBuffer() {
#pragma omp parallel for collapse(2)
for (UInt i = 0 ; i < size ; i++) {
for (UInt j = 0 ; j < reduced_size ; j++) {
UInt index = i * reduced_size + j;
spectral_buffer[index][0] = spectral(i, j).real();
spectral_buffer[index][1] = spectral(i, j).imag();
}
}
}
/* -------------------------------------------------------------------------- */
template class FFTransform<Real>;
__END_TAMAAS__