diff --git a/src/bem_grid.cpp b/src/bem_grid.cpp index c61e32f..43ea6f6 100644 --- a/src/bem_grid.cpp +++ b/src/bem_grid.cpp @@ -1,254 +1,308 @@ /** * * @author Lucas Frérot * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "bem_grid.hh" #include "grid.hh" #include "types.hh" #include "fft_plan_manager.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ BemGrid::BemGrid(Surface & surface): E(1.), nu(0.3), surface(surface), surface_normals(surface.sizes(), surface.lengths(), 3), - influence(surface.sizes(), surface.lengths(), 9), tractions(surface.sizes(), surface.lengths(), 3), displacements(surface.sizes(), surface.lengths(), 3), gaps(surface.sizes(), surface.lengths(), 3), true_tractions(NULL), true_displacements(NULL) { UInt hermitian_sizes[2] = {surface.sizes()[0], surface.sizes()[1]/2+1}; + influence.setNbComponents(9); spectral_tractions.setNbComponents(3); spectral_displacements.setNbComponents(3); + influence.resize(hermitian_sizes); spectral_tractions.resize(hermitian_sizes); spectral_displacements.resize(hermitian_sizes); } /* -------------------------------------------------------------------------- */ BemGrid::~BemGrid() { + FFTPlanManager::get().clean(); delete true_tractions; delete true_displacements; } /* -------------------------------------------------------------------------- */ const Grid2dReal & BemGrid::getTrueTractions() const { if (true_tractions == NULL) TAMAAS_EXCEPTION("True tractions were not allocated"); return *true_tractions; } const Grid2dReal & BemGrid::getTrueDisplacements() const { if (true_displacements == NULL) TAMAAS_EXCEPTION("True displacements were not allocated"); return *true_displacements; } /* -------------------------------------------------------------------------- */ +/// This function computes surface slopes in fourier space and fills the normals grid void BemGrid::computeSurfaceNormals() { + SurfaceComplex spectral_surface(surface.size(), surface.getL()); + GridHermitian spectral_slopes(spectral_surface.sizes(), + spectral_surface.lengths(), + 3); + FFTransform & plan1 = FFTPlanManager::get().createPlan(surface, spectral_surface); + FFTransform & plan2 = FFTPlanManager::get().createPlan(surface_normals, spectral_slopes); + + plan1.forwardTransform(); + + const UInt * n = spectral_surface.sizes(); + constexpr complex I(0, 1); + +#pragma omp parallel +{ +#pragma omp for collapse(2) + for (UInt i = 0 ; i <= n[0] / 2 ; i++) { + for (UInt j = 0 ; j < n[1] ; j++) { + Real q1 = 2*M_PI*j; + Real q2 = 2*M_PI*i; + spectral_slopes(i, j, 0) = -q1*I*spectral_surface(i, j); + spectral_slopes(i, j, 1) = -q2*I*spectral_surface(i, j); + spectral_slopes(i, j, 2) = 0; + } + } + +#pragma omp for collapse(2) + for (UInt i = n[0]/2 + 1 ; i < n[0] ; i++) { + for (UInt j = 0 ; j < n[1] ; j++) { + Real q1 = 2*M_PI*j; + Real q2 = 2*M_PI*(i-n[0]); + spectral_slopes(i, j, 0) = -q1*I*spectral_surface(i, j); + spectral_slopes(i, j, 1) = -q2*I*spectral_surface(i, j); + spectral_slopes(i, j, 2) = 0; + } + } +} // end #pragma omp parallel + + // Dirac + spectral_slopes(0, 0, 2) = n[0]*n[0]; + + plan2.backwardTransform(); + FFTPlanManager::get().destroyPlan(surface, spectral_surface); + FFTPlanManager::get().destroyPlan(surface_normals, spectral_slopes); } /* -------------------------------------------------------------------------- */ void BemGrid::computeDisplacementsFromTractions() { applyInfluenceFunctions(tractions, displacements); } void BemGrid::computeTractionsFromDisplacements() { applyInverseInfluenceFunctions(tractions, displacements); } /* -------------------------------------------------------------------------- */ void BemGrid::computeGaps(Grid &disp) { const UInt * n = surface.sizes(); #pragma omp parallel for collapse(2) for (UInt i = 0 ; i < n[0] ; i++) { for (UInt j = 0 ; j < n[1] ; j++) { gaps(i, j, 0) = disp(i, j, 0); gaps(i, j, 1) = disp(i, j, 1); gaps(i, j, 2) = disp(i, j, 2) - surface(i, j); } } } /* -------------------------------------------------------------------------- */ void BemGrid::applyInfluenceFunctions(Grid & input, Grid & output) { TAMAAS_ASSERT(input.dataSize() == output.dataSize(), "input and output have different size"); FFTransform & plan_1 = FFTPlanManager::get().createPlan(input, spectral_tractions); FFTransform & plan_2 = FFTPlanManager::get().createPlan(output, spectral_displacements); plan_1.forwardTransform(); - auto it_p = spectral_tractions.begin(3); - auto it_u = spectral_displacements.begin(3); + complex * F = influence.getInternalData(), * p = spectral_tractions.getInternalData(), + * u = spectral_displacements.getInternalData(); VectorProxy spectral_p(NULL, 3); VectorProxy spectral_u(NULL, 3); MatrixProxy spectral_F(NULL, 3, 3); - // OpenMP not very cool with this: must have only one iteration variable in for-loop -#pragma omp parallel for firstprivate(it_p, it_u, spectral_p, spectral_u, spectral_F) - for (auto it_F = influence.begin(9); it_F < influence.end(9); ++it_F) { - spectral_p.setPointer(*it_p); - spectral_u.setPointer(*it_u); - spectral_F.setPointer(*it_F); + const UInt N = influence.dataSize() / influence.getNbComponents(); + + // OpenMP not very cool with this: best to have a POD iteration variable +#pragma omp parallel for firstprivate(spectral_p, spectral_u, spectral_F) + for (UInt i = 0 ; i < N ; ++i) { + spectral_p.setPointer(p + i * spectral_p.dataSize()); + spectral_u.setPointer(u + i * spectral_u.dataSize()); + spectral_F.setPointer(F + i * spectral_F.dataSize()); spectral_u.mul(spectral_F, spectral_p); - ++it_p; - ++it_u; } plan_2.backwardTransform(); } /* -------------------------------------------------------------------------- */ void BemGrid::applyInverseInfluenceFunctions(Grid & __attribute__((unused)) input, Grid & __attribute__((unused)) output) { TAMAAS_EXCEPTION("applyInverseInfluenceFunctions not yet implemented"); } /* -------------------------------------------------------------------------- */ // warning: this function will suck your soul out void BemGrid::computeInfluence() { // Iterator over points of surface MatrixProxy it(influence.getInternalData(), 3, 3); // Imaginary unit - complex I(0, 1); + constexpr complex I(0, 1); const UInt size = this->surface.size(); this->influence = 0; // Influence functions matrix, ask Lucas #define INFLUENCE(q_1, q_2, q, E, nu) { \ 1./(E*q*q*q)*(-2*nu*nu*q_1*q_1 - 2*nu*q_2*q_2 + 2*q_1*q_1 + 2*q_2*q_2), \ -q_1*q_2/(E*q*q*q)*(1+nu)*2*nu, \ I*q_1/(E*q*q)*(1+nu)*(1-2*nu), \ -q_1*q_2/(E*q*q*q)*(1+nu)*2*nu, \ 1./(E*q*q*q)*(-2*nu*nu*q_2*q_2 - 2*nu*q_1*q_1 + 2*q_1*q_1 + 2*q_2*q_2), \ I*q_2/(E*q*q)*(1+nu)*(1-2*nu), \ -I*q_1/(E*q*q)*(1+nu)*(1-2*nu), \ -I*q_2/(E*q*q)*(1+nu)*(1-2*nu), \ 2./(E*q)*(1-nu*nu)} // Fill influence for single point #define FILL_INFLUENCE(k, l) do { \ - Real q_1 = 2*M_PI*(k); \ - Real q_2 = 2*M_PI*(l); \ + Real q_1 = 2*M_PI*(l); \ + Real q_2 = 2*M_PI*(k); \ Real q = std::sqrt(q_1*q_1+q_2*q_2); \ complex influence_1[] = INFLUENCE(q_1, q_2, q, E, nu); \ MatrixProxy inf_mat(influence_1, 3, 3); \ it.setPointer(&influence(i, j, 0)); \ it = inf_mat; \ } while(0) // Fill influence points for "symetric" areas #define FILL_SYM_INFLUENCE(k, l) do { \ - Real q_1 = 2*M_PI*(k); \ - Real q_2 = 2*M_PI*(l); \ + Real q_1 = 2*M_PI*(l); \ + Real q_2 = 2*M_PI*(k); \ Real q = std::sqrt(q_1*q_1+q_2*q_2); \ complex influence_1[] = INFLUENCE(q_1, q_2, q, E, nu); \ complex influence_2[] = INFLUENCE(q_2, q_1, q, E, nu); \ MatrixProxy inf_mat(influence_1, 3, 3); \ it.setPointer(&influence(i, j, 0)); \ it = inf_mat; \ inf_mat.setPointer(influence_2); \ it.setPointer(&influence(j, i, 0)); \ it = inf_mat; \ } while(0) +#pragma omp parallel firstprivate(it) // all the following loops are independent +{ + // Loops over parts of surface /* [ ][ ][ ][ ] [ ][x][ ][ ] [ ][ ][ ][ ] [ ][ ][ ][o] */ +#pragma omp for // can't collapse triangluar loops for (UInt i = 1 ; i < size/2 ; i++) { for (UInt j = i ; j < size/2 ; j++) { FILL_INFLUENCE(i, j); } } /* [ ][ ][ ][ ] [ ][ ][ ][o] [ ][ ][ ][ ] [ ][x][ ][ ] */ +#pragma omp for collapse(2) for (UInt i = size/2+1 ; i < size ; i++) { for (UInt j = 0 ; j < size/2 ; j++) { FILL_INFLUENCE(i, j-size); } } /* [ ][x][x][o] [x][ ][x][ ] [x][x][x][ ] [o][ ][ ][ ] */ +#pragma omp for for (UInt i = 1 ; i <= size/2 ; i++) { { UInt j = size/2; FILL_SYM_INFLUENCE(i, j); } { UInt j = 0; FILL_SYM_INFLUENCE(i, j); } } /* [ ][ ][ ][ ] [ ][ ][ ][ ] [ ][ ][ ][o] [ ][ ][x][ ] */ +#pragma omp for for (UInt i = size/2+1 ; i < size ; i++) { UInt j = size/2; FILL_INFLUENCE(i-size, j); } +} // End #pragma omp parallel /* State of surface now: [0][x][x][o] [x][x][x][o] [x][x][x][x] [o][x][x][o] Legend: x => explicitely filled o => hermitian symetry 0 => actually 0 */ } /* -------------------------------------------------------------------------- */ __END_TAMAAS__ diff --git a/src/fftransform.cpp b/src/fftransform.cpp index 51daf0b..ed14bac 100644 --- a/src/fftransform.cpp +++ b/src/fftransform.cpp @@ -1,55 +1,55 @@ /** * * @author Lucas Frérot * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "fftransform.hh" __BEGIN_TAMAAS__ template FFTransform::FFTransform(Grid & real, GridHermitian & spectral): real(real), spectral(spectral) {} /* -------------------------------------------------------------------------- */ template FFTransform::~FFTransform() {} /* -------------------------------------------------------------------------- */ template void FFTransform::normalize() { - real /= real.dataSize(); + real /= real.dataSize() / (Real)real.getNbComponents(); } /* -------------------------------------------------------------------------- */ template class FFTransform; __END_TAMAAS__ diff --git a/src/types.cpp b/src/types.cpp index 555eb0d..027b7f4 100644 --- a/src/types.cpp +++ b/src/types.cpp @@ -1,133 +1,145 @@ /** * * @author Lucas Frérot * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "types.hh" +#include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /* -------------------------------------------------------------------------- */ /* VectorProxy definitions */ /* -------------------------------------------------------------------------- */ template VectorProxy::VectorProxy(T * mem, UInt n) : Proxy() { this->data.wrapMemory(mem, n); this->n[0] = n; } /* -------------------------------------------------------------------------- */ template VectorProxy::~VectorProxy() {} /* -------------------------------------------------------------------------- */ template const T & VectorProxy::operator()(UInt i) const { return Grid::operator()(i, 0); } /* -------------------------------------------------------------------------- */ template T & VectorProxy::operator()(UInt i) { return Grid::operator()(i, 0); } /* -------------------------------------------------------------------------- */ template void VectorProxy::mul(const MatrixProxy &mat, const VectorProxy &vec) { const UInt * n_mat = mat.sizes(); const UInt * n_vec = vec.sizes(); TAMAAS_ASSERT(n_mat[1] == n_vec[0], "matrix and vector sizes do not match"); TAMAAS_ASSERT(n_mat[0] == this->n[0], "matrix and result vector sizes do not match"); // No parallelism here for (UInt i = 0 ; i < n_mat[0] ; ++i) { for (UInt j = 0 ; j < n_vec[0] ; ++j) { (*this)(i) += mat(i, j)*vec(j); } } } /* -------------------------------------------------------------------------- */ template T VectorProxy::dot(const VectorProxy &vec) { TAMAAS_ASSERT(this->n[0] == vec.sizes()[0], "vector sizes do not match"); T res = 0; // No parallelism here for (UInt i = 0 ; i < this->n[0] ; i++) res += (*this)(i) * vec(i); return res; } +/* -------------------------------------------------------------------------- */ + +template +T VectorProxy::l2norm() const { + T res = 0; + // No parallelism here + for (UInt i = 0 ; i < this->n[0] ; i++) + res += (*this)(i)*(*this)(i); + return std::sqrt(res); +} + /* -------------------------------------------------------------------------- */ /* MatrixProxy definitions */ /* -------------------------------------------------------------------------- */ template MatrixProxy::MatrixProxy(T * mem, UInt n, UInt m) : Proxy() { this->data.wrapMemory(mem, n*m); this->n[0] = n; this->n[1] = m; } /* -------------------------------------------------------------------------- */ template MatrixProxy::~MatrixProxy() {} /* -------------------------------------------------------------------------- */ template const T & MatrixProxy::operator()(UInt i, UInt j) const { return Grid::operator()(i, j, 0); } /* -------------------------------------------------------------------------- */ template T & MatrixProxy::operator()(UInt i, UInt j) { return Grid::operator()(i, j, 0); } /* -------------------------------------------------------------------------- */ #define TYPE_INSTANCIATE(type) template class VectorProxy; \ template class MatrixProxy TYPE_INSTANCIATE(Real); TYPE_INSTANCIATE(UInt); TYPE_INSTANCIATE(complex); TYPE_INSTANCIATE(int); TYPE_INSTANCIATE(bool); TYPE_INSTANCIATE(unsigned long); #undef TYPE_INSTANCIATE __END_TAMAAS__ diff --git a/src/types.hh b/src/types.hh index d59c26e..8c1c5a4 100644 --- a/src/types.hh +++ b/src/types.hh @@ -1,91 +1,104 @@ /** * * @author Lucas Frérot * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __TYPES_HH__ #define __TYPES_HH__ /* -------------------------------------------------------------------------- */ #include "tamaas.hh" #include "grid.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ template class Proxy : public Grid { public: /// Default constructor Proxy() : Grid() {} + /// Copy constructor + Proxy(Proxy & p) : Grid() { + if (p.n != this->n) + memcpy(this->n, p.n, dim * sizeof(UInt)); + this->nb_components = p.nb_components; + this->data.wrapMemory(p.getInternalData(), p.dataSize()); + } /// Destructor virtual ~Proxy() {} public: /// Set pointer void setPointer(T * ptr) {this->data.setPtr(ptr);} /// Get increment ptrdiff_t increment() const {return this->data.size();} }; // Forward declarations template class MatrixProxy; template class VectorProxy : public Proxy { public: /// Wrap on memory VectorProxy(T * mem, UInt n); + /// Copy constructor + VectorProxy(VectorProxy & v): Proxy(v) {} /// Destructor virtual ~VectorProxy(); public: /// Access const T & operator()(UInt i) const; /// Access (non-const) T & operator()(UInt i); /// Matrix-vector multiplication void mul(const MatrixProxy & mat, const VectorProxy & vec); /// Dot product T dot(const VectorProxy & vec); + /// L2 norm + T l2norm() const; }; /* -------------------------------------------------------------------------- */ template class MatrixProxy : public Proxy { public: /// Wrap on memory MatrixProxy(T * mem, UInt n, UInt m); + /// Copy constructor + MatrixProxy(MatrixProxy & v): Proxy(v) {} /// Destructor virtual ~MatrixProxy(); public: /// Access const T & operator()(UInt i, UInt j) const; /// Access (non-const) T & operator()(UInt i, UInt j); }; __END_TAMAAS__ #endif // __TYPES_HH__ diff --git a/tests/SConscript b/tests/SConscript index 24010b2..36f8042 100644 --- a/tests/SConscript +++ b/tests/SConscript @@ -1,23 +1,24 @@ from os.path import join Import('main_env') test_files = Split(""" run_tests.sh test_hertz_pressure.py test_westergaard.py test_fractal.py test_fractal.verified test_hertz_disp.py test_hertz_kato.py test_hertz_adhesion.py test_saturated_pressure.py test_fftransform.py +test_bem_grid.py """) 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/test_bem_grid.py b/tests/test_bem_grid.py new file mode 100644 index 0000000..6f1410f --- /dev/null +++ b/tests/test_bem_grid.py @@ -0,0 +1,114 @@ +#!/usr/bin/env python +# coding: utf-8 +# ----------------------------------------------------------------------------- +# @author Lucas Frérot +# +# @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 . +# ----------------------------------------------------------------------------- + +import sys +import tamaas as tm +import numpy as np + +from numpy.linalg import norm +from numpy.fft import fft2, ifft2 + +def make_profile(f, profile, k): + x = np.linspace(0, 1, profile.shape[0], endpoint=False) + y = np.linspace(0, 1, profile.shape[1], endpoint=False) + x, y = np.meshgrid(x, y) + profile[:,:] = f(2*k*np.pi*x) + +def sin_profile(profile, k): + make_profile(np.sin, profile, k) + +def cos_profile(profile, k): + make_profile(np.cos, profile, k) + + +def main(): + """Checking that all the pipes are clean""" + tm.initialize() + + size = 256 + k = 9 + + E = 0.91 + nu = 0.3 + mu = E / (2*(1+nu)) + + pressure = np.zeros((size, size)) + surface = np.zeros((size, size)) + westergaard_disp = np.zeros((size, size)) + + cos_profile(pressure, k) + cos_profile(surface, k) + bem = tm.BemGridPolonski(surface) + bem.setElasticity(E, nu) + bem.computeInfluence() + + bem.getTractions()[:, :, 2] = pressure[:, :] + bem.computeDisplacementsFromTractions() + + displacements = bem.getDisplacements() + + # Constructing vertical displacement solution + cos_profile(westergaard_disp, k) + westergaard_disp *= (1-nu**2)/(E*np.pi*k) + westergaard_norm = (1-nu**2)/(E*np.pi*k)*np.sqrt(0.5) * size + # Computing error + error = norm(displacements[:, :, 2] - westergaard_disp)/westergaard_norm + + if error > 1e-10: + print "Error in normal displacement: {}".format(error) + tm.finalize() + return 1 + + # Constructing horizontal displacement solution + sin_profile(westergaard_disp, k) + westergaard_disp *= (2*nu-1)/(4*np.pi*k*mu) + westergaard_norm = np.abs(2*nu-1)/(4*np.pi*k*mu)*np.sqrt(0.5)*size + error = norm(displacements[:, :, 0] - westergaard_disp)/westergaard_norm + + if error > 1e-10: + print "Error in horizontal displacement: {}".format(error) + tm.finalize() + return 1 + + # Looking at normal vectors + normals = np.zeros((size, size, 3)) + sin_profile(normals[:, :, 0], k) + normals[:, :, 0] *= 2*np.pi*k + normals[:, :, 2] = 1 + + bem.computeSurfaceNormals() + bem_normals = bem.getSurfaceNormals() + + error = norm(bem_normals - normals) + if error > 1e-9: + print "Error in surface normals: {}".format(error) + tm.finalize() + return 1 + + tm.finalize() + return 0 + +if __name__ == "__main__": + sys.exit(main())