Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/bem_grid.cpp
562a3d65fd33features/netcdf-mpi

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),
tractions(surface.sizes(), surface.lengths(), 3),
displacements(surface.sizes(), surface.lengths(), 3),
gaps(surface.sizes(), surface.lengths(), 3),
true_tractions(NULL),
true_displacements(NULL),
dump_freq(100)
{
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);
//computeSurfaceNormals();
}
/* -------------------------------------------------------------------------- */
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<Real> spectral_surface(surface.size(), surface.getL());
GridHermitian<Real, 2> spectral_slopes(spectral_surface.sizes(),
spectral_surface.lengths(),
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
const Int * n = reinterpret_cast<const Int * >(spectral_surface.sizes());
constexpr 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) {
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<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> & __attribute__((unused)) input,
Grid<Real, 2> & __attribute__((unused)) output) {
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
constexpr complex I(0, 1);
const Int size = this->surface.size();
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); \
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<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); \
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<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 += std::pow(std::norm(m.l2norm()), 2);
}
lipschitz = std::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;
const UInt * 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;
const UInt * 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