Homec4science
Diffusion tamaas (features/declarative) Browse

tamaas/src/bem/bem_grid.cpp
cd034e49e5c8features/declarative

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 "fft_plan_manager.hh"
#include "grid.hh"
#include "types.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
BemGrid::BemGrid(Surface<Real>& surface)
: E(1.), nu(0.3), surface(), surface_normals(surface.sizes(), 3),
tractions(surface.sizes(), 3), displacements(surface.sizes(), 3),
gaps(surface.sizes(), 3), true_tractions(NULL), true_displacements(NULL),
dump_freq(100) {
this->surface.wrap(surface);
auto hermitian_sizes =
GridHermitian<Real, 2>::hermitianDimensions(surface.sizes());
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);
// computeSurfaceNormals();
}
/* -------------------------------------------------------------------------- */
BemGrid::~BemGrid() {
FFTPlanManager::get().clean();
delete true_tractions;
delete true_displacements;
}
/* -------------------------------------------------------------------------- */
const Grid<Real, 2>& BemGrid::getTrueTractions() const {
if (true_tractions == NULL)
TAMAAS_EXCEPTION("True tractions were not allocated");
return *true_tractions;
}
const Grid<Real, 2>& 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<Real> spectral_surface(surface.size(), surface.getL());
GridHermitian<Real, 2> spectral_slopes(spectral_surface.sizes(), 3);
FFTransform<Real, 2>& plan1 =
FFTPlanManager::get().createPlan(surface, spectral_surface);
FFTransform<Real, 2>& plan2 =
FFTPlanManager::get().createPlan(surface_normals, spectral_slopes);
plan1.forwardTransform();
// Loops should be with signed integers
std::array<Int, 2> n;
std::copy(spectral_surface.sizes().begin(), spectral_surface.sizes().end(),
n.begin());
const Complex I(0, 1);
#pragma omp parallel
{
#pragma omp for collapse(2)
for (Int i = 0; i <= n[0] / 2; i++) {
for (Int 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 (Int i = n[0] / 2 + 1; i < n[0]; i++) {
for (Int 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);
// Make surface normal
VectorProxy<Real> sn(nullptr, 3);
const UInt N = surface_normals.getNbPoints();
#pragma omp parallel for firstprivate(sn)
for (UInt i = 0; i < N; i++) {
sn.setPointer(surface_normals.getInternalData() + i * sn.dataSize());
sn /= sn.l2norm();
}
}
/* -------------------------------------------------------------------------- */
void BemGrid::computeDisplacementsFromTractions() {
applyInfluenceFunctions(tractions, displacements);
}
void BemGrid::computeTractionsFromDisplacements() {
applyInverseInfluenceFunctions(tractions, displacements);
}
/* -------------------------------------------------------------------------- */
void BemGrid::computeGaps(Grid<Real, 2>& disp) {
auto 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<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();
Complex *F = influence.getInternalData(),
*p = spectral_tractions.getInternalData(),
*u = spectral_displacements.getInternalData();
VectorProxy<Complex> spectral_p(nullptr, 3), spectral_u(nullptr, 3);
MatrixProxy<Complex> spectral_F(nullptr, 3, 3);
spectral_displacements = 0;
const UInt N = influence.getNbPoints();
// 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);
}
plan_2.backwardTransform();
}
/* -------------------------------------------------------------------------- */
void BemGrid::applyInverseInfluenceFunctions(
Grid<Real, 2>& input[[gnu::unused]], Grid<Real, 2>& output[[gnu::unused]]) {
TAMAAS_EXCEPTION("applyInverseInfluenceFunctions not yet implemented");
}
/* -------------------------------------------------------------------------- */
// warning: this function will suck your soul out
// loops must be written with signed integers
void BemGrid::computeInfluence() {
// Iterator over points of surface
MatrixProxy<Complex> it(influence.getInternalData(), 3, 3);
// Imaginary unit
const Complex I(0, 1);
const Int size = this->surface.size();
const Real L_1 = 1.; // this->surface.lengths()[0];
const Real L_2 = 1.; // this->surface.lengths()[1];
this->influence = 0;
// Influence functions matrix, ask Lucas
#define INFLUENCE(q_1, q_2, q, E, nu) \
{ \
2. / (E * q * q * q) * (-nu * nu * q_1 * q_1 - nu * q_2 * q_2 + q * q), \
-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, \
2. / (E * q * q * q) * \
(-nu * nu * q_2 * q_2 - nu * q_1 * q_1 + q * q), \
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 * (l) / L_1; \
Real q_2 = 2 * M_PI * (k) / L_2; \
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 * (l) / L_1; \
Real q_2 = 2 * M_PI * (k) / L_2; \
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)
#pragma omp parallel firstprivate( \
it) // all the following loops are independent
{
// Loops over parts of surface
/*
[ ][ ][ ][ ]
[ ][x][ ][ ]
[ ][ ][ ][ ]
[ ][ ][ ][o]
*/
#pragma omp for schedule(dynamic) // can't collapse triangluar loops
for (Int i = 1; i < size / 2; i++) {
for (Int j = i; j < size / 2; j++) {
FILL_SYM_INFLUENCE(i, j);
}
}
/*
[ ][ ][ ][ ]
[ ][ ][ ][o]
[ ][ ][ ][ ]
[ ][x][ ][ ]
*/
#pragma omp for collapse(2)
for (Int i = size / 2 + 1; i < size; i++) {
for (Int j = 1; j < size / 2; j++) {
FILL_INFLUENCE(i - size, j);
}
}
/*
[ ][x][x][o]
[x][ ][x][ ]
[x][x][x][ ]
[ ][ ][ ][ ]
*/
#pragma omp for
for (Int i = 1; i <= size / 2; i++) {
{
Int j = size / 2;
FILL_SYM_INFLUENCE(i, j);
}
{
Int j = 0;
FILL_SYM_INFLUENCE(i, j);
}
}
/*
[ ][ ][ ][ ]
[ ][ ][ ][ ]
[ ][ ][ ][o]
[x][ ][x][ ]
*/
#pragma omp for
for (Int i = size / 2 + 1; i < size; i++) {
{
Int j = size / 2;
FILL_INFLUENCE(i - size, j);
}
{
Int j = 0;
FILL_INFLUENCE(i - size, j);
}
}
} // End #pragma omp parallel
/* State of surface now:
[0][x][x][o]
[x][x][x][o]
[x][x][x][o]
[x][x][x][o]
Legend: x => explicitely filled
o => hermitian symetry
0 => actually 0
*/
computeInfluenceLipschitz();
#undef INFLUENCE
#undef FILL_INFLUENCE
#undef FILL_SYM_INFLUENCE
}
/* -------------------------------------------------------------------------- */
/// Computing L_2 norm of influence functions
void BemGrid::computeInfluenceLipschitz() {
const UInt N = influence.getNbPoints();
VectorProxy<Complex> m(nullptr, 9);
Real _norm = 0.;
#pragma omp parallel for firstprivate(m) reduction(+ : _norm)
for (UInt i = 0; i < N; i++) {
m.setPointer(influence.getInternalData() + i * m.dataSize());
_norm += pow(norm(m.l2norm()), 2);
}
lipschitz = sqrt(_norm);
}
/* -------------------------------------------------------------------------- */
/* Friction related computations */
/* -------------------------------------------------------------------------- */
void BemGrid::projectOnFrictionCone(Grid<Real, 2>& field, Real mu) {
VectorProxy<Real> p_proxy(nullptr, 3);
const UInt N = field.getNbPoints();
#pragma omp parallel for firstprivate(p_proxy)
for (UInt i = 0; i < N; i++) {
p_proxy.setPointer(field.getInternalData() + i * p_proxy.dataSize());
projectOnFrictionCone(p_proxy, mu);
}
}
/* -------------------------------------------------------------------------- */
void BemGrid::enforceMeanValue(Grid<Real, 2>& field,
const Grid<Real, 1>& target) {
Real red_0 = 0, red_1 = 0, red_2 = 0;
auto n = field.sizes();
#pragma omp parallel for collapse(2) reduction(+ : red_0, red_1, red_2)
for (UInt i = 0; i < n[0]; i++) {
for (UInt j = 0; j < n[1]; j++) {
red_0 += field(i, j, 0);
red_1 += field(i, j, 1);
red_2 += field(i, j, 2);
}
}
red_0 /= n[0] * n[1];
red_1 /= n[0] * n[1];
red_2 /= n[0] * n[1];
Real p1_data[3] = {red_0, red_1, red_2};
// Loop for pressure shifting
VectorProxy<Real> p1(p1_data, 3);
VectorProxy<Real> p0((Real*)target.getInternalData(), 3);
VectorProxy<Real> p(nullptr, 3);
const UInt N = field.getNbPoints();
#pragma omp parallel for firstprivate(p0, p1)
for (UInt i = 0; i < N; i++) {
p.setPointer(field.getInternalData() + i * p.dataSize());
p -= p1;
p += p0;
}
}
/* -------------------------------------------------------------------------- */
void BemGrid::computeMeanValueError(Grid<Real, 2>& field,
const Grid<Real, 1>& target,
Grid<Real, 1>& error_vec) {
Real p0 = 0, p1 = 0, p2 = 0;
auto n = field.sizes();
#pragma omp parallel for collapse(2) reduction(+ : p0, p1, p2)
for (UInt i = 0; i < n[0]; i++) {
for (UInt j = 0; j < n[1]; j++) {
p0 += field(i, j, 0);
p1 += field(i, j, 1);
p2 += field(i, j, 2);
}
}
p0 /= n[0] * n[1];
p1 /= n[0] * n[1];
p2 /= n[0] * n[1];
Real error[3] = {p0 - target(0), p1 - target(1), p2 - target(2)};
VectorProxy<Real> e(error, 3);
error_vec = e;
}
/* -------------------------------------------------------------------------- */
Real BemGrid::computeMeanValueError(Grid<Real, 2>& field,
const Grid<Real, 1>& target) {
Real e[3] = {0.};
VectorProxy<Real> error(e, 3);
computeMeanValueError(field, target, error);
Real norm =
target(0) * target(0) + target(1) * target(1) + target(2) * target(2);
return error.l2norm() / std::sqrt(norm);
}
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
c4science · Help