diff --git a/src/model/hooke.cpp b/src/model/hooke.cpp
index e018d2f..266a90c 100644
--- a/src/model/hooke.cpp
+++ b/src/model/hooke.cpp
@@ -1,94 +1,182 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2022 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 .
*
*/
/* -------------------------------------------------------------------------- */
#include "hooke.hh"
#include "influence.hh"
#include "loop.hh"
#include "model.hh"
#include "static_types.hh"
+#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
template
void Hooke::apply(GridBase& strain, GridBase& stress) const {
constexpr auto dim = model_type_traits::dimension;
const influence::ElasticHelper elasticity(
this->model->getShearModulus(), this->model->getPoissonRatio());
// Checking incompressibility
if (elasticity.nu == 0.5)
TAMAAS_EXCEPTION("Incompressibility error");
// Checking number of components
if (strain.getNbComponents() == dim * dim) {
Loop::loop(
[elasticity] CUDA_LAMBDA(MatrixProxy sigma,
MatrixProxy epsilon) {
sigma = elasticity(epsilon);
},
range>(stress),
range>(strain));
}
// in case we have voigt notation
else if (strain.getNbComponents() == voigt_size::value) {
Loop::loop(
[elasticity] CUDA_LAMBDA(SymMatrixProxy sigma,
SymMatrixProxy epsilon) {
sigma = elasticity(epsilon);
},
range>(stress),
range>(strain));
}
else
TAMAAS_EXCEPTION("Strain components do not match dimension");
}
/* -------------------------------------------------------------------------- */
template
std::pair Hooke::matvecShape() const {
auto shape = this->model->getDiscretization();
auto N =
std::accumulate(shape.begin(), shape.end(), 6, std::multiplies());
return std::make_pair(N, N);
}
template
GridBase Hooke::matvec(GridBase& X) const {
constexpr auto dim = model_type_traits::dimension;
X.setNbComponents(voigt_size::value);
GridBase x(X.getGlobalNbPoints(), X.getNbComponents());
apply(X, x);
return x;
}
/* -------------------------------------------------------------------------- */
-template class Hooke;
-template class Hooke;
-template class Hooke;
-template class Hooke;
-template class Hooke;
-template class Hooke;
+template
+HookeField::HookeField(Model* model) : Hooke(model) {
+ (*this)["mu"] = allocateGrid(*model, 1);
+ (*this)["nu"] = allocateGrid(*model, 1);
+
+ this->template field("mu") = model->getShearModulus();
+ this->template field("nu") = model->getPoissonRatio();
+}
+
+template
+struct HookeFieldHelper {
+ void operator()(MatrixProxy sigma,
+ MatrixProxy epsilon, const Real& mu,
+ const Real& nu) {
+ elasticity.mu = mu;
+ elasticity.nu = nu;
+ sigma = elasticity(epsilon);
+ }
+
+ influence::ElasticHelper elasticity{0, 0};
+};
+
+template
+struct SymHookeFieldHelper {
+ void operator()(SymMatrixProxy sigma,
+ SymMatrixProxy epsilon, const Real& mu,
+ const Real& nu) {
+ elasticity.mu = mu;
+ elasticity.nu = nu;
+ sigma = elasticity(epsilon);
+ }
+
+ influence::ElasticHelper elasticity{0, 0};
+};
+
+template
+void HookeField::apply(GridBase& strain,
+ GridBase& stress) const {
+ constexpr auto dim = model_type_traits::dimension;
+
+ // Checking incompressibility
+ auto incompressible =
+ Loop::reduce([](Real nu) -> UInt { return nu >= 0.5; },
+ this->template field("nu"));
+
+ if (incompressible > 0)
+ TAMAAS_EXCEPTION("Incompressibility error");
+
+ std::vector> hookes;
+
+ const auto& mu = this->template field("mu");
+ const auto& nu = this->template field("nu");
+
+ std::transform(mu.begin(), mu.end(), nu.begin(), std::back_inserter(hookes),
+ [](const Real& mu, const Real& nu) {
+ return influence::ElasticHelper(mu, nu);
+ });
+
+ // Checking number of components
+ if (strain.getNbComponents() == dim * dim) {
+ Loop::loop(
+ [](MatrixProxy stress,
+ MatrixProxy strain,
+ const influence::ElasticHelper& elasticity) {
+ stress = elasticity(strain);
+ },
+ range>(stress),
+ range>(strain), hookes);
+ }
+
+ // in case we have voigt notation
+ else if (strain.getNbComponents() == voigt_size::value) {
+ Loop::loop(
+ [](SymMatrixProxy stress,
+ SymMatrixProxy strain,
+ const influence::ElasticHelper& elasticity) {
+ stress = elasticity(strain);
+ },
+ range>(stress),
+ range>(strain), hookes);
+ }
+
+ else
+ TAMAAS_EXCEPTION("Strain components do not match dimension");
+}
+
+/* -------------------------------------------------------------------------- */
+#define INSTANTIATE_HOOKE(r, _, type) \
+ template class Hooke; \
+ template class HookeField;
+
+BOOST_PP_SEQ_FOR_EACH(INSTANTIATE_HOOKE, ~, TAMAAS_MODEL_TYPES);
+#undef INSTANTIATE_HOOKE
} // namespace tamaas
diff --git a/src/model/hooke.hh b/src/model/hooke.hh
index 864b576..cc8031b 100644
--- a/src/model/hooke.hh
+++ b/src/model/hooke.hh
@@ -1,63 +1,72 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2022 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 .
*
*/
/* -------------------------------------------------------------------------- */
#ifndef HOOKE_HH
#define HOOKE_HH
/* -------------------------------------------------------------------------- */
#include "integral_operator.hh"
#include "model_type.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/**
* @brief Applies Hooke's law of elasticity
*/
template
class Hooke : public IntegralOperator {
public:
using IntegralOperator::IntegralOperator;
/// Type of underlying model
model_type getType() const override { return type; }
/// Hooke's law computes stresses ~ Neumann
IntegralOperator::kind getKind() const override {
return IntegralOperator::dirac;
}
/// Does not update
void updateFromModel() override {}
/// Apply Hooke's tensor
void apply(GridBase& strain, GridBase& stress) const override;
/// LinearOperator shape
std::pair matvecShape() const override;
/// LinearOperator interface
GridBase matvec(GridBase& X) const override;
};
+template
+class HookeField : public Hooke {
+public:
+ HookeField(Model* model);
+
+ /// Apply Hooke's tensor
+ void apply(GridBase& strain, GridBase& stress) const override;
+};
+
} // namespace tamaas
#endif // HOOKE_HH