Homec4science
Diffusion tamaas (ci-fix-static-checks) Browse

tamaas/src/bem_grid.cpp
303235694413ci-fix-static-checks

Last Change
Raw File

bem_grid.cpp
View Options

/**
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Tamaas is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "bem_grid.hh"
#include "grid.hh"
#include "types.hh"
#include "fft_plan_manager.hh"
#include <cmath>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
BemGrid::BemGrid(Surface<Real> & 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};
spectral_tractions.setNbComponents(3);
spectral_displacements.setNbComponents(3);
spectral_tractions.resize(hermitian_sizes);
spectral_displacements.resize(hermitian_sizes);
}
/* -------------------------------------------------------------------------- */
BemGrid::~BemGrid() {
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;
}
/* -------------------------------------------------------------------------- */
void BemGrid::computeSurfaceNormals() {
}
/* -------------------------------------------------------------------------- */
void BemGrid::computeDisplacementsFromTractions() {
applyInfluenceFunctions(tractions, displacements);
}
void BemGrid::computeTractionsFromDisplacements() {
applyInverseInfluenceFunctions(tractions, displacements);
}
/* -------------------------------------------------------------------------- */
void BemGrid::computeGaps(Grid<Real, 2> &disp) {
#pragma omp parallel for collapse(2)
for (UInt i = 0 ; i < surface.sizes()[0] ; i++) {
for (UInt j = 0 ; j < surface.sizes()[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<Real, 2> & input, Grid<Real, 2> & output) {
TAMAAS_ASSERT(input.dataSize() == output.dataSize(), "input and output have different size");
FFTransform<Real, 2> & plan_1 = FFTPlanManager::get().createPlan(input, spectral_tractions);
FFTransform<Real, 2> & 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);
VectorProxy<complex> spectral_p(NULL, 3);
VectorProxy<complex> spectral_u(NULL, 3);
MatrixProxy<complex> 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);
spectral_u.mul(spectral_F, spectral_p);
++it_p;
++it_u;
}
plan_2.backwardTransform();
}
/* -------------------------------------------------------------------------- */
void BemGrid::applyInverseInfluenceFunctions(Grid<Real, 2> & __attribute__((unused)) input,
Grid<Real, 2> & __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<complex> it(influence.getInternalData(), 3, 3);
// Imaginary unit
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 = std::sqrt(q_1*q_1+q_2*q_2); \
complex influence_1[] = INFLUENCE(q_1, q_2, q, E, nu); \
MatrixProxy<complex> 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 = 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<complex> 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)
// Loops over parts of surface
/*
[ ][ ][ ][ ]
[ ][x][ ][ ]
[ ][ ][ ][ ]
[ ][ ][ ][o]
*/
for (UInt i = 1 ; i < size/2 ; i++) {
for (UInt j = i ; j < size/2 ; j++) {
FILL_INFLUENCE(i, j);
}
}
/*
[ ][ ][ ][ ]
[ ][ ][ ][o]
[ ][ ][ ][ ]
[ ][x][ ][ ]
*/
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][ ][ ][ ]
*/
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][ ]
*/
for (UInt i = size/2+1 ; i < size ; i++) {
UInt j = size/2;
FILL_INFLUENCE(i-size, j);
}
/* 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__
c4science · Help