Homec4science
Diffusion tamaas (master) Browse

tamaas/src/model/dcfft.cpp
5187505f4535master

Last Change
Raw File

dcfft.cpp
View Options

/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2023 Lucas Frérot
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "dcfft.hh"
#include "model.hh"
#include <cmath>
/* -------------------------------------------------------------------------- */
namespace tamaas {
/// Square patch pressure solution, Johnson (1985) page 54
Real boussinesq(Real x, Real y, Real a, Real b) {
Real plus_sqrt = std::sqrt((x + a) * (x + a) + (y + b) * (y + b));
Real minu_sqrt = std::sqrt((x - a) * (x - a) + (y - b) * (y - b));
Real pm_sqrt = std::sqrt((x + a) * (x + a) + (y - b) * (y - b));
Real mp_sqrt = std::sqrt((x - a) * (x - a) + (y + b) * (y + b));
return 1. / M_PI *
((x + a) * std::log((y + b + plus_sqrt) / (y - b + pm_sqrt)) +
(y + b) * std::log((x + a + plus_sqrt) / (x - a + mp_sqrt)) +
(x - a) * std::log((y - b + minu_sqrt) / (y + b + mp_sqrt)) +
(y - b) * std::log((x - a + minu_sqrt) / (x + a + pm_sqrt)));
}
DCFFT::DCFFT(Model* model) : Westergaard(model) {
auto bdisc = model->getBoundaryDiscretization();
decltype(bdisc) extended_shape = {bdisc[0] * 2, bdisc[1] * 2};
auto extended_hermitian =
GridHermitian<Real, bdim>::hermitianDimensions(extended_shape);
buffer.resize(extended_hermitian);
extended_buffer.resize(extended_shape);
influence->resize(extended_hermitian);
initInfluence();
}
void DCFFT::initInfluence() {
auto system_size = model->getBoundarySystemSize();
auto bdisc = model->getBoundaryDiscretization();
Vector<Real, 2> metric{
{system_size[0] / bdisc[0], system_size[1] / bdisc[1]}};
Vector<Real, 2> pixel = metric * 0.5;
auto xshape = extended_buffer.sizes();
Grid<Real, bdim> rspace_influence(xshape, 1);
for (UInt m = 0; m < xshape[0]; ++m) {
for (UInt n = 0; n < xshape[1]; ++n) {
Vector<Int, 2> idx{{static_cast<Int>(m), static_cast<Int>(n)}};
Vector<Int, 2> over_half{{m > xshape[0] / 2, n > xshape[1] / 2}};
// Technique from https://doi.org/10.1088/2051-672X/ac860a page 13
for (UInt k = 0; k < 2; k++) {
idx(k) =
(1 - over_half(k)) * idx(k) - over_half(k) * (xshape[k] - idx(k));
}
Vector<Real, 2> x{{idx(0) * metric(0), idx(1) * metric(1)}};
rspace_influence(m, n) = boussinesq(x(0), x(1), pixel(0), pixel(1));
}
}
engine->forward(rspace_influence, *influence);
}
void DCFFT::apply(GridBase<Real>& input, GridBase<Real>& output) const {
const auto& shape = model->getBoundaryDiscretization();
Grid<Real, bdim> in{shape, 1, input.view()}, out{shape, 1, output.view()};
auto E_star = model->getHertzModulus();
extended_buffer = 0;
for (UInt i = 0; i < shape[0]; ++i)
for (UInt j = 0; j < shape[1]; ++j)
extended_buffer(i, j) = in(i, j);
Westergaard<model_type::basic_2d, IntegralOperator::neumann>::apply(
extended_buffer, extended_buffer);
for (UInt i = 0; i < shape[0]; ++i)
for (UInt j = 0; j < shape[1]; ++j)
out(i, j) = extended_buffer(i, j) / E_star;
}
} // namespace tamaas
c4science · Help