Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F84219460
mindlin.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Sep 21, 10:55
Size
7 KB
Mime Type
text/x-c++
Expires
Mon, Sep 23, 10:55 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20961559
Attached To
rTAMAAS tamaas
mindlin.cpp
View Options
/**
* @file
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2017 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 "mindlin.hh"
#include "elasto_plastic_model.hh"
#include "influence.hh"
#include "kelvin_helper.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
namespace detail {
template <model_type type, typename T>
class MindlinBoussinesqHelper {
using boussinesq_t = T;
using trait = model_type_traits<type>;
static constexpr UInt dim = trait::dimension;
static constexpr UInt bdim = trait::boundary_dimension;
public:
MindlinBoussinesqHelper(const T& boussinesq,
const influence::ElasticHelper<dim>& el, Real y_3,
Real cutoff)
: boussinesq(boussinesq), el(el), y_3(y_3), cutoff(cutoff) {}
inline void applyIntegral(GridHermitian<Real, bdim>& out,
GridHermitian<Real, bdim>& surface_gradients,
GridHermitian<Real, bdim>& source_stresses,
const Grid<Real, bdim>& wavevectors) const {
using k_out = typename KelvinTrait<T>::out_t;
using m_type = MatrixProxy<Complex, dim, dim>;
using v_type = VectorProxy<const Real, bdim>;
Loop::loop(
[&](typename KelvinTrait<T>::out_t out_local,
MatrixProxy<Complex, dim, dim> surface_gradient,
MatrixProxy<Complex, dim, dim> source_stress,
VectorProxy<const Real, bdim> q) {
// No cutoff
// if (-q.l2norm() * std::abs(y_3) < std::log(cutoff))
// return;
Vector<Real, dim> normal{{{0, 0, 1}}};
auto traction =
el(surface_gradient) * normal; // computing strains from gradu
traction -= source_stress * normal; // subtracting source stresses
out_local +=
boussinesq.applyU0(traction, q) *
influence::KelvinIntegrator<0>::g0<true>(q.l2norm() * y_3);
out_local +=
(boussinesq.applyU1(traction, q) *
influence::KelvinIntegrator<0>::g1<true>(q.l2norm() * y_3));
},
range<k_out>(out), range<m_type>(surface_gradients),
range<m_type>(source_stresses), range<v_type>(wavevectors));
}
inline void applyConstantTerm(typename KelvinTrait<T>::out_t& out,
Vector<Complex, dim>& /*source*/,
const influence::ElasticHelper<dim>& /*el*/) const {
out = 0;
}
protected:
const T& boussinesq;
const influence::ElasticHelper<dim>& el;
const Real y_3, cutoff;
};
template <>
inline void
MindlinBoussinesqHelper<model_type::volume_2d, influence::Boussinesq<3, 1>>::
applyConstantTerm(
typename KelvinTrait<influence::Boussinesq<3, 1>>::out_t& out,
Vector<Complex, dim>& source,
const influence::ElasticHelper<dim>& el) const {
out = 0;
out(2, 0) = -source(0) / el.mu;
out(2, 1) = -source(1) / el.mu;
out(2, 2) = -source(2) / (el.lambda + 2 * el.mu);
}
} // namespace detail
template <model_type type, UInt order>
void Mindlin<type, order>::applyIf(GridBase<Real>& source, GridBase<Real>& out,
filter_t pred) const {
Real nu = this->model->getPoissonRatio(), mu = this->model->getShearModulus();
constexpr UInt derivative = order - 2;
influence::Kelvin<trait::dimension, derivative> kelvin(mu, nu);
influence::Kelvin<trait::dimension, 2> kelvin_strain(mu, nu);
influence::Boussinesq<trait::dimension, derivative - 1> boussinesq(mu, nu);
influence::ElasticHelper<trait::dimension> elasticity(mu, nu);
auto apply = [&](UInt node, decltype(this->source_buffers)& source_buffers,
decltype(this->disp_buffer)& out_buffer) {
const Real L = this->model->getSystemSize().front();
const UInt N = this->model->getDiscretization().front();
const Real dl = L / (N - 1);
const Real xi = node * dl;
// this should only be used when node == xi == 0
detail::LinearElement<type, decltype(kelvin_strain)> strain_element{
thrust::make_pair(node, xi), kelvin_strain, source_buffers,
surface_strains, this->wavevectors};
detail::LinearElement<type, decltype(kelvin)> element{
thrust::make_pair(node, xi), kelvin, source_buffers, out_buffer,
this->wavevectors};
out_buffer = 0;
// Simple condition to ensure correct integration
auto filter_out_elem = [&pred](UInt e) {
return not pred(e) and not pred(e + 1);
};
/// Compute surface strains only once (dirty)
if (node == 0) {
surface_strains = 0;
for (UInt e : Loop::range(N - 1)) {
if (filter_out_elem(e))
continue;
const auto node_positions = thrust::make_pair(e * dl, (e + 1) * dl);
strain_element.integrateElement(e, node_positions);
}
}
// Apply Kelvin operator
for (UInt e : Loop::range(N - 1)) {
if (filter_out_elem(e))
continue;
const auto node_positions = thrust::make_pair(e * dl, (e + 1) * dl);
element.integrateElement(e, node_positions);
}
// Correcting for the tractions on the surface
detail::MindlinBoussinesqHelper<type, decltype(boussinesq)> helper(
boussinesq, elasticity, xi, 1e-3);
helper.applyIntegral(out_buffer, surface_strains, source_buffers.front(),
this->wavevectors);
typename detail::KelvinTrait<decltype(boussinesq)>::out_t out_fundamental(
out_buffer(0));
typename detail::KelvinTrait<decltype(kelvin)>::source_t source_fundamental(
source_buffers[node](0));
Vector<Real, trait::dimension> e3;
e3 = 0;
e3(trait::dimension - 1) = -1;
auto in_fundamental = source_fundamental * e3;
// Uniform shift in case of gradient computation
helper.applyConstantTerm(out_fundamental, in_fundamental, elasticity);
};
this->fourierApplyIf(apply, source, out, pred);
}
/* -------------------------------------------------------------------------- */
/* Template instanciation */
/* -------------------------------------------------------------------------- */
template class Mindlin<model_type::volume_2d, 3>;
template class Mindlin<model_type::volume_2d, 4>;
/* -------------------------------------------------------------------------- */
} // namespace tamaas
Event Timeline
Log In to Comment