diff --git a/.gitignore b/.gitignore
index 9460480f..7f247251 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,6 @@
*~
.sconf_temp
.sconsign.dblite
*.log
build*
-site_scons/site_tools/swig.pyc
+*.pyc
diff --git a/SConstruct b/SConstruct
index bf278af6..d9ad04c4 100644
--- a/SConstruct
+++ b/SConstruct
@@ -1,107 +1,108 @@
import os,sys
#path = '/usr/lib/scons/'
#sys.path.append(path)
#from SCons import *
#from SCons.Environment import *
#from SCons.Variables import *
#from SCons.Util import *
#script_dir = GetOption('file')
#if script_dir == []:
# script_dir = '.'
#else:
# script_dir = script_dir[0]
colors = {}
colors['cyan'] = ''
colors['purple'] = ''
colors['blue'] = ''
colors['green'] = ''
colors['yellow'] = ''
colors['red'] = ''
colors['end'] = ''
#colors['cyan'] = '\033[96m'
#colors['purple'] = '\033[95m'
#colors['blue'] = '\033[94m'
#colors['green'] = '\033[92m'
#colors['yellow'] = '\033[93m'
#colors['red'] = '\033[91m'
#colors['end'] = '\033[0m'
main_env = Environment(CXX= 'g++', ENV = os.environ)
vars = Variables('build-setup.conf')
vars.Add(EnumVariable('build_type', 'Build type', 'release',
allowed_values=('release', 'profiling', 'debug'),
ignorecase=2))
vars.Add('prefix', 'Prefix where to install', '/usr/local')
vars.Add(BoolVariable('timer', 'Activate the timer possibilities', False))
vars.Add(BoolVariable('verbose', 'Activate verbosity', False))
vars.Update(main_env)
Help(vars.GenerateHelpText(main_env))
build_type = main_env['build_type']
build_dir = 'build-' + main_env['build_type']
print "Building in ", build_dir
timer_flag = main_env['timer']
verbose = main_env['verbose']
if not verbose:
main_env['SHCXXCOMSTR'] = unicode('{0}[Compiling] {1}$SOURCE{2}'.format(colors['green'],colors['blue'],colors['end']))
main_env['SHLINKCOMSTR'] = unicode('{0}[Linking] {1}$TARGET{2}'.format(colors['purple'],colors['blue'],colors['end']))
main_env['SWIGCOMSTR'] = unicode('{0}[Swig] {1}$SOURCE{2}'.format(colors['yellow'],colors['blue'],colors['end']))
main_env.AppendUnique(CPPPATH=['#/src', '#/python'])
main_env.AppendUnique(LIBS='gomp')
# Flags and options
main_env.AppendUnique(CXXFLAGS=['-std=c++11',
'-Wall',
'-fopenmp'])
main_env.AppendUnique(SHLINKFLAGS=['-fopenmp'])
if timer_flag:
main_env.AppendUnique(CPPDEFINES=['USING_TIMER'])
cxxflags_dict = {
"debug":Split("-g -O0"),
"profiling":Split("-g -pg -O2"),
"release":Split("-O3")
}
shlinkflags_dict = {
"debug":[],
"profiling":['-pg'],
"release":[]
}
main_env.AppendUnique(CXXFLAGS=cxxflags_dict[build_type])
main_env.AppendUnique(SHLINKFLAGS=shlinkflags_dict[build_type])
main_env['LIBPATH'] = [os.path.abspath(os.path.join(build_dir,'src'))]
main_env['RPATH'] = "$LIBPATH"
# save new values
vars.Save('build-setup.conf', main_env)
Export('main_env')
SConscript('src/SConscript',variant_dir=os.path.join(build_dir,'src'),duplicate=False)
SConscript('python/SConscript',variant_dir=os.path.join(build_dir,'python'),duplicate=False)
+SConscript('tests/SConscript',variant_dir=os.path.join(build_dir,'tests'),duplicate=False)
# saving the env file
try:
env_file = open(os.path.join(build_dir,'tamaas_environement.sh'),'w')
env_file.write("""
export PYTHONPATH=$PYTHONPATH:{0}/python
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:{0}/src
""".format(os.path.abspath(build_dir)))
env_file.close()
except:
pass
diff --git a/python-examples/hertz.py b/python-examples/hertz.py
index b84716bd..ef3c2640 100644
--- a/python-examples/hertz.py
+++ b/python-examples/hertz.py
@@ -1,146 +1,123 @@
#!/usr/bin/python
# -*- coding: utf-8 -*-
##*
#
# @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/>.
#
#
################################################################
import tamaas
def dumpHertzPrediction(file, bem, pressure, r, E ):
A = tamaas.SurfaceStatistics.computeContactArea(bem.getTractions())
Aratio = tamaas.SurfaceStatistics.computeContactAreaRatio(bem.getTractions())
radius = np.sqrt(A/np.pi)
#L = bem.getSurface().getL()
L = 1.
load = pressure * L*L
radius_hertz = pow(0.75*load*r/E,1./3.)
p0_hertz = 1./np.pi*pow(6.*E*E*load/r/r,1./3.)
p0 = tamaas.SurfaceStatistics.computeMaximum(bem.getTractions())
n = bem.getTractions().shape[0]
computed_load = bem.getTractions().sum()/n/n*L*L
file.write("{0}\t{1}\t{2}\t{3}\t{4}\t{5}\t{6}\t{7}\t{8}\t{9}\n".format(pressure,load,computed_load,Aratio,A,radius,radius_hertz,radius_hertz/radius,p0,p0_hertz))
################################################################
def main(argv):
parser = argparse.ArgumentParser(description='Hertz test for the Boundary element method of Stanley and Kato')
parser.add_argument("--N",type=int, help="Surface size", required=True)
parser.add_argument("--r",type=float, help="radius of hertz sphere", required=True)
- parser.add_argument("--steps",type=int, help="number of steps within which the pressure is applied.", required=True)
+ parser.add_argument("--steps",type=int, help="number of steps within which the pressure is applied.", required=True)
parser.add_argument("--precision",type=float, help="relative precision, convergence if | (f_i - f_{i-1})/f_i | < Precision.", required=True)
parser.add_argument("--e_star",type=float, help="effective elastic modulus" , required=True)
parser.add_argument("--L",type=float, help="size of the surface" , required=True)
parser.add_argument("--max_pressure",type=float, help="maximal load requested" , required=True)
parser.add_argument("--plot_surface",type=bool, help="request output of text files containing the contact pressures on the surface" , default=False)
parser.add_argument("--nthreads",type=int, help="request a number of threads to use via openmp to compute" , default=1)
args = parser.parse_args()
arguments = vars(args)
n = arguments["N"]
r = arguments["r"]
max_pressure = arguments["max_pressure"]
Ninc = arguments["steps"]
epsilon = arguments["precision"]
Estar = arguments["e_star"]
L = arguments["L"]
plot_surface = arguments["plot_surface"]
nthreads = arguments["nthreads"]
pressure = 0.
dp = max_pressure/float(Ninc)
s = np.zeros((n,n))
print s.shape
for i in range(0,n):
for j in range(0,n):
x = 1.*i - n/2
y = 1.*j - n/2
d = (x*x + y*y)*1./n/n*L*L
if d < r*r: s[i,j] = - r + np.sqrt(r*r-d)
else: s[i,j] = - r
print "\n::: DATA ::::::::::::::::::::::::::::::\n"
print " [N] {0}\n".format(n)
print " [r] {0}\n".format(r)
print " [Pext] {0}\n".format(max_pressure)
print " [Steps] {0}\n".format(Ninc)
print " [Precision] {0}\n".format(epsilon)
bem = tamaas.BemPolonski(s)
bem.setEffectiveModulus(Estar)
file = open("hertz-prediction",'w')
file.write("#pressure\tload\tcomputed-load\tarea_ratio\tarea\tradius\tradius-hertz\tradius/radius-hertz\tp0\tp0-hertz\n")
for i in range(0,Ninc):
pressure += dp
bem.computeEquilibrium(epsilon,pressure)
A = tamaas.SurfaceStatistics.computeContactAreaRatio(bem.getTractions())
-# std::stringstream fname
-# fname << n << "-" << r << "-"
-# << pressure
-# bem.getTractions().dumpToParaview(fname.str())
-# fname << "-disp"
-# bem.getTrueDisplacements().dumpToParaview(fname.str())
-#
-# std::cout << "* " << i << "/" << Ninc << ", pressure = " << pressure << " area " << A << std::endl
dumpHertzPrediction(file,bem,pressure,r,Estar)
-# if (plot_surface) {
-# std::stringstream prefix
-# prefix << "-A-" << A
-# << "-pressure-" << pressure
-# << "-step-" << i
-#
-# std::string filename = parser.makeFilename("pressures" + prefix.str(),"plot")
-# bem.getTractions().dumpToTextFile<false>(filename)
-# std::string filename_disp = parser.makeFilename("displacement"+ prefix.str(),"plot")
-# bem.getDisplacements().dumpToTextFile<false>(filename_disp)
-# std::string filename_gap = parser.makeFilename("gap"+ prefix.str(),"plot")
-# bem.getGap().dumpToTextFile<false>(filename_gap)
-# }
-
-
if A == 1.0:
file.close()
break
tamaas.dumpTimes()
################################################################
import sys
import argparse
import numpy as np
main(sys.argv)
diff --git a/src/SConscript b/src/SConscript
index 004baf7d..f8589c76 100644
--- a/src/SConscript
+++ b/src/SConscript
@@ -1,79 +1,80 @@
Import('main_env')
import distutils.sysconfig
import os
fftw_include = ""
fftw_library = ""
if 'FFTW_INCLUDE' in main_env['ENV']:
fftw_include = main_env['ENV']['FFTW_INCLUDE']
if 'FFTW_LIBRARY' in main_env['ENV']:
fftw_library = main_env['ENV']['FFTW_LIBRARY']
env = main_env.Clone(
tools = [fftw],
FFTW_LIBRARY_WISH = ['main', 'omp'],
FFTW_INCLUDE_DIR = fftw_include,
FFTW_LIBRARY_DIR = fftw_library,
)
# Lib roughcontact
generator_list = """
surface_generator.cpp
surface_generator_voss.cpp
surface_generator_crenel.cpp
surface_generator_ellipsoid.cpp
surface_generator_filter.cpp
surface_generator_filter_bessel.cpp
surface_generator_filter_fft.cpp
""".split()
#env.SharedLibrary('Generator',generator_list)
# Lib SURFACE
surface_list = """
map_2d.cpp
map_2d_square.cpp
surface.cpp
surface_timer.cpp
tamaas.cpp
""".split()
#env.SharedLibrary('Surface',surface_list)
# Lib PERCOLATION
percolation_list = """
cluster_grow.cpp
contact_area.cpp
contact_cluster.cpp
contact_cluster_collection.cpp
""".split()
#env.SharedLibrary('Percolation',percolation_list)
# BEM PERCOLATION
bem_list = """
bem_kato.cpp
bem_polonski.cpp
bem_gigi.cpp
bem_gigipol.cpp
bem_penalty.cpp
bem_uzawa.cpp
bem_fft_base.cpp
functional.cpp
elastic_energy_functional.cpp
adhesion_functional.cpp
fftransform.cpp
+fft_plan_manager.cpp
""".split()
#env.SharedLibrary('BEM',bem_list)
rough_contact_list = generator_list + surface_list + percolation_list + bem_list
current_files = set(os.listdir(Dir('.').srcnode().abspath))
unexpected_files = list(current_files - set(rough_contact_list + ['SConscript']))
unexpected_headers = [ f for f in unexpected_files if os.path.splitext(f)[1] == '.hh' ]
unexpected_files = [ f for f in unexpected_files if not os.path.splitext(f)[1] == '.hh' ]
#print "Unexpected files : {}".format(unexpected_files)
#print "Include path : {}".format([str(x) for x in env['CPPPATH']])
#print unexpected_headers
env.SharedLibrary('Tamaas',rough_contact_list)
diff --git a/src/bem_fft_base.cpp b/src/bem_fft_base.cpp
index 31931797..8638fe37 100644
--- a/src/bem_fft_base.cpp
+++ b/src/bem_fft_base.cpp
@@ -1,305 +1,284 @@
/**
*
* @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 "fft_plan_manager.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_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;
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;
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::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;
+void BemFFTBase::computeDisplacementsFromTractions(){
- 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;
+ STARTTIMER("computeDisplacementsFromTractions");
-#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;
+ this->applyInfluenceFunctions(surface_tractions,surface_displacements);
-
- STOPTIMER("computeSpectralDisplacement");
+ STOPTIMER("computeDisplacementFromTractions");
}
/* -------------------------------------------------------------------------- */
-
-void BemFFTBase::computeSpectralTractions(){
-
- STARTTIMER("computeSpectralTraction");
-
- pressure_transform.forwardTransform();
-
- STOPTIMER("computeSpectralTraction");
-
+void BemFFTBase::computeTractionsFromDisplacements(){
+
+ STARTTIMER("computeDisplacementsFromTractions");
+
+ this->applyInverseInfluenceFunctions(surface_displacements,surface_tractions);
+
+ STOPTIMER("computeDisplacementFromTractions");
+
}
/* -------------------------------------------------------------------------- */
+void BemFFTBase::applyInfluenceFunctions(Surface<Real> & input, Surface<Real> & output){
-void BemFFTBase::computeDisplacements(){
+ STARTTIMER("applyInfluenceFunctions");
- STARTTIMER("computeDisplacement");
+ FFTransform<Real> & plan1 = FFTPlanManager::get().createPlan(input,surface_spectral_output);
+
+ plan1.forwardTransform();
- this->computeSpectralTractions();
- this->computeSpectralDisplacements();
- this->computeRealDisplacementsFromSpectral();
+ surface_spectral_output *= surface_spectral_influence_disp;
- STOPTIMER("computeDisplacement");
+ FFTransform<Real> & plan2 = FFTPlanManager::get().createPlan(output,surface_spectral_output);
-}
+ plan2.backwardTransform();
-/* -------------------------------------------------------------------------- */
-
-void BemFFTBase::computeRealTractionsFromSpectral(){
+ STOPTIMER("applyInfluenceFunctions");
- STARTTIMER("computeRealTractionsFromSpectral");
+}
- pressure_transform.backwardTransform();
+/* -------------------------------------------------------------------------- */
- STOPTIMER("computeRealDisplacementFromSpectral");
+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;
-void BemFFTBase::computeRealDisplacementsFromSpectral(){
+ 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);
+ }
- STARTTIMER("computeRealDisplacementFromSpectral");
+ FFTransform<Real> & plan2 = FFTPlanManager::get().createPlan(output,surface_spectral_output);
- displacement_transform.backwardTransform();
+ plan2.backwardTransform();
- STOPTIMER("computeRealDisplacementFromSpectral");
+ 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) = 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::setFunctional(Functional *new_funct) {
delete this->functional;
this->functional = new_funct;
}
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
diff --git a/src/bem_fft_base.hh b/src/bem_fft_base.hh
index fe111cfd..fa2d3812 100644
--- a/src/bem_fft_base.hh
+++ b/src/bem_fft_base.hh
@@ -1,167 +1,157 @@
/**
*
* @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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef BEM_FFT_BASE_H
#define BEM_FFT_BASE_H
/* -------------------------------------------------------------------------- */
#include "bem_interface.hh"
#include <cmath>
#include "surface_statistics.hh"
#include "fftransform.hh"
#include "functional.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
class BemFFTBase : public BemInterface {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
BemFFTBase(Surface<Real> & p);
virtual ~BemFFTBase();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
+ //! compute the displacement from traction input
+ virtual void computeDisplacementsFromTractions();
+ //! compute the tractions from displacement input
+ virtual void computeTractionsFromDisplacements();
+ //! apply influence funtions
+ virtual void applyInfluenceFunctions(Surface<Real> & input, Surface<Real> & output);
+ //! apply inverse influence funtions
+ virtual void applyInverseInfluenceFunctions(Surface<Real> & input, Surface<Real> & output);
//! compute the influcence coefficients of pressure over displacement
virtual void computeSpectralInfluenceOverDisplacement();
- //! compute real tractions from fourier represnetation (IFFT)
- virtual void computeRealTractionsFromSpectral();
- //! compute real displacements from fourier represnetation (IFFT)
- virtual void computeRealDisplacementsFromSpectral();
- //! compute the displacement in real space
- virtual void computeDisplacements();
- //! compute the displacement fourier space
- virtual void computeSpectralDisplacements();
- //! compute the tractions in fourier space
- virtual void computeSpectralTractions();
- //! compute the gradient the objective function (with respect to traction)
+ //! compute the displacements
virtual void computeTrueDisplacements();
//! compute the gap
virtual void computeGaps();
//! get the number of iterations
virtual UInt getNbIterations() const = 0;
//! get the convergence steps of iterations
virtual const std::vector<Real> & getConvergenceIterations() const = 0;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
//! return the influcence coefficients
virtual const Surface<Real> & getSpectralInfluenceOverDisplacement();
//! return the actual tractions
virtual const Surface<Real> & getTractions() const;
//! return the actual displacements
virtual const Surface<Real> & getDisplacements()const;
- //! return the actual displacements
- virtual const SurfaceComplex<Real> & getSpectralDisplacements()const;
- //! return the actual displacements
- virtual const SurfaceComplex<Real> & getSpectralTractions()const;
//! get the computed true displacement
const Surface<Real> & getTrueDisplacements(){return true_displacements;};
//! get the computed gap
const Surface<Real> & getGap(){return gap;};
//! set the number of threads to be used to perform the calculations
void setNumberOfThreads(UInt nthreads);
//! get the number of threads
UInt getNumberOfThreads() const { return nthreads; }
//! set the elastic effective modulus
void setEffectiveModulus(Real e_star);
//! set the maximal number of iterations
void setMaxIterations(UInt max_iter){this->max_iterations = max_iter;};
//! set dump freq
void setDumpFreq(Real dump_freq){this->dump_freq = dump_freq;};
//! set functional
void setFunctional(::tamaas::Functional * new_funct);
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
//! Functional for optimization
Functional * functional;
//! surface tractions
Surface<Real> surface_tractions;
//! surface displacement
Surface<Real> surface_displacements;
- //! surface tractions
- SurfaceComplex<Real> surface_spectral_tractions;
- //! surface displacement
- SurfaceComplex<Real> surface_spectral_displacements;
- //! FF Transform for pressures
- FFTransform<Real> pressure_transform;
- //! FF Transform for displacements
- FFTransform<Real> displacement_transform;
- //! surface influcence coefficient over displacement
+ //! surface input in spectral
+ SurfaceComplex<Real> surface_spectral_input;
+ //! surface output in spectral
+ SurfaceComplex<Real> surface_spectral_output;
+ //! surface influence coefficient over displacement
Surface<Real> surface_spectral_influence_disp;
//! shifted displacement in order to get the true displacement
Surface<Real> true_displacements;
//! old displacements (for stopping criterion)
Surface<Real> old_displacements;
//! computing of the gap
Surface<Real> gap;
//! Effective young modulus
Real e_star;
//! Number of threads to use during calculations
UInt nthreads;
//! maximum of iterations
UInt max_iterations;
//! number of iterations
UInt nb_iterations;
//! number of iterations
std::vector<Real> convergence_iterations;
//! frequency at which to dump iteration info
UInt dump_freq;
/* ------------------------------------------------------------------------ */
/* temporary trick about some obscure scaling factor */
/* ------------------------------------------------------------------------ */
/*
0 2 = 2*1
1 4 = 2*2
2 10 = 2*5
3 20 = 2*10
4 34 = 2*12
*/
//#define tmp_coeff (2.*n/M_PI)
//#define tmp_coeff (n/M_PI)
const Real tmp_coeff;
};
__END_TAMAAS__
#endif /* BEM_FFT_BASE_H */
diff --git a/src/bem_gigi.cpp b/src/bem_gigi.cpp
index b29fc5a8..b62f5d51 100644
--- a/src/bem_gigi.cpp
+++ b/src/bem_gigi.cpp
@@ -1,356 +1,339 @@
/**
*
* @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_gigi.hh"
#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <cmath>
/* -------------------------------------------------------------------------- */
#define TIMER
#include "surface_timer.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
Real BemGigi::computeEquilibrium(Real epsilon,
Real mean_displacement) {
this->computeSpectralInfluenceOverDisplacement();
Real Rold = 1.;
Real f = 1e300;
this->search_direction = 0.;
this->true_displacements =0;
this->true_displacements = mean_displacement;
this->computeGaps();
this->optimizeToMeanDisplacement(mean_displacement);
this->computeGaps();
convergence_iterations.clear();
nb_iterations = 0;
max_iterations =50;
Real R=0;
while (f> epsilon && nb_iterations++ < max_iterations) {
this->computeGaps();
this->functional->computeGradF();
this->updateSearchDirection(R / Rold);
Rold = R;
Real alpha = this->computeOptimalStep();
std::cout << "alpha vaut "<< alpha<< std::endl;
this->old_displacements=this->true_displacements;
this->updateDisplacements(alpha);
this->computeGaps();
this->optimizeToMeanDisplacement(mean_displacement); //espace admissible
this->computeGaps();
f=computeStoppingCriterion();
}
this->computePressures(mean_displacement);
return f;
}
/* -------------------------------------------------------------------------- */
Real BemGigi::computeStoppingCriterion() {
Real crit=0.;
Real disp_norm = 0.;
UInt n = surface.size();
UInt size = n*n;
#pragma omp parallel for reduction(+:crit, disp_norm)
for (UInt i = 0; i < size; ++i) {
crit += (true_displacements(i)-old_displacements(i))*(true_displacements(i)-old_displacements(i));
disp_norm += (true_displacements(i)*true_displacements(i));
}
return fabs(crit) / disp_norm;
}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
void BemGigi::optimizeToMeanDisplacement(Real imposed_mean) {
Real target_value = imposed_mean - SurfaceStatistics::computeAverage(surface, 0);
UInt n = surface.size();
UInt size = n * n;
// Initial guesses for upper and lower bound
Real step_min = -10;
Real step_max = 10;
// Gaps for upper and lower bound
Real gap_min = positiveGapAverage(step_min);
Real gap_max = positiveGapAverage(step_max);
UInt max_expansion = 8;
// Expand bounds if necessary
for (UInt i = 0 ; gap_max < target_value && i < max_expansion ; i++, step_max *= 10)
gap_max = positiveGapAverage(step_max);
for (UInt i = 0 ; gap_min > target_value && i < max_expansion ; i++, step_min *= 10)
gap_min = positiveGapAverage(step_min);
Real g = 0.;
Real epsilon = 1e-12;
Real step = 0;
while (fabs(g - target_value) > epsilon) {
step = (step_min + step_max) / 2.;
g = positiveGapAverage(step);
if (g > target_value) step_max = step;
else if (g < target_value) step_min = step;
else {
step_max = step;
step_min = step;
}
}
step = (step_min + step_max) / 2.;
#pragma omp parallel for
for (UInt i = 0 ; i < size ; i++) {
gap(i) += step;
if (gap(i) < 0) gap(i) = 0;
true_displacements(i) = gap(i) + surface(i);
}
}
/* -------------------------------------------------------------------------- */
Real BemGigi::positiveGapAverage(Real shift) {
UInt n = surface.size();
Real res = 0;
#pragma omp parallel for reduction(+: res)
for (UInt i = 0 ; i < n * n ; i++) {
Real shifted_gap = gap(i) + shift;
res += shifted_gap * (shifted_gap > 0);
}
return res / (n * n);
}
/* -------------------------------------------------------------------------- */
void BemGigi::updateSearchDirection(Real factor) {
STARTTIMER("updateSearchDirection");
UInt n = surface.size();
UInt size = n*n;
Surface<Real> & gradF = this->functional->getGradF();
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
this->search_direction(i)= gradF(i);
}
STOPTIMER("updateSearchDirection");
}
/* -------------------------------------------------------------------------- */
Real BemGigi::computeOptimalStep() {
STARTTIMER("computeOptimalStep");
search_direction.FFTTransform(spectral_search_direction, nthreads);
UInt n = surface.size();
UInt size = n*n;
spectral_search_direction(0) = 0;
#pragma omp parallel for
for (UInt i = 1; i < size; ++i) {
spectral_search_direction(i) /= this->surface_spectral_influence_disp(i);
}
// /!\ does not contain the spectral search direction anymore
spectral_search_direction.FFTITransform(nthreads);
// Real average = SurfaceStatistics::computeAverage(spectral_search_direction.real(), 0);
// spectral_search_direction -= average;
Surface<Real> & gradF = this->functional->getGradF();
Real numerator = 0., denominator = 0.;
#pragma omp parallel for reduction(+: numerator, denominator)
for (UInt i = 0; i < size; ++i) {
numerator += gradF(i) * search_direction(i);
denominator += spectral_search_direction(i).real() * search_direction(i);
}
Real alpha = numerator / denominator;
STOPTIMER("computeOptimalStep");
return alpha;
}
/* -------------------------------------------------------------------------- */
Real BemGigi::computeTau() {
STARTTIMER("computeOptimalStep");
search_direction.FFTTransform(spectral_search_direction, nthreads);
UInt n = surface.size();
UInt size = n*n;
spectral_search_direction(0) = 0;
#pragma omp parallel for
for (UInt i = 1; i < size; ++i) {
spectral_search_direction(i) /= this->surface_spectral_influence_disp(i);
}
// /!\ does not contain the spectral search direction anymore
spectral_search_direction.FFTITransform(nthreads);
// Real average = SurfaceStatistics::computeAverage(spectral_search_direction.real(), 0);
// spectral_search_direction -= average;
Surface<Real> & gradF = this->functional->getGradF();
Real numerator = 0., denominator = 0.;
#pragma omp parallel for reduction(+: numerator, denominator)
for (UInt i = 0; i < size; ++i) {
if (this->gap(i) > 0) {
numerator += gradF(i) * search_direction(i);
denominator += spectral_search_direction(i).real() * search_direction(i);}
}
Real tau = numerator / denominator;
STOPTIMER("computeOptimalStep");
return tau;
}
/* -------------------------------------------------------------------------- */
void BemGigi::updateDisplacements(Real alpha) {
STARTTIMER("updateDisplacements");
UInt n = surface.size();
UInt size = n*n;
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
this->true_displacements(i) -= alpha*this->search_direction(i);
}
STOPTIMER("updateDisplacements");
}
/* -------------------------------------------------------------------------- */
void BemGigi::emptyOverlap() {
STARTTIMER("emptyoverlap");
UInt n = surface.size();
UInt size = n*n;
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
if (gap(i) < 0)
this->true_displacements(i) = this->surface(i);
}
STOPTIMER("emptyoverlap");
}
/* -------------------------------------------------------------------------- */
void BemGigi::enforceMeanDisplacement(Real mean_displacement) {
STARTTIMER("enforceMeanDisplacement");
UInt n = surface.size();
UInt size = n*n;
Real average = SurfaceStatistics::computeAverage(this->true_displacements, 0);
Real factor = mean_displacement / average;
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
this->true_displacements(i) *= factor;
}
STOPTIMER("enforceMeanDisplacement");
}
/* -------------------------------------------------------------------------- */
void BemGigi::computePressures(Real mean_displacement) {
- UInt n = surface.size();
- UInt size = n*n;
-
-
- this->true_displacements.FFTTransform(this->surface_spectral_tractions, nthreads);
- this->surface_spectral_tractions(0) = 0;
-
-#pragma omp parallel for
- for (UInt i = 1 ; i < size ; ++i) {
- this->surface_spectral_tractions(i) /= this->surface_spectral_influence_disp(i);
- }
-
- this->surface_spectral_tractions.FFTITransform(this->surface_tractions, nthreads);
-
-
-
-
-
+ this->computeTractionsFromDisplacements();
}
__END_TAMAAS__
diff --git a/src/bem_gigipol.cpp b/src/bem_gigipol.cpp
index 2ac91cf4..3456d7a4 100644
--- a/src/bem_gigipol.cpp
+++ b/src/bem_gigipol.cpp
@@ -1,415 +1,408 @@
/**
*
* @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_gigipol.hh"
#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <cmath>
/* -------------------------------------------------------------------------- */
#define TIMER
#include "surface_timer.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
Real BemGigipol::computeEquilibrium(Real epsilon,
Real mean_displacement) {
this->computeSpectralInfluenceOverDisplacement();
this->surface_t = 0.;
this->surface_r = 0.;
this->surface_spectral_r = 0.;
this->search_direction = 0.;
this->pold = 0.;
this->surface_displacements = 0.;
Real delta = 0.;
Real Gold = 1.;
Real f = 1e300;
Real fPrevious = 1e300;
UInt n = surface.size();
UInt size = n*n;
this->true_displacements =0;
this->true_displacements = mean_displacement;
this->computeGaps();
this->optimizeToMeanDisplacement(mean_displacement);
convergence_iterations.clear();
nb_iterations = 0;max_iterations=500;
while(f > epsilon && nb_iterations < max_iterations) {
fPrevious = f;
this->functional->computeGradF();
this->search_direction= this->functional->getGradF();
Real gbar = this->computeMeanPressuresInNonContact();
this->search_direction -= gbar;
Real G = this->computeG();
std::cout << "G vaut "<< G<< std::endl;
this->updateT(G,Gold,delta);
Real tau = this->computeTau();
this->old_displacements=this->true_displacements;
Gold = G;
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
this->true_displacements(i)-=tau*this->surface_t(i);
}
//Projection on admissible space
this->computeGaps();
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
if (this->gap(i)<0){
this->true_displacements(i)=this->surface(i);}
}
this->computeGaps();
delta = this->updateDisplacements(tau, mean_displacement);
// std::cout << "delta vaut "<< delta << std::endl;
//Projection on admissible space
this->computeGaps();
this->enforceMeanDisplacement(mean_displacement);
this->computeGaps();
this->computePressures(mean_displacement);
f=this->computeStoppingCriterion();
std::cout << "f vaut "<< f << std::endl;
if (nb_iterations % dump_freq ==1000){
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;
}
this->computePressures(mean_displacement);
return f;
}
/* -------------------------------------------------------------------------- */
Real BemGigipol::computeStoppingCriterion() {
Real crit=0.;
UInt n = surface.size();
UInt size = n*n;
#pragma omp parallel for reduction(+:crit)
for (UInt i = 0; i < size; ++i) {
crit += std::abs(this->gap(i)*(this->surface_tractions(i)));
;
}
return crit ;
}
/* -------------------------------------------------------------------------- */
Real BemGigipol::computeG(){
STARTTIMER("computeG");
unsigned int n = surface.size();
unsigned int size = n*n;
Real res = 0.;
#pragma omp parallel for reduction(+: res)
for (unsigned int i = 0; i < size; ++i) {
if (this->gap(i) > 0.) {
Real val = this->search_direction(i);
res += val*val;}
}
STOPTIMER("computeG");
return res;
}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
Real BemGigipol::computeMeanPressuresInNonContact(){
STARTTIMER("computeMeanPressuresInNonContact");
unsigned int n = surface.size();
unsigned int size = n*n;
Real res = 0.;
UInt nb_contact = 0;
#pragma omp parallel for reduction(+: nb_contact,res)
for (unsigned int i = 0; i < size; ++i) {
if (this->gap(i) == 0. ) continue;
++nb_contact;
res += this->search_direction(i);
}
res /= nb_contact;
STOPTIMER("computeMeanPressuresInNonContact");
return res;
}
/* -------------------------------------------------------------------------- */
void BemGigipol::optimizeToMeanDisplacement(Real imposed_mean) {
Real target_value = imposed_mean - SurfaceStatistics::computeAverage(surface, 0);
UInt n = surface.size();
UInt size = n * n;
// Initial guesses for upper and lower bound
Real step_min = -10;
Real step_max = 10;
// Gaps for upper and lower bound
Real gap_min = positiveGapAverage(step_min);
Real gap_max = positiveGapAverage(step_max);
UInt max_expansion = 8;
// Expand bounds if necessary
for (UInt i = 0 ; gap_max < target_value && i < max_expansion ; i++, step_max *= 10)
gap_max = positiveGapAverage(step_max);
for (UInt i = 0 ; gap_min > target_value && i < max_expansion ; i++, step_min *= 10)
gap_min = positiveGapAverage(step_min);
Real g = 0.;
Real epsilon = 1e-12;
Real step = 0;
while (fabs(g - target_value) > epsilon) {
step = (step_min + step_max) / 2.;
g = positiveGapAverage(step);
if (g > target_value) step_max = step;
else if (g < target_value) step_min = step;
else {
step_max = step;
step_min = step;
}
}
step = (step_min + step_max) / 2.;
#pragma omp parallel for
for (UInt i = 0 ; i < size ; i++) {
gap(i) += step;
if (gap(i) < 0) gap(i) = 0;
true_displacements(i) = gap(i) + surface(i);
}
}
/* -------------------------------------------------------------------------- */
Real BemGigipol::positiveGapAverage(Real shift) {
UInt n = surface.size();
Real res = 0;
#pragma omp parallel for reduction(+: res)
for (UInt i = 0 ; i < n * n ; i++) {
Real shifted_gap = gap(i) + shift;
res += shifted_gap * (shifted_gap > 0);
}
return res / (n * n);
}
/* -------------------------------------------------------------------------- */
void BemGigipol::updateT(Real G, Real Gold, Real delta){
STARTTIMER("updateT");
unsigned int n = surface.size();
unsigned int size = n*n;
Real factor = delta*G/Gold;
#pragma omp parallel for
for (unsigned int i = 0; i < size; ++i) {
if (this->gap(i) == 0.) this->surface_t(i) = 0.;
else{
this->surface_t(i) *= factor;
this->surface_t(i) += this->search_direction(i);
}
}
STOPTIMER("updateT");
}
/* -------------------------------------------------------------------------- */
Real BemGigipol::computeTau() {
STARTTIMER("computeOptimalStep");
surface_spectral_r = surface_t;
surface_spectral_r.FFTTransform(nthreads);
unsigned int n = surface.size();
unsigned int size = n*n;
surface_spectral_r(0)=0;
#pragma omp parallel for
for (unsigned int i = 1; i < size; ++i) {
surface_spectral_r(i) /= this->surface_spectral_influence_disp(i);
}
surface_spectral_r.FFTITransform(surface_r,nthreads);
Real rbar = 0;
UInt nb_contact = 0;
#pragma omp parallel for reduction(+: nb_contact,rbar)
for (unsigned int i = 0; i < size; ++i) {
if (this->gap(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 (unsigned int i = 0; i < size; ++i) {
if (this->gap(i) == 0) continue;
tau_sum1 += this->search_direction(i)*surface_t(i);
tau_sum2 += surface_r(i)*surface_t(i);
}
Real tau = tau_sum1/tau_sum2;
STOPTIMER("computeTau");
return tau;
}
/* -------------------------------------------------------------------------- */
Real BemGigipol::updateDisplacements(Real tau, Real mean_displacement) {
STARTTIMER("updateDisplacements");
unsigned int n = surface.size();
unsigned int size = n*n;
//compute number of interpenetration without contact
UInt nb_iol = 0;
#pragma omp parallel for reduction(+: nb_iol)
for (unsigned int i = 0; i < size; ++i) {
if ( this->search_direction(i) <0. && this->gap(i) == 0.){
this->true_displacements(i) -= tau* this->search_direction(i);
++nb_iol;
}
}
Real delta = 0;
if (nb_iol > 0) delta = 0.;
else delta = 1.;
return delta;
STOPTIMER("updateDisplacements");
}
/* -------------------------------------------------------------------------- */
void BemGigipol::enforceMeanDisplacement(Real mean_displacement) {
STARTTIMER("enforceMeanDisplacement");
UInt n = surface.size();
UInt size = n*n;
Real moyenne_surface=SurfaceStatistics::computeAverage(this->surface, 0);
Real average = SurfaceStatistics::computeAverage(this->true_displacements, 0);
Real factor = (mean_displacement-moyenne_surface) / (average-moyenne_surface);
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
this->true_displacements(i) = factor*(this->true_displacements(i)-this->surface(i)) +this->surface(i);
}
STOPTIMER("enforceMeanDisplacement");
}
/* -------------------------------------------------------------------------- */
void BemGigipol::computePressures(Real mean_displacement) {
UInt n = surface.size();
UInt size = n*n;
Real rho = this->functional->getParameter("rho");
Real surface_energy = this->functional->getParameter("surface_energy");
- this->true_displacements.FFTTransform(this->surface_spectral_tractions, nthreads);
- this->surface_spectral_tractions(0) = 0;
-
-#pragma omp parallel for
- for (UInt i = 1 ; i < size ; ++i) {
- this->surface_spectral_tractions(i) /= this->surface_spectral_influence_disp(i);
- }
-
- this->surface_spectral_tractions.FFTITransform(this->surface_tractions, nthreads);
-
+ this->applyInverseInfluenceFunctions(this->true_displacements, this->surface_tractions);
Real mini=3000;
for (UInt i = 0; i < size; ++i) {
- if (this->surface_tractions(i)+surface_energy/rho * exp(- (gap(i)) / rho)<mini )
- mini= this->surface_tractions(i)+surface_energy/rho * exp(- (gap(i)) / rho);
+ if (this->surface_tractions(i)+surface_energy/rho * exp(- (gap(i)) / rho)<mini )
+ mini= this->surface_tractions(i)+surface_energy/rho * exp(- (gap(i)) / rho);
}
-this->surface_tractions-=mini;
+ this->surface_tractions-=mini;
}
+/* -------------------------------------------------------------------------- */
+
__END_TAMAAS__
diff --git a/src/bem_gigipol.hh b/src/bem_gigipol.hh
index 890a46b9..fec1b379 100644
--- a/src/bem_gigipol.hh
+++ b/src/bem_gigipol.hh
@@ -1,110 +1,110 @@
/**
*
* @author Valentine Rey <valentine.rey@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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef BEM_GIGIPOL_H
#define BEM_GIGIPOL_H
/* -------------------------------------------------------------------------- */
#include "bem_fft_base.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
class BemGigipol : public BemFFTBase {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
BemGigipol(Surface<Real> & p) :
BemFFTBase(p),
search_direction(p.size(), p.getL()),
+ spectral_search_direction(p.size(), p.getL()),
+ mean_disp_search_direction(p.size(), p.getL()),
surface_t(p.size(),p.getL()),
surface_r(p.size(),p.getL()),
surface_spectral_r(p.size(),p.getL()),
- spectral_search_direction(p.size(), p.getL()),
- mean_disp_search_direction(p.size(), p.getL()),
- pold(p.size(),p.getL()){
+ pold(p.size(),p.getL()){
e_star = 1.;
max_iterations = 10000;
}
virtual ~BemGigipol(){}
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
Real computeEquilibrium(Real epsilon, Real mean_displacement);
Real computeStoppingCriterion();
Real computeG();
Real computeMeanPressuresInNonContact();
void optimizeToMeanDisplacement(Real imposed_mean);
Real positiveGapAverage(Real shift);
void updateT(Real G, Real Gold, Real delta);
Real updateDisplacements(Real tau, Real mean_displacement);
void enforceMeanDisplacement(Real mean_displacement);
void computePressures(Real mean_displacement);
Real computeTau() ;
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
//! get the number of iterations
UInt getNbIterations() const{return this->nb_iterations;};
//! get the convergence steps of iterations
const std::vector<Real> & getConvergenceIterations() const{return this->convergence_iterations;};
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
//! exploration direction for the Gigi algo
Surface<Real> search_direction;
//! projected exploration direction for the Gigi algo (fourier version)
SurfaceComplex<Real> spectral_search_direction;
//! exploration direction for optimization of mean displacement
Surface<Real> mean_disp_search_direction;
//! exploration direction for the polonski algo
Surface<Real> surface_t;
//! projected exploration direction for the polonski algo
Surface<Real> surface_r;
//! projected exploration direction for the polonski algo (fourier version)
SurfaceComplex<Real> surface_spectral_r;
//! previous pressure
Surface<Real> pold;
};
__END_TAMAAS__
#endif /* BEM_GIGIPOL_H */
diff --git a/src/bem_interface.hh b/src/bem_interface.hh
index 20c5b260..104eb577 100644
--- a/src/bem_interface.hh
+++ b/src/bem_interface.hh
@@ -1,100 +1,102 @@
/**
*
* @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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef BEM_INTERFACE_H
#define BEM_INTERFACE_H
#include "surface.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
class BemInterface {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
BemInterface(Surface<Real> & surface):
surface(surface),
spectral_surface(surface.size(),surface.getL()) {
this->computeSpectralSurface();
};
virtual ~BemInterface(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
protected:
- //! compute the displacement from the pressure field
- virtual void computeDisplacements() = 0;
+ //! compute the displacements from the tractions
+ virtual void computeDisplacementsFromTractions(){SURFACE_FATAL("TO BE IMPLEMENTED")};
+ //! compute the tractions from the displacements
+ virtual void computeTractionsFromDisplacements(){SURFACE_FATAL("TO BE IMPLEMENTED")};
public:
//! compute the spectral surface
virtual void computeSpectralSurface(){
surface.FFTTransform(spectral_surface);
};
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
//! return the actual tractions
virtual const Surface<Real> & getTractions() const = 0;
//! return the actual displacements
virtual const Surface<Real> & getDisplacements() const = 0;
//! return the manipulated surface
const Surface<Real> & getSurface() const{return this->surface;};
//! return the manipulated surface in fourier space
const SurfaceComplex<Real> & getSpectralSurface() const{return this->spectral_surface;};
//! return the manipulated surface
void setSurface(Surface<Real> & surface){this->surface = surface;};
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
//! the surface profile we work with
Surface<Real> & surface;
//! the surface profile in fourier space
SurfaceComplex<Real> spectral_surface;
};
__END_TAMAAS__
#endif /* BEM_INTERFACE_H */
diff --git a/src/bem_kato.cpp b/src/bem_kato.cpp
index f45d90da..51332293 100644
--- a/src/bem_kato.cpp
+++ b/src/bem_kato.cpp
@@ -1,253 +1,253 @@
/**
*
* @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_kato.hh"
#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <cmath>
/* -------------------------------------------------------------------------- */
#define TIMER
#include "surface_timer.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
Real BemKato::linescan(Real scale,Real pressure){
updatePressure(scale);
shiftPressure(pressure);
truncatePressure();
this->functional->computeF();
Real res = this->functional->getF();
return res;
}
/* -------------------------------------------------------------------------- */
Real BemKato::linesearch(Real & hmax, Real fold,
Real pressure, int search_flag){
if (search_flag){
Real h = hmax;
// Real fold = bem.computeF();
//if (fold == 0) fold = 1e300;
Real f = linescan(h,pressure);
if (f < fold) return f;
while (f > fold){
h *= 0.5;
if (h < 1e-3) throw 1;
f = linescan(h,pressure);
}
f = linescan(h,pressure);
// if (hmax / h > 10) hmax /=2;
return f;
}
return linescan(hmax,pressure);
}
/* -------------------------------------------------------------------------- */
Real BemKato::computeEquilibrium(Real epsilon,
Real pressure){
// UInt n = surface.size();
// UInt size = n*n;
- this->computeDisplacements();
+ this->computeDisplacementsFromTractions();
this->functional->computeGradF();
this->backupTractions();
Real f = 1.;
Real fPrevious = 1e300;
this->nb_iterations = 0;
this->convergence_iterations.clear();
while(f > epsilon && this->nb_iterations < this->max_iterations) {
fPrevious = f;
- this->computeDisplacements();
+ this->computeDisplacementsFromTractions();
this->functional->computeGradF();
this->backupTractions();
try {
f = linescan(1.,pressure);
}
catch (int e){
std::cout << " line search failed " << std::endl;
f = linescan(1,pressure);
nb_iterations = 0;
break;
}
if (nb_iterations % dump_freq == 0) {
// std::cout << std::scientific << std::setprecision(10) << nb_iterations << " " << f << " " << f-fold << " " << ((f-fold)/forigin) << std::endl;
std::cout << std::scientific << std::setprecision(10)
<< nb_iterations << " "
<< f << " " << f-fPrevious
<< std::endl;
}
convergence_iterations.push_back(f);
++nb_iterations;
}
this->computeTrueDisplacements();
this->computeGaps();
return f;
}
/* -------------------------------------------------------------------------- */
void BemKato::updatePressure(Real scale){
STARTTIMER("updatePressure");
UInt n = surface.size();
UInt size = n*n;
const Surface<Real> & gradF = this->functional->getGradF();
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
this->surface_tractions(i) =
this->surface_tractions_backup(i) - gradF(i) * scale;
}
STOPTIMER("updatePressure");
}
/* -------------------------------------------------------------------------- */
void BemKato::backupTractions(){
STARTTIMER("switchPressure");
UInt n = surface.size();
UInt size = n*n;
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
this->surface_tractions_backup(i) = this->surface_tractions(i);
}
STOPTIMER("switchPressure");
}
/* -------------------------------------------------------------------------- */
Real BemKato::positivePressure(Real step){
STARTTIMER("positivePressure");
UInt n = surface.size();
UInt size = n*n;
Real p_tot = 0.0;
#pragma omp parallel for reduction(+:p_tot)
for (UInt i = 0; i < size; ++i) {
Real sh_press = this->surface_tractions(i) + step;
if (sh_press > max_pressure) p_tot += max_pressure;
else p_tot += sh_press*(sh_press > 0);
}
STOPTIMER("positivePressure");
return p_tot/n/n;
}
/* -------------------------------------------------------------------------- */
void BemKato::shiftPressure(Real required_pressure){
Real step_min = -10;
Real step_max = 10;
STARTTIMER("shiftPressureInitialGuess");
Real p_max = positivePressure(step_max);
Real p_min = positivePressure(step_min);
for(UInt i = 0; p_max < required_pressure && i < 8; ++i, step_max*=10) {
p_max = positivePressure(step_max);
}
for(UInt i = 0; p_min > required_pressure && i < 8; ++i, step_min*=10) {
p_min = positivePressure(step_min);
}
Real p = positivePressure(0.0);
Real epsilon = 1e-12;
STOPTIMER("shiftPressureInitialGuess");
STARTTIMER("shiftPressureDichotomy");
while (fabs(step_min - step_max) > epsilon){
Real step = (step_min + step_max)/2;
p = positivePressure(step);
if (p > required_pressure) step_max = step;
else if (p < required_pressure) step_min = step;
else {
step_max = step;
step_min = step;
}
}
STOPTIMER("shiftPressureDichotomy");
STARTTIMER("shiftPressureTrueShift");
UInt n = surface.size();
UInt size = n*n;
//shift the pressure so that satisfies the constraint
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
this->surface_tractions(i) += (step_max+step_min)/2;
}
STOPTIMER("shiftPressureTrueShift");
}
/* -------------------------------------------------------------------------- */
void BemKato::truncatePressure(){
STARTTIMER("truncatePressure");
UInt n = surface.size();
UInt size = n*n;
//shift the pressure so that satisfies the constraint
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
if (this->surface_tractions(i) > max_pressure)
this->surface_tractions(i) = max_pressure;
else this->surface_tractions(i) *= (this->surface_tractions(i) > 0);
}
STOPTIMER("truncatePressure");
}
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
diff --git a/src/bem_penalty.cpp b/src/bem_penalty.cpp
index 44e58331..b46bae61 100644
--- a/src/bem_penalty.cpp
+++ b/src/bem_penalty.cpp
@@ -1,241 +1,218 @@
/**
*
* @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_penalty.hh"
#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <cmath>
/* -------------------------------------------------------------------------- */
#define TIMER
#include "surface_timer.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
Real BemPenalty::computeEquilibrium(Real epsilon,
Real mean_displacement, Real penalization_parameter) {
this->computeSpectralInfluenceOverDisplacement();
this->search_direction = 0.;
this->true_displacements =0;
this->true_displacements = SurfaceStatistics::computeMaximum(this->surface);
this->computeGaps();
Real f = 1e300;
std::cout << "moyenne deplacement "<< SurfaceStatistics::computeAverage(this->true_displacements, 0)<< std::endl;
convergence_iterations.clear();
nb_iterations = 0;
max_iterations =5000;
while (f> epsilon && nb_iterations++ < max_iterations) {
this->functional->computeGradF();
this->computeSearchDirection(mean_displacement,penalization_parameter);
Real alpha=0.0001;
this->old_displacements=this->true_displacements;
this->updateUnknown( alpha, mean_displacement);
this->computeGaps();
f = computeStoppingCriterion();
if (nb_iterations % dump_freq == 0) {
std::cout << std::scientific << std::setprecision(10)
<< nb_iterations << " "
<< f << std::fixed << std::endl;
}
convergence_iterations.push_back(f);
}
this->computePressures(mean_displacement);
return f;
}
/* -------------------------------------------------------------------------- */
Real BemPenalty::computeStoppingCriterion() {
Real crit=0.;
Real disp_norm = 0.;
UInt n = surface.size();
UInt size = n*n;
#pragma omp parallel for reduction(+:crit, disp_norm)
for (UInt i = 0; i < size; ++i) {
crit += this->search_direction(i)*this->search_direction(i);
disp_norm += (true_displacements(i)*true_displacements(i));
}
return crit / disp_norm;
}
/* -------------------------------------------------------------------------- */
void BemPenalty::computeSearchDirection(Real mean_displacement,Real penalization_parameter) {
STARTTIMER("computeOptimalStep");
UInt n = surface.size();
UInt size = n*n;
Surface<Real> & gradF = this->functional->getGradF();
#pragma omp parallel for
for (UInt i = 1; i < size; ++i) {
this->search_direction(i)=gradF(i);
if (gap(i)<0)
{
this->search_direction(i)+=penalization_parameter*gap(i);
}
}
}
/* -------------------------------------------------------------------------- */
Real BemPenalty::computeOptimalStep() {
STARTTIMER("computeOptimalStep");
search_direction.FFTTransform(spectral_search_direction, nthreads);
UInt n = surface.size();
UInt size = n*n;
spectral_search_direction(0) = 0;
#pragma omp parallel for
for (UInt i = 1; i < size; ++i) {
spectral_search_direction(i) /= this->surface_spectral_influence_disp(i);
}
// /!\ does not contain the spectral search direction anymore
spectral_search_direction.FFTITransform(nthreads);
// Real average = SurfaceStatistics::computeAverage(spectral_search_direction.real(), 0);
//spectral_search_direction -= average;
Real numerator = 0., denominator = 0.;
#pragma omp parallel for reduction(+: numerator, denominator)
for (UInt i = 0; i < size; ++i) {
numerator += search_direction(i) * search_direction(i);
denominator += spectral_search_direction(i).real() * search_direction(i);
}
Real alpha = numerator / denominator;
STOPTIMER("computeOptimalStep");
return alpha;
}
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
void BemPenalty::updateUnknown(Real alpha, Real mean_displacement) {
STARTTIMER("updateDisplacements");
UInt n = surface.size();
UInt size = n*n;
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
this->true_displacements(i) -= alpha*this->search_direction(i);
}
Real moyenne=SurfaceStatistics::computeAverage(this->true_displacements, 0);
for (UInt i = 0; i < size; ++i) {
this->true_displacements(i) =this->true_displacements(i)-moyenne+mean_displacement;
}
STOPTIMER("updateDisplacements");
}
/* -------------------------------------------------------------------------- */
void BemPenalty::computePressures(Real mean_displacement) {
- UInt n = surface.size();
- UInt size = n*n;
-
- Real moy_surface=SurfaceStatistics::computeAverage(this->surface, 0);
-
- this->true_displacements.FFTTransform(this->surface_spectral_tractions, nthreads);
- this->surface_spectral_tractions(0) = 0;
-
-#pragma omp parallel for
- for (UInt i = 1 ; i < size ; ++i) {
- this->surface_spectral_tractions(i) /= this->surface_spectral_influence_disp(i);
- }
-
- this->surface_spectral_tractions.FFTITransform(this->surface_tractions, nthreads);
-
-
-
-
-
-
-
-
-
-
+ this->computeTractionsFromDisplacements();
}
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
diff --git a/src/bem_polonski.cpp b/src/bem_polonski.cpp
index a5fd79db..26ec1d31 100644
--- a/src/bem_polonski.cpp
+++ b/src/bem_polonski.cpp
@@ -1,314 +1,314 @@
/**
*
* @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->computeDisplacementsFromTractions();
this->computeGaps();
Real gbar = this->computeMeanGapsInContact();
this->gap -= gbar;
Real G = this->computeG();
//std::cout << "G vaut "<< G<< std::endl;
this->updateT(G,Gold,delta);
Real tau = this->computeTau();
pold = this->surface_tractions;
Gold = G;
delta = this->updateTractions(tau);
//Projection on admissible space
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_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/bem_uzawa.cpp b/src/bem_uzawa.cpp
index 093f9ce7..97130b2b 100644
--- a/src/bem_uzawa.cpp
+++ b/src/bem_uzawa.cpp
@@ -1,320 +1,303 @@
/**
*
* @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_uzawa.hh"
#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <cmath>
/* -------------------------------------------------------------------------- */
#define TIMER
#include "surface_timer.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
Real BemUzawa::computeEquilibrium(Real epsilon,
- Real mean_displacement, Real penalization_parameter) {
+ Real mean_displacement, Real penalization_parameter) {
this->computeSpectralInfluenceOverDisplacement();
- this->search_direction = 0.;
- this->true_displacements =0;
- this->true_displacements =mean_displacement;
+ this->search_direction = 0.;
+ this->true_displacements = 0;
+ this->true_displacements = mean_displacement;
this->computeGaps();
- this->multipliers =0;
+ this->multipliers = 0;
- Real f_2 = 1e300;
- Real new_pen = 1e300;
- Real old_pen = 1e300;
+ Real f_2 = 1e300;
+ Real new_pen = 1e300;
+ Real old_pen = 1e300;
- convergence_iterations.clear();
- UInt nb_iterations_1 = 0;
- UInt nb_iterations_2 = 0;
+ convergence_iterations.clear();
+ UInt nb_iterations_1 = 0;
+ UInt nb_iterations_2 = 0;
- Real max_iterations_1 =500;
- Real max_iterations_2 =1000;
+ Real max_iterations_1 = 500;
+ Real max_iterations_2 = 1000;
while (new_pen> epsilon && nb_iterations_1++ < max_iterations_1) {
- std::cout << "iter ext "<< nb_iterations_1 << std::endl;
+ std::cout << "iter ext "<< nb_iterations_1 << std::endl;
while (f_2> 1e-12 && nb_iterations_2++ < max_iterations_2) {
- std::cout << "iter int "<< nb_iterations_2 << std::endl;
- this->functional->computeGradF();
- this->computeSearchDirection(mean_displacement,penalization_parameter);
+ std::cout << "iter int "<< nb_iterations_2 << std::endl;
+ this->functional->computeGradF();
+ this->computeSearchDirection(mean_displacement,penalization_parameter);
- Real alpha=1/(10*penalization_parameter);
- this->old_displacements=this->true_displacements;
- this->updateUnknown( alpha, mean_displacement);
- this->computeGaps();
+ Real alpha=1/(10*penalization_parameter);
+ this->old_displacements=this->true_displacements;
+ this->updateUnknown( alpha, mean_displacement);
+ this->computeGaps();
- f_2 = computeStoppingCriterion();
+ f_2 = computeStoppingCriterion();
- if (nb_iterations_2 % dump_freq == 0) {
- std::cout << std::scientific << std::setprecision(10)
- << nb_iterations << " "
- << f_2 << std::fixed << std::endl;
+ if (nb_iterations_2 % dump_freq == 0) {
+ std::cout << std::scientific << std::setprecision(10)
+ << nb_iterations << " "
+ << f_2 << std::fixed << std::endl;
+ }
}
- }
this->updateMultipliers(penalization_parameter);
old_pen=computeInterpenetration(this->old_displacements);
new_pen=computeInterpenetration(this->true_displacements);
std::cout << "old penetration is "<< old_pen << std::endl;
std::cout << "new penetration is "<< new_pen << std::endl;
penalization_parameter=this->updatePenalization(penalization_parameter, old_pen, new_pen);
//to avoid ill conditioned system
if (penalization_parameter>1.0e8) new_pen=0.;
nb_iterations_2 = 0;
- f_2 = 1e300;
+ f_2 = 1e300;
std::cout << "penalization is "<< penalization_parameter << std::endl;
- }
+ }
this->computeGaps();
- this->computePressures(mean_displacement);
+ this->computePressures();
return new_pen;
}
/* -------------------------------------------------------------------------- */
Real BemUzawa::computeStoppingCriterion() {
- Real crit=0.;
- Real disp_norm = 0.;
+ Real crit=0.;
+ Real disp_norm = 0.;
- UInt n = surface.size();
- UInt size = n*n;
+ UInt n = surface.size();
+ UInt size = n*n;
- #pragma omp parallel for reduction(+:crit, disp_norm)
- for (UInt i = 0; i < size; ++i) {
+#pragma omp parallel for reduction(+:crit, disp_norm)
+ for (UInt i = 0; i < size; ++i) {
- crit += (this->search_direction(i))*(this->search_direction(i));
- disp_norm += (true_displacements(i)*true_displacements(i));
- }
+ crit += (this->search_direction(i))*(this->search_direction(i));
+ disp_norm += (true_displacements(i)*true_displacements(i));
+ }
- return crit / disp_norm;
+ return crit / disp_norm;
}
/* -------------------------------------------------------------------------- */
void BemUzawa::computeSearchDirection(Real mean_displacement,Real penalization_parameter) {
STARTTIMER("computeOptimalStep");
UInt n = surface.size();
UInt size = n*n;
- Surface<Real> & gradF = this->functional->getGradF();
+ Surface<Real> & gradF = this->functional->getGradF();
#pragma omp parallel for
for (UInt i = 1; i < size; ++i) {
- this->search_direction(i)=gradF(i);
+ this->search_direction(i)=gradF(i);
- if (this->multipliers(i)+penalization_parameter*gap(i)<=0){
+ if (this->multipliers(i)+penalization_parameter*gap(i)<=0){
- this->search_direction(i)=this->search_direction(i)+ this->multipliers(i)+penalization_parameter*gap(i);
+ this->search_direction(i)=this->search_direction(i)+ this->multipliers(i)+penalization_parameter*gap(i);
}
-}
+ }
}
/* -------------------------------------------------------------------------- */
Real BemUzawa::computeOptimalStep() {
STARTTIMER("computeOptimalStep");
search_direction.FFTTransform(spectral_search_direction, nthreads);
UInt n = surface.size();
UInt size = n*n;
spectral_search_direction(0) = 0;
#pragma omp parallel for
- for (UInt i = 1; i < size; ++i) {
+ for (UInt i = 1; i < size; ++i) {
spectral_search_direction(i) /= this->surface_spectral_influence_disp(i);
}
- // /!\ does not contain the spectral search direction anymore
+ // /!\ does not contain the spectral search direction anymore
spectral_search_direction.FFTITransform(nthreads);
- // Real average = SurfaceStatistics::computeAverage(spectral_search_direction.real(), 0);
-//spectral_search_direction -= average;
+ // Real average = SurfaceStatistics::computeAverage(spectral_search_direction.real(), 0);
+ //spectral_search_direction -= average;
- Real numerator = 0., denominator = 0.;
+ Real numerator = 0., denominator = 0.;
#pragma omp parallel for reduction(+: numerator, denominator)
- for (UInt i = 0; i < size; ++i) {
- numerator += search_direction(i) * search_direction(i);
+ for (UInt i = 0; i < size; ++i) {
+ numerator += search_direction(i) * search_direction(i);
denominator += spectral_search_direction(i).real() * search_direction(i);
- }
+ }
- Real alpha = numerator / denominator;
+ Real alpha = numerator / denominator;
STOPTIMER("computeOptimalStep");
- return alpha;
+ return alpha;
}
/* -------------------------------------------------------------------------- */
void BemUzawa::updateMultipliers(Real penalization_parameter) {
STARTTIMER("updateMultipliers");
UInt n = surface.size();
UInt size = n*n;
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
- if (this->multipliers(i)+penalization_parameter*gap(i)>0.){
- this->multipliers(i) = 0;
- }
- else{
- this->multipliers(i)+=penalization_parameter*gap(i);
- }
+ if (this->multipliers(i)+penalization_parameter*gap(i)>0.){
+ this->multipliers(i) = 0;
+ }
+ else{
+ this->multipliers(i)+=penalization_parameter*gap(i);
+ }
}
STOPTIMER("updateMultipliers");
}
/* -------------------------------------------------------------------------- */
Real BemUzawa::updatePenalization(Real penalization_parameter, Real old_pen, Real new_pen) {
STARTTIMER("updatePenalization");
- if (new_pen>0.25*old_pen){
- penalization_parameter*=5;
+ if (new_pen>0.25*old_pen){
+ penalization_parameter*=5;
}
- return penalization_parameter;
+ return penalization_parameter;
STOPTIMER("updatePenalization");
}
/* -------------------------------------------------------------------------- */
Real BemUzawa::computeInterpenetration( Surface<Real> & displacements ) {
STARTTIMER("updateMultipliers");
UInt n = surface.size();
UInt size = n*n;
Real res=0.;
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
- if (0>this->true_displacements(i)-surface(i)){
- res+=(this->true_displacements(i)-surface(i))*(this->true_displacements(i)-surface(i));
- }
+ if (0>this->true_displacements(i)-surface(i)){
+ res+=(this->true_displacements(i)-surface(i))*(this->true_displacements(i)-surface(i));
+ }
-}
+ }
- return sqrt(res);
+ return sqrt(res);
STOPTIMER("updateMultipliers");
}
/* -------------------------------------------------------------------------- */
void BemUzawa::updateUnknown(Real alpha, Real mean_displacement) {
STARTTIMER("updateDisplacements");
UInt n = surface.size();
UInt size = n*n;
#pragma omp parallel for
for (UInt i = 0; i < size; ++i) {
this->true_displacements(i) -= alpha*this->search_direction(i);
}
-Real moyenne=SurfaceStatistics::computeAverage(this->true_displacements, 0);
+ Real moyenne=SurfaceStatistics::computeAverage(this->true_displacements, 0);
-for (UInt i = 0; i < size; ++i) {
- this->true_displacements(i) =this->true_displacements(i)-moyenne+mean_displacement;
+ for (UInt i = 0; i < size; ++i) {
+ this->true_displacements(i) =this->true_displacements(i)-moyenne+mean_displacement;
-}
+ }
STOPTIMER("updateDisplacements");
}
/* -------------------------------------------------------------------------- */
-void BemUzawa::computePressures(Real mean_displacement) {
- UInt n = surface.size();
- UInt size = n*n;
-
- Real moy_surface=SurfaceStatistics::computeAverage(this->surface, 0);
-
- this->true_displacements.FFTTransform(this->surface_spectral_tractions, nthreads);
- this->surface_spectral_tractions(0) = 0;
-
-#pragma omp parallel for
- for (UInt i = 1 ; i < size ; ++i) {
- this->surface_spectral_tractions(i) /= this->surface_spectral_influence_disp(i);
- }
-
- this->surface_spectral_tractions.FFTITransform(this->surface_tractions, nthreads);
-
-
-
-
+void BemUzawa::computePressures() {
+ this->applyInverseInfluenceFunctions(this->true_displacements, this->surface_tractions);
}
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
diff --git a/src/bem_uzawa.hh b/src/bem_uzawa.hh
index 474a720c..2ec22d23 100644
--- a/src/bem_uzawa.hh
+++ b/src/bem_uzawa.hh
@@ -1,102 +1,100 @@
/**
*
* @author Valentine Rey <valentine.rey@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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef BEM_UZAWA_H
#define BEM_UZAWA_H
/* -------------------------------------------------------------------------- */
#include "bem_fft_base.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
class BemUzawa : public BemFFTBase {
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
BemUzawa(Surface<Real> & p) :
BemFFTBase(p),
search_direction(p.size(), p.getL()),
spectral_search_direction(p.size(), p.getL()),
multipliers(p.size(),p.getL()),
surface_r(p.size(),p.getL()),
surface_spectral_r(p.size(),p.getL()),
mean_disp_search_direction(p.size(), p.getL()) {
e_star = 1.;
max_iterations = 10000;
};
virtual ~BemUzawa(){};
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
//! compute the equilibrium situation
Real computeEquilibrium(Real epsilon,Real pressure,Real penalization_parameter);
void computeSearchDirection(Real mean_displacement,Real penalization_parameter);
- void updateMultipliers(Real penalization_parameter);
- Real updatePenalization(Real penalization_parameter, Real old_pen, Real new_pen);
- Real computeInterpenetration( Surface<Real> & displacements ) ;
+ void updateMultipliers(Real penalization_parameter);
+ Real updatePenalization(Real penalization_parameter, Real old_pen, Real new_pen);
+ Real computeInterpenetration( Surface<Real> & displacements ) ;
void updateUnknown(Real tau, Real mean_displacement);
Real computeStoppingCriterion();
- void computePressures(Real mean_displacement);
- Real computeOptimalStep();
-
-
+ void computePressures();
+ Real computeOptimalStep();
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
//! get the number of iterations
UInt getNbIterations() const{return this->nb_iterations;};
//! get the convergence steps of iterations
const std::vector<Real> & getConvergenceIterations() const{return this->convergence_iterations;};
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
//! exploration direction for the Gigi algo
Surface<Real> search_direction;
//! projected exploration direction for the Gigi algo (fourier version)
SurfaceComplex<Real> spectral_search_direction;
//! exploration direction for the polonski algo
Surface<Real> multipliers;
//! projected exploration direction for the polonski algo
Surface<Real> surface_r;
//! projected exploration direction for the polonski algo (fourier version)
SurfaceComplex<Real> surface_spectral_r;
//! exploration direction for optimization of mean displacement
Surface<Real> mean_disp_search_direction;
};
__END_TAMAAS__
#endif /* BEM_Uzawa_H */
diff --git a/src/fft_plan_manager.cpp b/src/fft_plan_manager.cpp
new file mode 100644
index 00000000..2e15d15a
--- /dev/null
+++ b/src/fft_plan_manager.cpp
@@ -0,0 +1,89 @@
+/**
+ *
+ * @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 "fft_plan_manager.hh"
+/* -------------------------------------------------------------------------- */
+
+
+__BEGIN_TAMAAS__
+/* -------------------------------------------------------------------------- */
+
+FFTPlanManager::FFTPlanManager(){
+}
+/* -------------------------------------------------------------------------- */
+
+FFTPlanManager::~FFTPlanManager(){
+}
+/* -------------------------------------------------------------------------- */
+
+
+FFTPlanManager & FFTPlanManager::get(){
+ if (FFTPlanManager::singleton == NULL)
+ FFTPlanManager::singleton = new FFTPlanManager();
+
+ return *FFTPlanManager::singleton;
+}
+/* -------------------------------------------------------------------------- */
+
+FFTransform<Real> & FFTPlanManager::createPlan(Surface<Real> & input,
+ SurfaceComplex<Real> & output){
+
+ auto index = std::make_pair(const_cast<Real*>(input.getInternalData()),
+ const_cast<complex*>(output.getInternalData()));
+
+ auto it = plans.find(index);
+ auto end = plans.end();
+
+ if (it == end){
+ plans[index] = new FFTransform<Real>(input.size(),input,output);
+ }
+
+ return *plans[index];
+}
+
+/* -------------------------------------------------------------------------- */
+
+void FFTPlanManager::destroyPlan(Surface<Real> & input,
+ SurfaceComplex<Real> & output){
+
+ auto index = std::make_pair(const_cast<Real*>(input.getInternalData()),
+ const_cast<complex*>(output.getInternalData()));
+
+ auto it = plans.find(index);
+ auto end = plans.end();
+
+ if (it != end){
+ delete (plans[index]);
+ plans.erase(index);
+ }
+}
+
+/* -------------------------------------------------------------------------- */
+FFTPlanManager * FFTPlanManager::singleton = NULL;
+/* -------------------------------------------------------------------------- */
+
+
+
+__END_TAMAAS__
diff --git a/src/bem_interface.hh b/src/fft_plan_manager.hh
similarity index 56%
copy from src/bem_interface.hh
copy to src/fft_plan_manager.hh
index 20c5b260..90d5bde7 100644
--- a/src/bem_interface.hh
+++ b/src/fft_plan_manager.hh
@@ -1,100 +1,77 @@
/**
*
* @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/>.
*
*/
-
/* -------------------------------------------------------------------------- */
-#ifndef BEM_INTERFACE_H
-#define BEM_INTERFACE_H
-
-#include "surface.hh"
+#ifndef FFT_PLAN_MANAGER_H
+#define FFT_PLAN_MANAGER_H
+/* -------------------------------------------------------------------------- */
+#include <map>
+#include "fftransform.hh"
/* -------------------------------------------------------------------------- */
-
__BEGIN_TAMAAS__
+/* -------------------------------------------------------------------------- */
+class FFTPlanManager {
-
-class BemInterface {
-
+
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
-public:
-
- BemInterface(Surface<Real> & surface):
- surface(surface),
- spectral_surface(surface.size(),surface.getL()) {
- this->computeSpectralSurface();
- };
- virtual ~BemInterface(){};
-
- /* ------------------------------------------------------------------------ */
- /* Methods */
- /* ------------------------------------------------------------------------ */
-
protected:
- //! compute the displacement from the pressure field
- virtual void computeDisplacements() = 0;
-public:
+ FFTPlanManager();
- //! compute the spectral surface
- virtual void computeSpectralSurface(){
- surface.FFTTransform(spectral_surface);
- };
+public:
+ virtual ~FFTPlanManager();
+public:
+
/* ------------------------------------------------------------------------ */
- /* Accessors */
+ /* Methods */
/* ------------------------------------------------------------------------ */
-public:
-
- //! return the actual tractions
- virtual const Surface<Real> & getTractions() const = 0;
- //! return the actual displacements
- virtual const Surface<Real> & getDisplacements() const = 0;
- //! return the manipulated surface
- const Surface<Real> & getSurface() const{return this->surface;};
- //! return the manipulated surface in fourier space
- const SurfaceComplex<Real> & getSpectralSurface() const{return this->spectral_surface;};
- //! return the manipulated surface
- void setSurface(Surface<Real> & surface){this->surface = surface;};
-
+ static FFTPlanManager & get();
+ FFTransform<Real> & createPlan(Surface<Real> & input, SurfaceComplex<Real> & output);
+ void destroyPlan(Surface<Real> & input, SurfaceComplex<Real> & output);
+
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
-protected:
+private:
+
+ std::map<std::pair<Real*,complex*>,FFTransform<Real> *> plans;
- //! the surface profile we work with
- Surface<Real> & surface;
- //! the surface profile in fourier space
- SurfaceComplex<Real> spectral_surface;
+ static FFTPlanManager * singleton;
};
+/* -------------------------------------------------------------------------- */
__END_TAMAAS__
+/* -------------------------------------------------------------------------- */
+
-#endif /* BEM_INTERFACE_H */
+#endif /* FFT_PLAN_MANAGER_H */
diff --git a/src/map_2d.cpp b/src/map_2d.cpp
index 9ff44ae9..2650acb8 100644
--- a/src/map_2d.cpp
+++ b/src/map_2d.cpp
@@ -1,246 +1,246 @@
/**
*
* @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 "map_2d.hh"
#include <fstream>
#include <sstream>
#include <fftw3.h>
#include <omp.h>
/* -------------------------------------------------------------------------- */
#define TIMER
#include "surface_timer.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
template <typename T>
void Map2d<T>::rescaleHeights(Real C){
/* compute actual RMS */
for (UInt i = 0 ; i < n[0] ; ++i)
for (UInt j = 0 ; j < n[1] ; ++j){
this->at(i,j) *= C;
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::setGridSize(UInt n0, UInt n1){
n[0] = n0;
n[1] = n1;
data.resize(n[0]*n[1]);
data.assign(n[0]*n[1],T(0.));
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::setGridSize(const UInt n[2]){
setGridSize(n[0],n[1]);
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::operator=(const T & e){
for (UInt i = 0 ; i < n[0]*n[1] ; ++i){
this->at(i) = e;
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
Map2d<T>::Map2d(UInt n0, UInt n1, Real L0, Real L1) {
setGridSize(n0,n1);
this->L[0] = L0;
this->L[1] = L1;
this->plan = NULL;
generateFFTPlan();
}
/* -------------------------------------------------------------------------- */
template <typename T>
Map2d<T>::Map2d(const UInt * n, const Real *L) {
setGridSize(n[0],n[1]);
this->L[0] = L[0];
this->L[1] = L[1];
this->plan = NULL;
generateFFTPlan();
}
/* -------------------------------------------------------------------------- */
template <typename T>
Map2d<T>::~Map2d() {
destroyFFTPlan();
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::generateFFTPlan() {
// this should not do anything
}
/* -------------------------------------------------------------------------- */
template <>
void Map2d<complex>::generateFFTPlan() {
fftw_complex * in = (fftw_complex*) &data[0];
this->plan = fftw_plan_dft_2d(n[0], n[1], in, in, FFTW_FORWARD, FFTW_ESTIMATE);
}
/* -------------------------------------------------------------------------- */
template<typename T>
void Map2d<T>::destroyFFTPlan() {
if (this->plan) {
fftw_destroy_plan(this->plan);
this->plan = NULL;
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::FFTTransform(int nthreads){
SURFACE_FATAL("impossible situation");
throw;
}
/* -------------------------------------------------------------------------- */
template <typename T>
-void Map2d<T>::FFTTransform(Map2d<std::complex<Real> > & output, int nthreads){
+void Map2d<T>::FFTTransform(Map2d<std::complex<Real> > & output, int nthreads)const{
output = *this;
output.FFTTransform(nthreads);
}
/* -------------------------------------------------------------------------- */
template <typename T>
-void Map2d<T>::FFTITransform(Map2d<std::complex<Real> > & output, int nthreads){
+void Map2d<T>::FFTITransform(Map2d<std::complex<Real> > & output, int nthreads)const{
output = *this;
output.FFTITransform(nthreads);
}
/* -------------------------------------------------------------------------- */
template <typename T>
-void Map2d<T>::FFTITransform(Map2d<Real> & output, int nthreads){
+void Map2d<T>::FFTITransform(Map2d<Real> & output, int nthreads) const{
const UInt * n = this->sizes();
const Real * Ls = this->getLs();
Map2d<std::complex<Real> > _output(n,Ls);
_output = *this;
_output.FFTITransform(nthreads);
for (UInt i = 0; i < n[0]*n[1]; i++) {
output(i) = _output(i).real();
}
}
/* -------------------------------------------------------------------------- */
template <>
void Map2d<complex>::FFTTransform(int nthreads){
STARTTIMER("FFTW");
fftw_execute(this->plan);
STOPTIMER("FFTW");
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::FFTITransform(int nthreads){
SURFACE_FATAL("impossible situation");
throw;
}
/* -------------------------------------------------------------------------- */
template <>
void Map2d<complex>::FFTITransform(int nthreads){
#pragma omp parallel for
for (UInt i = 0; i < n[0]*n[1]; ++i) {
data[i] = std::conj(data[i]);
}
FFTTransform(nthreads);
#pragma omp parallel for
for (UInt i = 0; i < n[0]*n[1]; ++i) {
data[i] = std::conj(data[i]);
data[i] *= 1./(n[0]*n[1]);
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::dumpToTextFile(std::string filename) {
UInt n1 = this->size(0);
UInt n2 = this->size(1);
std::cout << "Writing surface file " << filename << std::endl;
std::ofstream file(filename.c_str());
if (!file.is_open()) SURFACE_FATAL("cannot open file " << filename);
file.precision(15);
for (UInt i = 0; i < n1; i++) {
for (UInt j=0; j< n2; j++) {
file << i << " " << j << " ";
this->writeValue(file,this->at(i,j));
file << std::endl;
}
file << std::endl;
}
file.close();
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::writeValue(std::ostream & os, const T & val) {
os << val;
}
/* -------------------------------------------------------------------------- */
template <>
void Map2d<complex>::writeValue(std::ostream & os, const complex & val) {
os << std::scientific << val.real();
os << " " << val.imag();
}
/* -------------------------------------------------------------------------- */
template class Map2d<complex>;
template class Map2d<Real>;
template class Map2d<int>;
template class Map2d<unsigned int>;
template class Map2d<bool>;
template class Map2d<unsigned long>;
__END_TAMAAS__
diff --git a/src/map_2d.hh b/src/map_2d.hh
index a7487eca..6dd42708 100644
--- a/src/map_2d.hh
+++ b/src/map_2d.hh
@@ -1,482 +1,591 @@
/**
*
* @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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __MAP_2D_HH__
#define __MAP_2D_HH__
/* -------------------------------------------------------------------------- */
#include <vector>
#include <iostream>
#include <fftw3.h>
#include "tamaas.hh"
#include <fstream>
#include "my_vector.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
template <typename T>
class Map2d {
/* ------------------------------------------------------------------------ */
/* Typedefs */
/* ------------------------------------------------------------------------ */
// class iterator;
/* ------------------------------------------------------------------------ */
/* Constructors/Destructors */
/* ------------------------------------------------------------------------ */
public:
Map2d(UInt n1, UInt n2, Real L0 = 1., Real L1 = 1.);
Map2d(const UInt * n, const Real * L);
+
virtual ~Map2d();
/* ------------------------------------------------------------------------ */
/* Methods */
/* ------------------------------------------------------------------------ */
public:
template <bool periodic>
std::vector<std::pair<UInt,UInt> > getNeighborIndexes(UInt i, UInt j);
template <bool periodic>
std::vector<std::pair<UInt,UInt> > getNeighborIndexes(std::pair<UInt,UInt> i);
// //! return an iterator over the neighbors
// iterator beginNeighborhood();
// //! return the end iterator over the neighbors
// iterator endNeighborhood();
//! generate fft plan
void generateFFTPlan();
//! destroy fft plan
void destroyFFTPlan();
//! perform fft transform
- void FFTTransform(Map2d<std::complex<Real> > & output,int nthreads = 1);
+ void FFTTransform(Map2d<std::complex<Real> > & output,int nthreads = 1) const;
//! perform inverse fft transform
- void FFTITransform(Map2d<std::complex<Real> > & output, int nthreads = 1);
+ void FFTITransform(Map2d<std::complex<Real> > & output, int nthreads = 1) const;
//! perform inverse fft transform
- void FFTITransform(Map2d<Real> & output, int nthreads = 1);
+ void FFTITransform(Map2d<Real> & output, int nthreads = 1) const;
//! perform fft transform
void FFTTransform(int nthreads = 1);
// //! inverse fft transform
void FFTITransform(int nthreads = 1);
//compute compute component by component product and put result in object
template <typename T1, typename T2>
void fullProduct(Map2d<T1> & s1,Map2d<T2> &s2);
//! set the grid size
void setGridSize(const UInt n[2]);
//! set the grid size
void setGridSize(UInt n0,UInt n1);
//! rescale height profile
void rescaleHeights(Real factor);
//! Write heights array to text file usable with plot programs
void dumpToTextFile(std::string filename);
+ //! copy method
+ template <template<typename> class surf, typename T2>
+ void copy(const surf<T2> & s);
// set all the members of the map to a given value
void operator=(const T & e);
-
+ // multiply the map by scalar
+ template <typename T2>
+ void operator*=(const T2 & e);
+ // divide the map by scalar
+ template <typename T2>
+ void operator/=(const T2 & e);
+ // multiply the map by scalar
+ void operator*=(const T & e);
+ // divide the map by scalar
+ void operator/=(const T & e);
+ // multiply the map with another map, component by component
+ template <template<typename> class surf, typename T2>
+ void operator*=(const surf<T2> & e);
+ // divide the map with another map, component by component
+ template <template<typename> class surf, typename T2>
+ void operator/=(const surf<T2> & e);
+ // copy the map passed as parameter
template <template<typename> class surf, typename T2>
void operator=(const surf<T2> & s);
+ // temporary hack
+ void operator=(const Map2d<complex> & s);
+ // copy the map passed as const parameter
template <template<typename> class surf, typename T2>
void operator=(surf<T2> & s);
/* -------------------------------------------------------------------------- */
private:
//! write an entry to a stringstream
void writeValue(std::ostream & os, const T & val);
/* ------------------------------------------------------------------------ */
/* Accessors */
/* ------------------------------------------------------------------------ */
public:
// access to the i,j-th value in the surface
inline void operator+=(const Map2d<T> & s);
// access to the i,j-th value in the surface
inline void operator+=(const T & val);
// access to the i,j-th value in the surface
inline void operator-=(const T & val);
// access to the i,j-th value in the surface
inline T & operator()(UInt i,UInt j);
// access to the i,j-th value in the surface
inline T & operator()(std::pair<UInt,UInt> i);
// access to the i,j-th value in the surface
inline const T & operator()(UInt i,UInt j)const;
// access to the i,j-th value in the surface
inline const T & operator()(std::pair<UInt,UInt> i)const;
// access to the i,j-th value in the surface
inline T & at(UInt i,UInt j);
// access to the i-th value in the array
inline T & at(UInt i);
// access to the i,j-th value in the surface
inline T & at(std::pair<UInt,UInt> i);
// access to the i,j-th value in the surface
inline const T & at(UInt i,UInt j) const;
// access to the i,j-th value in the surface
inline const T & at(std::pair<UInt,UInt> i) const;
// access to the i,j-th value in the surface
inline const T & at(UInt i) const;
// access to the i-th value in the array
inline T & operator()(int i);
// access to the i-th value in the array
// inline T & operator[](std::pair<UInt,UInt> i);
// access to the i-th value in the array
inline const T & operator()(int i) const;
//! get lateral number of discretization points
inline UInt size(UInt dir)const {return n[dir];};
//! get both lateral number of discretization points
inline const UInt* sizes()const {return n;};
//! get the lateral size of the sample
inline const Real * getLs() const {return L;};
//! get the lateral size of the sample
inline Real getL(UInt dir) const {return L[dir];};
//! get the lateral size of the sample
inline void setL(Real L[2]){this->L[0] = L[0];this->L[1] = L[1];};
//! get the lateral size of the sample
inline void setL(Real L0,Real L1){this->L[0] = L0;this->L[1] = L1;};
const T * getInternalData() const{return data.getPtr();};
//! Write heights array to text file usable with plot programs
template <bool consider_zeros=true>
void dumpToTextFile(std::string filename) const;
/* ------------------------------------------------------------------------ */
/* Class Members */
/* ------------------------------------------------------------------------ */
protected:
//! number of points in each direction
UInt n[2];
//! size in each direction
Real L[2];
//! internal storage of data
MyVector<T> data;
//! fftw plan
fftw_plan plan;
};
/* -------------------------------------------------------------------------- */
//#define _at(I,J,N) I+J*N[1]
#define _at(I,J,N) I*N[1]+J
/* -------------------------------------------------------------------------- */
template <typename T>
inline T & Map2d<T>::at(UInt i,UInt j){
if (_at(i,j,n) >= n[0]*n[1]) SURFACE_FATAL("out of bound access "
<< i << " " << j << " "
<< n[0] << " " << n[1] << " "
<< _at(i,j,n));
return data[_at(i,j,n)];
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline const T & Map2d<T>::at(UInt i,UInt j)const {
if (_at(i,j,n) >= n[0]*n[1]) SURFACE_FATAL("out of bound access "
<< i << " " << j << " "
<< n[0] << " " << n[1] << " "
<< _at(i,j,n));
return data[_at(i,j,n)];
}
/* -------------------------------------------------------------------------- */
#undef _at
/* -------------------------------------------------------------------------- */
template <typename T>
inline T & Map2d<T>::at(UInt i){
return data[i];
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline const T & Map2d<T>::at(UInt i) const{
return data[i];
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline T & Map2d<T>::operator()(UInt i,UInt j){
return this->at(i,j);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline T & Map2d<T>::at(std::pair<UInt,UInt> i){
return this->at(i.first,i.second);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline const T & Map2d<T>::at(std::pair<UInt,UInt> i)const {
return this->at(i.first,i.second);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline T & Map2d<T>::operator()(std::pair<UInt,UInt> i){
return this->at(i.first,i.second);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline const T & Map2d<T>::operator()(UInt i,UInt j)const{
return this->at(i,j);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline const T & Map2d<T>::operator()(std::pair<UInt,UInt> i)const {
return this->at(i.first,i.second);
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline T & Map2d<T>::operator()(int i){
return data[i];
}
/* -------------------------------------------------------------------------- */
template <typename T>
inline const T & Map2d<T>::operator()(int i) const{
return data[i];
}
/* -------------------------------------------------------------------------- */
// template <typename T>
// inline T & Map2d<T>::operator[](std::pair<UInt,UInt> c){
// return this->at(c.first,c.second);
// }
// /* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::operator-=(const T & e){
for (UInt i = 0 ; i < n[0]*n[1] ; ++i){
this->at(i) -= e;
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::operator+=(const T & e){
for (UInt i = 0 ; i < n[0]*n[1] ; ++i){
this->at(i) += e;
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
void Map2d<T>::operator+=(const Map2d<T> & s){
for (UInt i = 0 ; i < n[0]*n[1] ; ++i){
this->at(i) += s.at(i);
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
template <bool periodic>
std::vector<std::pair<UInt,UInt> > Map2d<T>::getNeighborIndexes(UInt i, UInt j){
std::vector<std::pair<UInt,UInt> > vres;
std::pair<int,int> res;
for (int ii = -1; ii < 2; ++ii) {
res.first = ii+static_cast<int>(i);
if (res.first < 0 && !periodic) continue;
if (res.first >= (int)n[0] && !periodic) continue;
if (res.first < 0 && periodic) res.first += n[0];
if (res.first >= (int)n[0] && periodic) res.first -= n[0];
for (int jj = -1; jj < 2; ++jj) {
res.second = jj+static_cast<int>(j);
if (res.second < 0 && !periodic) continue;
if (res.second >= (int)n[1] && !periodic) continue;
if (res.second < 0 && periodic) res.second += n[1];
if (res.second >= (int)n[1] && periodic) res.second -= n[1];
if (ii == 0 && jj == 0) continue;
{
UInt i = res.first;
UInt j = res.second;
if (i+j*n[0] >= n[0]*n[1]) SURFACE_FATAL("out of bound access "
<< i << " " << j << " "
<< n[0] << " " << n[1] << " "
<< i+j*n[0]);
}
vres.push_back(res);
}
}
return vres;
}
/* -------------------------------------------------------------------------- */
template <typename T>
template <bool periodic>
std::vector<std::pair<UInt,UInt> > Map2d<T>
::getNeighborIndexes(std::pair<UInt,UInt> i){
return this->getNeighborIndexes<periodic>(i.first,i.second);
}
/* -------------------------------------------------------------------------- */
template <typename T>
template <typename T1, typename T2>
void Map2d<T>::fullProduct(Map2d<T1> & s1,Map2d<T2> & s2){
auto n1 = s1.sizes();
auto n2 = s2.sizes();
if (n1[0] != n2[0] ||
n1[1] != n2[1]) SURFACE_FATAL("fullproduct of non compatible matrices");
if (n[0] != n1[0] ||
n[1] != n1[1]) SURFACE_FATAL("fullproduct of non compatible matrices")
for (UInt i = 0 ; i < n[0] ; ++i)
for (UInt j = 0 ; j < n[1] ; ++j){
this->at(i,j) = s1(i,j)* s2(i,j);
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
template <bool consider_zeros>
void Map2d<T>::dumpToTextFile(std::string filename) const {
UInt m = this->n[0];
UInt n = this->n[1];
std::cout << "Writing surface file " << filename << std::endl;
std::ofstream file(filename.c_str());
if (!file.is_open()) SURFACE_FATAL("cannot open file " << filename);
file.precision(15);
for (UInt i = 0; i < m; i++) {
for (UInt j=0; j< n; j++) {
if (!consider_zeros &&
std::norm(this->at(i,j)) < zero_threshold*zero_threshold)
continue;
file << i << " " << j << " ";
file << std::scientific << this->at(i,j);
file << std::endl;
}
file << std::endl;
}
file.close();
}
// template <typename T>
// class Map2d<T>::iterator {
// public:
// iterator(std::vector<T> & data,
// std::queue<std::pair<UInt,UInt> > & indices):
// data(data),indices(indices){
// }
// inline void operator++(){
// indices.pop();
// }
// inline T & operator *(){
// return data[indices.front];
// }
// bool operator != (Map2d<T>::iterator & it){
// return (indices.front == it.indices.front);
// }
// private:
// std::vector<T> & data;
// std::queue<std::pair<UInt,UInt> > indices;
// bool finished;
// };
/* -------------------------------------------------------------------------- */
-
template <typename T>
template <template<typename> class surf, typename T2>
-void Map2d<T>::operator=(const surf<T2> & s){
+void Map2d<T>::copy(const surf<T2> & s){
this->setGridSize(s.sizes());
UInt m = this->sizes()[0];
UInt n = this->sizes()[1];
#pragma omp parallel for
for (UInt i = 0; i < m*n; i++) {
this->at(i) = s(i);
}
}
/* -------------------------------------------------------------------------- */
template <typename T>
template <template<typename> class surf, typename T2>
void Map2d<T>::operator=(surf<T2> & s){
- this->setGridSize(s.sizes());
- UInt m = this->sizes()[0];
- UInt n = this->sizes()[1];
+ this->copy(s);
+}
+/* -------------------------------------------------------------------------- */
+template <typename T>
+template <template<typename> class surf, typename T2>
+void Map2d<T>::operator=(const surf<T2> & s){
+ this->copy(s);
+}
+
+/* -------------------------------------------------------------------------- */
+
+template <typename T>
+void Map2d<T>::operator=(const Map2d<complex> & s){
+ SURFACE_FATAL("Cannot copy a complex surface into incompatible type");
+}
+template <>
+inline void Map2d<complex>::operator=(const Map2d<complex> & s){
+ this->copy(s);
+}
+
+/* -------------------------------------------------------------------------- */
+template <typename T>
+template <typename T2>
+void Map2d<T>::operator*=(const T2 & e){
+
+ UInt sz = this->n[0]*this->n[1];
#pragma omp parallel for
- for (UInt i = 0; i < m*n; i++) {
- this->at(i) = s(i);
+ for (unsigned int i = 0; i < sz; ++i) {
+ this->at(i) *= e;
+ }
+
+}
+/* -------------------------------------------------------------------------- */
+
+template <typename T>
+template <typename T2>
+void Map2d<T>::operator/=(const T2 & e){
+
+ UInt sz = this->n[0]*this->n[1];
+#pragma omp parallel for
+ for (unsigned int i = 0; i < sz; ++i) {
+ this->at(i) /= e;
+ }
+
+}
+/* -------------------------------------------------------------------------- */
+
+template <typename T>
+void Map2d<T>::operator*=(const T & e){
+
+ UInt sz = this->n[0]*this->n[1];
+#pragma omp parallel for
+ for (unsigned int i = 0; i < sz; ++i) {
+ this->at(i) *= e;
+ }
+}
+
+/* -------------------------------------------------------------------------- */
+
+template <typename T>
+void Map2d<T>::operator/=(const T & e){
+
+ UInt sz = this->n[0]*this->n[1];
+#pragma omp parallel for
+ for (unsigned int i = 0; i < sz; ++i) {
+ this->at(i) /= e;
+ }
+
+}
+
+/* -------------------------------------------------------------------------- */
+
+template <typename T>
+template <template<typename> class surf, typename T2>
+void Map2d<T>::operator*=(const surf<T2> & e){
+
+ UInt sz = this->n[0]*this->n[1];
+#pragma omp parallel for
+ for (unsigned int i = 0; i < sz; ++i) {
+ this->at(i) *= e.at(i);
+ }
+
+}
+/* -------------------------------------------------------------------------- */
+
+template <typename T>
+template <template<typename> class surf, typename T2>
+void Map2d<T>::operator/=(const surf<T2> & e){
+
+ UInt sz = this->n[0]*this->n[1];
+#pragma omp parallel for
+ for (unsigned int i = 0; i < sz; ++i) {
+ this->at(i) /= e.at(i);
}
}
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
#endif /* __MAP_2D_HH__ */
diff --git a/tests/SConscript b/tests/SConscript
new file mode 100644
index 00000000..644f79dc
--- /dev/null
+++ b/tests/SConscript
@@ -0,0 +1,18 @@
+from os.path import join
+
+Import('main_env')
+
+test_files = Split("""
+run_tests.sh
+test_hertzian.py
+test_westergaard.py
+test_fractal.py
+test_fractal.verified
+""")
+
+src_dir = "#/tests"
+build_dir = 'build-' + main_env['build_type'] + '/tests'
+
+for file in test_files:
+ source = join(src_dir, file)
+ Command(file, source, Copy("$TARGET", "$SOURCE"))
diff --git a/tests/run_tests.sh b/tests/run_tests.sh
index 89bad37c..2aebe223 100755
--- a/tests/run_tests.sh
+++ b/tests/run_tests.sh
@@ -1,32 +1,33 @@
#!/bin/bash
n_tests=0
passed_tests=0
for test_file in test_*.py; do
((n_tests++))
echo "Testing $test_file"
echo -n "------------ "
python $test_file >/dev/null 2>&1
result=$?
if [ $result -ne 0 ]; then
echo "failure"
else
echo "success"
((passed_tests++))
fi
done
python test_fractal.py 2>/dev/null | tail -n 19 > test_fractal.out 2>&1
diff test_fractal.out test_fractal.verified >/dev/null 2>&1
result=$?
if [ $result -ne 0 ]; then
+ echo ""
echo "test_fractal.py diff failed"
((passed_tests--))
fi
rm test_fractal.out
failed_tests=$((n_tests - passed_tests))
echo ""
echo "Test results :"
echo " - passed $passed_tests/$n_tests"
echo " - failed $failed_tests/$n_tests"
diff --git a/tests/test_fftransform.py_no b/tests/test_fftransform.py_no
index bd579d0b..883d8e1e 100644
--- a/tests/test_fftransform.py_no
+++ b/tests/test_fftransform.py_no
@@ -1,55 +1,58 @@
#!/usr/bin/env python
# coding: utf-8
# -----------------------------------------------------------------------------
# @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/>.
# -----------------------------------------------------------------------------
import sys
-import tamaas as tm
import numpy as np
import matplotlib.pyplot as plt
+from util_test import import_tamaas
+
def main():
+ tm = import_tamaas()
+
size = 512
# Generate random rough surface
SG = tm.SurfaceGeneratorFilterFFT()
SG.getGridSize().assign(size)
SG.getHurst().assign(0.8)
SG.getRMS().assign(1.);
SG.getQ0().assign(4);
SG.getQ1().assign(4);
SG.getQ2().assign(32);
SG.getRandomSeed().assign(156);
SG.Init()
surface = SG.buildSurface()
# Compute fourier transform through numpy
numpy_spectral = np.fft.fft2(surface)
# Compute fourier tranform through tamaas
tamaas_spectral = np.zeros((size, size))
#transform = tm.FFTransformReal(size, surface, transform)
return 0
if __name__ == '__main__':
sys.exit(main())
diff --git a/tests/test_fractal.py b/tests/test_fractal.py
index 750b7e19..993ae2fa 100644
--- a/tests/test_fractal.py
+++ b/tests/test_fractal.py
@@ -1,100 +1,101 @@
#!/usr/bin/python
# -*- coding: utf-8 -*-
##*
#
# @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/>.
#
#
################################################################
-import sys
-from tamaas import *
+import sys, imp
+
+import tamaas as tm
def main():
#generate surface
- SG = SurfaceGeneratorFilterFFT()
+ SG = tm.SurfaceGeneratorFilterFFT()
SG.getGridSize().assign(512)
SG.getHurst().assign(0.8)
SG.getRMS().assign(1.);
SG.getQ0().assign(4);
SG.getQ1().assign(4);
SG.getQ2().assign(32);
SG.getRandomSeed().assign(156);
SG.Init()
a = SG.buildSurface()
#compute and print surface statistics
class Stats:
def __getitem__(self,key):
return self.__dict__[key]
stats = Stats()
stats.size = SG.getGridSize()
stats.hurst = SG.getHurst().value()
stats.rms = SG.getRMS()
stats.k0 = SG.getQ0()
stats.k1 = SG.getQ1().value()
stats.k2 = SG.getQ2().value()
stats.seed = SG.getRandomSeed()
- stats.rms_spectral = SurfaceStatistics.computeSpectralStdev(a);
- stats.rms_slopes_spectral = SurfaceStatistics.computeSpectralRMSSlope(a);
- stats.rms_geometric = SurfaceStatistics.computeStdev(SurfaceForPythonReal(a));
- stats.rms_slopes_geometric = SurfaceStatistics.computeRMSSlope(a);
- stats.moments = SurfaceStatistics.computeMoments(a);
+ stats.rms_spectral = tm.SurfaceStatistics.computeSpectralStdev(a);
+ stats.rms_slopes_spectral = tm.SurfaceStatistics.computeSpectralRMSSlope(a);
+ stats.rms_geometric = tm.SurfaceStatistics.computeStdev(tm.SurfaceForPythonReal(a));
+ stats.rms_slopes_geometric = tm.SurfaceStatistics.computeRMSSlope(a);
+ stats.moments = tm.SurfaceStatistics.computeMoments(a);
stats.m0 = stats['rms_spectral']**2
stats.m2 = stats.moments[0]
stats.m4 = stats.moments[1]
stats.alpha = stats['m0']*stats['m4']/(stats['m2']**2)
stats.L = 1.
- stats.m0prime = SurfaceStatistics.computeAnalyticFractalMoment(0,stats.k1,stats.k2,stats.hurst,1. , stats.L)
+ stats.m0prime = tm.SurfaceStatistics.computeAnalyticFractalMoment(0,stats.k1,stats.k2,stats.hurst,1. , stats.L)
stats.moment_A = stats.m0/stats.m0prime
- stats.analytic_m0 = SurfaceStatistics.computeAnalyticFractalMoment(0,stats.k1,stats.k2,stats.hurst,stats.moment_A,stats.L);
- stats.analytic_m2 = SurfaceStatistics.computeAnalyticFractalMoment(2,stats.k1,stats.k2,stats.hurst,stats.moment_A,stats.L);
- stats.analytic_m4 = SurfaceStatistics.computeAnalyticFractalMoment(4,stats.k1,stats.k2,stats.hurst,stats.moment_A,stats.L);
+ stats.analytic_m0 = tm.SurfaceStatistics.computeAnalyticFractalMoment(0,stats.k1,stats.k2,stats.hurst,stats.moment_A,stats.L);
+ stats.analytic_m2 = tm.SurfaceStatistics.computeAnalyticFractalMoment(2,stats.k1,stats.k2,stats.hurst,stats.moment_A,stats.L);
+ stats.analytic_m4 = tm.SurfaceStatistics.computeAnalyticFractalMoment(4,stats.k1,stats.k2,stats.hurst,stats.moment_A,stats.L);
stats.analytic_alpha = stats.analytic_m0*stats.analytic_m4/(stats.analytic_m2*stats.analytic_m2);
print """
[N] {size}
[rms] {rms}
[rmsSpectral] {rms_spectral}
[rmsSlopeSpectral] {rms_slopes_spectral}
[rmsSlopeGeometric] {rms_slopes_geometric}
[Hurst] {hurst}
[k1] {k1}
[k2] {k2}
[moment A] {moment_A}
[m0] {m0}
[analytic m0] {analytic_m0}
[m2] {m2}
[analytic m2] {analytic_m2}
[m4] {m4}
[analytic m4] {analytic_m4}
[alpha] {alpha}
[analytic_alpha] {analytic_alpha}
[seed] {seed}
""".format(**stats.__dict__)
return 0
if __name__ == "__main__":
sys.exit(main())
diff --git a/tests/test_hertzian.py b/tests/test_hertzian.py
index 7164ca17..d2bec8ec 100644
--- a/tests/test_hertzian.py
+++ b/tests/test_hertzian.py
@@ -1,121 +1,122 @@
#!/usr/bin/env python
# coding: utf-8
# -----------------------------------------------------------------------------
# @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/>.
# -----------------------------------------------------------------------------
import sys
-import tamaas as tm
import numpy as np
import matplotlib.pyplot as plt
+import tamaas as tm
+
def plotSurface(surface):
fig = plt.figure()
ax = fig.add_subplot(111)
img = ax.imshow(surface)
fig.colorbar(img)
def constructHertzProfile(size, curvature):
radius = 1. / curvature
x = np.linspace(-0.5, 0.5, size)
y = np.linspace(-0.5, 0.5, size)
x, y = np.meshgrid(x, y)
surface = np.sqrt(radius**2 - x**2 - y**2)
surface -= surface.mean()
return surface.copy()
def computeHertzDisplacement(e_star, contact_size, max_pressure, size):
x = np.linspace(-0.5, 0.5, size)
y = np.linspace(-0.5, 0.5, size)
x, y = np.meshgrid(x, y)
disp = np.pi * max_pressure / (4 * contact_size * e_star) * (2 * contact_size**2 - (x**2 + y**2))
return disp.copy()
def main():
grid_size = 1024
curvature = 0.1
effective_modulus = 1.
load = 0.0001
surface = constructHertzProfile(grid_size, curvature)
bem = tm.BemPolonski(surface)
bem.setEffectiveModulus(effective_modulus)
bem.computeEquilibrium(1e-12, load)
tractions = bem.getTractions()
displacements = bem.getDisplacements()
# Testing contact area against Hertz solution for solids of revolution
contact_area = tm.SurfaceStatistics.computeContactArea(tractions)
hertz_contact_size = (3 * load / (4 * curvature * effective_modulus))**(1. / 3.)
hertz_area = np.pi * hertz_contact_size**2
area_error = np.abs(hertz_area - contact_area) / hertz_area
print "Area error: {}".format(area_error)
# Testing maximum pressure
max_pressure = tractions.max()
hertz_max_pressure = (6 * load * effective_modulus**2 * curvature ** 2)**(1. / 3.) / np.pi
pressure_error = np.abs(hertz_max_pressure - max_pressure) / hertz_max_pressure
print "Max pressure error: {}".format(pressure_error)
# Testing displacements
hertz_displacements = computeHertzDisplacement(effective_modulus,
hertz_contact_size,
hertz_max_pressure,
grid_size)
# Selecing only the points that are in contact
contact_indexes = [(i, j, tractions[i, j] > 0) for i in range(grid_size) for j in range(grid_size)]
contact_indexes = map(lambda x: x[0:2], filter(lambda x: x[2], contact_indexes))
# Displacements of bem are centered around the mean of the whole surface
# and Hertz displacements are not centered, so we need to compute mean
# on the contact zone for both arrays
bem_mean = 0.
hertz_mean = 0.
for index in contact_indexes:
bem_mean += displacements[index]
hertz_mean += hertz_displacements[index]
bem_mean /= len(contact_indexes)
hertz_mean /= len(contact_indexes)
# Correction applied when computing error
correction = hertz_mean - bem_mean
# Computation of error of displacement in contact zone
error = 0.
hertz_norm = 0.
for index in contact_indexes:
error += (hertz_displacements[index] - displacements[index] - correction)**2
hertz_norm += (hertz_displacements[index] - hertz_mean)**2
displacement_error = np.sqrt(error / hertz_norm)
print "Displacement error (in contact zone): {}".format(displacement_error)
if area_error > 1e-3 or pressure_error > 3e-3 or displacement_error > 1e-4:
return 1
return 0
if __name__ == "__main__":
sys.exit(main())
diff --git a/tests/test_westergaard.py b/tests/test_westergaard.py
index e0a20fd6..88d7d3ea 100644
--- a/tests/test_westergaard.py
+++ b/tests/test_westergaard.py
@@ -1,76 +1,77 @@
#!/usr/bin/env python
# coding: utf-8
# -----------------------------------------------------------------------------
# @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/>.
# -----------------------------------------------------------------------------
import sys
-import tamaas as tm
import numpy as np
import matplotlib.pyplot as plt
+import tamaas as tm
+
def constructSinProfile(size, mode, amplitude):
x = np.linspace(0, 1, size)
y = np.linspace(0, 1, size)
x, y = np.meshgrid(x, y)
surface = amplitude * np.sin(2*np.pi*x*mode)
return surface.copy()
def main():
grid_size = 256
mode = 1
delta = 0.1
effective_modulus = 1.
p_star = np.pi * effective_modulus * delta * mode # Full contact load
load = 0.4 * p_star
surface = constructSinProfile(grid_size, mode, delta)
bem = tm.BemPolonski(surface)
bem.setEffectiveModulus(effective_modulus)
bem.computeEquilibrium(1e-12, load)
tractions = bem.getTractions()
displacements = bem.getDisplacements()
# Testing contact area against Westergaard solution
contact_area = tm.SurfaceStatistics.computeContactArea(tractions)
westergaard_contact_size = 2. / (mode * np.pi) * np.arcsin(np.sqrt(load / p_star))
area_error = np.abs(contact_area - mode * westergaard_contact_size) / (mode * westergaard_contact_size)
print "Area error: {}".format(area_error)
# Testing displacements at full contact
load = 1.1 * p_star
bem.computeEquilibrium(1e-12, load)
contact_area = tm.SurfaceStatistics.computeContactArea(tractions)
print "Area at load > full contact : {}".format(contact_area)
# Testing pressure distribution at full contact
westergaard_pressure = constructSinProfile(grid_size, mode, p_star)
tractions_amplitude = tractions.max() - tractions.mean()
return 0
if __name__ == "__main__":
sys.exit(main())