diff --git a/python/wrap/model.cpp b/python/wrap/model.cpp
index 9457459..41ba139 100644
--- a/python/wrap/model.cpp
+++ b/python/wrap/model.cpp
@@ -1,280 +1,284 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 "model.hh"
#include "adhesion_functional.hh"
#include "functional.hh"
#include "integral_operator.hh"
#include "model_dumper.hh"
#include "model_extensions.hh"
#include "model_factory.hh"
#include "numpy.hh"
#include "residual.hh"
#include "wrap.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace wrap {
using namespace py::literals;
/// Wrap functional classes
void wrapFunctionals(py::module& mod) {
py::class_ func(
mod, "Functional");
func.def(py::init<>())
.def("computeF", &functional::Functional::computeF)
.def("computeGradF", &functional::Functional::computeGradF);
py::class_ adh(mod, "AdhesionFunctional",
func);
adh.def_property("parameters", &functional::AdhesionFunctional::getParameters,
&functional::AdhesionFunctional::setParameters)
// legacy wrapper code
.def("setParameters", &functional::AdhesionFunctional::setParameters);
py::class_(
mod, "ExponentialAdhesionFunctional", adh)
.def(py::init&>(), "surface"_a);
py::class_(
mod, "MaugisAdhesionFunctional", adh)
.def(py::init&>(), "surface"_a);
py::class_(
mod, "SquaredExponentialAdhesionFunctional", adh)
.def(py::init&>(), "surface"_a);
}
template
std::unique_ptr> instanciateFromNumpy(numpy& num) {
std::unique_ptr> result = nullptr;
switch (num.ndim()) {
case 2:
result = std::make_unique>>(num);
return result;
case 3:
result = std::make_unique>>(num);
return result;
case 4:
result = std::make_unique>>(num);
return result;
default:
TAMAAS_EXCEPTION("instanciateFromNumpy expects the last dimension of numpy "
"array to be the number of components");
}
}
/// Wrap IntegralOperator
void wrapIntegralOperator(py::module& mod) {
py::class_(mod, "IntegralOperator")
.def("apply",
[](IntegralOperator& op, numpy input, numpy output) {
auto in = instanciateFromNumpy(input);
auto out = instanciateFromNumpy(output);
op.apply(*in, *out);
})
.def("updateFromModel", &IntegralOperator::updateFromModel)
.def("getModel", &IntegralOperator::getModel,
py::return_value_policy::reference)
.def("getKind", &IntegralOperator::getKind)
.def("getType", &IntegralOperator::getType);
py::enum_(mod, "integration_method")
.value("linear", integration_method::linear)
.value("cutoff", integration_method::cutoff);
}
/// Wrap BEEngine classes
void wrapBEEngine(py::module& mod) {
py::class_(mod, "BEEngine")
.def("solveNeumann", &BEEngine::solveNeumann)
.def("solveDirichlet", &BEEngine::solveDirichlet)
.def("getModel", &BEEngine::getModel, py::return_value_policy::reference)
.def("registerNeumann", &BEEngine::registerNeumann)
.def("registerDirichlet", &BEEngine::registerDirichlet);
}
/// Wrap Models
void wrapModelClass(py::module& mod) {
py::enum_(mod, "model_type")
.value("basic_1d", model_type::basic_1d)
.value("basic_2d", model_type::basic_2d)
.value("surface_1d", model_type::surface_1d)
.value("surface_2d", model_type::surface_2d)
.value("volume_1d", model_type::volume_1d)
.value("volume_2d", model_type::volume_2d);
py::class_(mod, "Model")
.def_property_readonly("type", &Model::getType)
.def("setElasticity", &Model::setElasticity, "E"_a, "nu"_a)
.def_property("E", &Model::getYoungModulus, &Model::setYoungModulus)
.def_property("nu", &Model::getPoissonRatio, &Model::setPoissonRatio)
.def("getHertzModulus", &Model::getHertzModulus)
.def("getYoungModulus", &Model::getYoungModulus)
.def("getShearModulus", &Model::getShearModulus)
.def("getPoissonRatio", &Model::getPoissonRatio)
.def("getTraction", (GridBase & (Model::*)()) & Model::getTraction,
py::return_value_policy::reference_internal)
.def("getDisplacement",
(GridBase & (Model::*)()) & Model::getDisplacement,
py::return_value_policy::reference_internal)
.def("getSystemSize", &Model::getSystemSize)
.def("getDiscretization", &Model::getDiscretization)
.def("getBoundarySystemSize", &Model::getBoundarySystemSize)
.def("getBoundaryDiscretization", &Model::getBoundaryDiscretization)
.def("solveNeumann", &Model::solveNeumann)
.def("solveDirichlet", &Model::solveDirichlet)
.def("dump", &Model::dump)
.def("addDumper", &Model::addDumper, "dumper"_a, py::keep_alive<1, 2>())
.def("setDumper",
[](Model&, const ModelDumper&) {
TAMAAS_EXCEPTION("setDumper() is not a member of Model; use "
"addDumper() instead");
})
.def("getBEEngine", &Model::getBEEngine,
py::return_value_policy::reference_internal)
.def("getIntegralOperator", &Model::getIntegralOperator,
"operator_name"_a, py::return_value_policy::reference_internal)
.def("registerField",
[](Model& m, std::string name, numpy field) {
auto f = instanciateFromNumpy(field);
m.registerField(name, std::move(f));
},
"field_name"_a, "field"_a, py::keep_alive<1, 3>())
.def("getField",
(GridBase & (Model::*)(const std::string&)) & Model::getField,
"field_name"_a, py::return_value_policy::reference_internal)
.def("getFields", &Model::getFields)
.def("applyElasticity",
[](Model& model, numpy stress, numpy strain) {
auto out = instanciateFromNumpy(stress);
auto in = instanciateFromNumpy(strain);
model.applyElasticity(*out, *in);
})
.def("__repr__",
[](const Model& m) {
std::stringstream ss;
ss << m;
return ss.str();
})
.def("__getitem__",
[](const Model& m, std::string key) -> const GridBase& {
try {
return m.getField(key);
} catch (std::out_of_range&) {
throw py::key_error();
}
},
py::return_value_policy::reference_internal)
.def("__contains__",
[](const Model& m, std::string key) {
const auto fields = m.getFields();
return std::find(fields.begin(), fields.end(), key) !=
fields.end();
},
py::keep_alive<0, 1>())
.def("__iter__",
[](const Model& m) {
const auto& fields = m.getFieldsMap();
return py::make_key_iterator(fields.cbegin(), fields.cend());
},
py::keep_alive<0, 1>());
py::class_>(
mod, "ModelDumper")
.def(py::init<>())
.def("dump", &ModelDumper::dump, "model"_a)
.def("__lshift__",
[](ModelDumper& dumper, Model& model) { dumper << model; });
}
/// Wrap factory for models
void wrapModelFactory(py::module& mod) {
py::class_(mod, "ModelFactory")
.def_static("createModel", &ModelFactory::createModel, "model_type"_a,
"system_size"_a, "discretization"_a)
.def_static("createResidual", &ModelFactory::createResidual,
"model_type"_a, "sigma_y"_a, "hardening"_a = 0.)
.def_static("registerVolumeOperators",
&ModelFactory::registerVolumeOperators, "model"_a);
}
/// Wrap residual class
void wrapResidual(py::module& mod) {
// TODO adapt to n-dim
py::class_(mod, "Residual")
.def(py::init())
.def("computeResidual",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.computeResidual(*in);
})
.def("computeStress",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.computeStress(*in);
})
.def("updateState",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.updateState(*in);
})
.def("computeResidualDisplacement",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.computeResidualDisplacement(*in);
})
.def("applyTangent",
[](Residual& res, numpy& output, numpy& input,
numpy& current_strain_inc) {
auto out = instanciateFromNumpy(output);
auto in = instanciateFromNumpy(input);
auto inc = instanciateFromNumpy(current_strain_inc);
res.applyTangent(*out, *in, *inc);
},
"output"_a, "input"_a, "current_strain_increment"_a)
.def("getVector", &Residual::getVector,
py::return_value_policy::reference_internal)
.def("getPlasticStrain", &Residual::getPlasticStrain,
py::return_value_policy::reference_internal)
.def("getStress", &Residual::getStress,
py::return_value_policy::reference_internal)
.def("setIntegrationMethod", &Residual::setIntegrationMethod, "method"_a,
- "cutoff"_a = 1e-12);
+ "cutoff"_a = 1e-12)
+ .def_property("yield_stress", &Residual::getYieldStress,
+ &Residual::setYieldStress)
+ .def_property("hardening_modulus", &Residual::getHardeningModulus,
+ &Residual::setHardeningModulus);
}
void wrapModel(py::module& mod) {
wrapBEEngine(mod);
wrapModelClass(mod);
wrapModelFactory(mod);
wrapFunctionals(mod);
wrapResidual(mod);
wrapIntegralOperator(mod);
}
} // namespace wrap
} // namespace tamaas
diff --git a/python/wrap/model_extensions.hh b/python/wrap/model_extensions.hh
index 1a90d10..6b9d39b 100644
--- a/python/wrap/model_extensions.hh
+++ b/python/wrap/model_extensions.hh
@@ -1,116 +1,132 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 "functional.hh"
#include "model.hh"
#include "model_dumper.hh"
#include "residual.hh"
#include "wrap.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
namespace functional {
namespace wrap {
/// Class for extension of Functional in python
class PyFunctional : public Functional {
public:
using Functional::Functional;
// Overriding pure virtual functions
Real computeF(GridBase& variable, GridBase& dual) const override {
PYBIND11_OVERLOAD_PURE(Real, Functional, computeF, variable, dual);
}
void computeGradF(GridBase& variable,
GridBase& gradient) const override {
PYBIND11_OVERLOAD_PURE(void, Functional, computeGradF, variable, gradient);
}
};
} // namespace wrap
} // namespace functional
/* -------------------------------------------------------------------------- */
namespace wrap {
/* -------------------------------------------------------------------------- */
class PyModelDumper : public ModelDumper {
public:
using ModelDumper::ModelDumper;
void dump(const Model& model) override {
PYBIND11_OVERLOAD_PURE(void, ModelDumper, dump, model);
}
};
/* -------------------------------------------------------------------------- */
class PyResidual : public Residual {
public:
using Residual::Residual;
void computeResidual(GridBase& strain_increment) override {
PYBIND11_OVERLOAD_PURE(void, Residual, computeResidual, strain_increment);
}
void applyTangent(GridBase& output, GridBase& input,
GridBase& current_strain_increment) override {
PYBIND11_OVERLOAD_PURE(void, Residual, applyTangent, output, input,
current_strain_increment);
}
void computeStress(GridBase& strain_increment) override {
PYBIND11_OVERLOAD_PURE(void, Residual, computeStress, strain_increment);
}
void updateState(GridBase& converged_strain_increment) override {
PYBIND11_OVERLOAD_PURE(void, Residual, updateState,
converged_strain_increment);
}
const GridBase& getVector() const override {
PYBIND11_OVERLOAD_PURE(const GridBase&, Residual, getVector);
}
const GridBase& getPlasticStrain() const override {
PYBIND11_OVERLOAD_PURE(const GridBase&, Residual, getPlasticStrain);
}
const GridBase& getStress() const override {
PYBIND11_OVERLOAD_PURE(const GridBase&, Residual, getStress);
}
void computeResidualDisplacement(GridBase& strain_increment) override {
PYBIND11_OVERLOAD_PURE(void, Residual, computeResidualDisplacement,
strain_increment);
}
void setIntegrationMethod(integration_method method, Real cutoff) override {
PYBIND11_OVERLOAD_PURE(void, Residual, setIntegrationMethod, method,
cutoff);
}
+
+ Real getYieldStress() const override {
+ PYBIND11_OVERLOAD_PURE(Real, Residual, getYieldStress);
+ }
+
+ Real getHardeningModulus() const override {
+ PYBIND11_OVERLOAD_PURE(Real, Residual, getHardeningModulus);
+ }
+
+ void setYieldStress(Real val) override {
+ PYBIND11_OVERLOAD_PURE(void, Residual, setYieldStress, val);
+ }
+
+ void setHardeningModulus(Real val) override {
+ PYBIND11_OVERLOAD_PURE(void, Residual, setHardeningModulus, val);
+ }
};
/* -------------------------------------------------------------------------- */
} // namespace wrap
/* -------------------------------------------------------------------------- */
} // namespace tamaas
diff --git a/src/model/elasto_plastic/isotropic_hardening.hh b/src/model/elasto_plastic/isotropic_hardening.hh
index f255cf8..f6d521e 100644
--- a/src/model/elasto_plastic/isotropic_hardening.hh
+++ b/src/model/elasto_plastic/isotropic_hardening.hh
@@ -1,165 +1,175 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 ISOTROPIC_HARDENING_HH
#define ISOTROPIC_HARDENING_HH
/* -------------------------------------------------------------------------- */
#include "grid.hh"
#include "influence.hh"
#include "model.hh"
#include "model_type.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
template
class IsotropicHardening {
using trait = model_type_traits;
static constexpr UInt dim = trait::dimension;
public:
/// Constructor
IsotropicHardening(Model* model, Real sigma_0, Real h);
/// Compute plastic strain increment with radial return algorithm
template
void computePlasticIncrement(Grid& increment,
const Grid& strain,
const Grid& strain_increment);
/// Compute stress
template
void computeStress(Grid& stress, const Grid& strain,
const Grid& strain_increment);
void applyTangentIncrement(Grid& output,
const Grid& input,
const Grid& strain,
const Grid& strain_increment) const;
Real hardening(Real p) const { return sigma_0 + h * p; }
- TAMAAS_ACCESSOR(h, Real, HardeningModulus);
- TAMAAS_ACCESSOR(sigma_0, Real, YieldStress);
+ Real getHardeningModulus() const { return h; }
+ Real getYieldStress() const { return sigma_0; }
+ void setHardeningModulus(Real h_) {
+ if (h_ < 0)
+ TAMAAS_EXCEPTION("Hardening modulus should be positive");
+ h = h_;
+ }
+ void setYieldStress(Real sigma_0_) {
+ if (sigma_0_ < 0)
+ TAMAAS_EXCEPTION("Yield stress should be positive");
+ sigma_0 = sigma_0_;
+ }
const GridBase& getPlasticStrain() const { return *plastic_strain; }
GridBase& getPlasticStrain() { return *plastic_strain; }
protected:
Model* model;
Real sigma_0; /// < initial yield stress
Real h; /// < hardening modulus
std::shared_ptr> plastic_strain, cumulated_plastic_strain;
};
/* -------------------------------------------------------------------------- */
template <>
template
void IsotropicHardening::computePlasticIncrement(
Grid& increment, const Grid& strain,
const Grid& strain_increment) {
influence::ElasticHelper elasticity(this->model->getShearModulus(),
this->model->getPoissonRatio());
using pmatrix = SymMatrixProxy;
using cpmatrix = SymMatrixProxy;
Loop::loop(
[&elasticity,
this](auto dep, // < plastic strain increment
auto epsilon, // < total strain
auto delta_epsilon, // < strain increment
auto ep, // < plastic strain
Real& p /* < cumulated plastic strain */) {
// Trial elastic stress
auto sigma_tr = elasticity(epsilon - ep + delta_epsilon);
// Zero increment by default
dep = 0;
// Deviatoric trial stress
decltype(sigma_tr) dev;
dev.deviatoric(sigma_tr);
// Von-Mises stress
Real von_mises = std::sqrt(1.5) * dev.l2norm();
// Plasticity condition
if (von_mises - this->hardening(p) > 0) {
// radial return
Real dp = (von_mises - this->hardening(p)) / (3 * elasticity.mu + h);
dev *= 3 * dp / (2 * von_mises); // saving memory (dev is delta ep)
dep = dev;
if (update) { // constexpr if when C++17
p += dp;
ep += dep;
}
}
},
range(increment), range(strain),
range(strain_increment), range(*this->plastic_strain),
*this->cumulated_plastic_strain);
}
/* -------------------------------------------------------------------------- */
template <>
template
void IsotropicHardening::computeStress(
Grid& stress, const Grid& strain,
const Grid& strain_increment) {
using pmatrix = SymMatrixProxy;
using cpmatrix = SymMatrixProxy;
const influence::ElasticHelper elasticity(
this->model->getShearModulus(), this->model->getPoissonRatio());
this->template computePlasticIncrement(stress, strain,
strain_increment);
Loop::loop(
[&elasticity](auto sigma, // < stress
auto epsilon, // < total strain
auto delta_epsilon, // < strain increment
auto ep // < plastic strain
) {
// !!! sigma contains delta_ep
if (not update) {
- // dep is not contained in ep
+ // dep is not contained in ep
sigma *= -1;
sigma += epsilon;
} else {
- // dep is contained in ep, can discard value of sigma
+ // dep is contained in ep, can discard value of sigma
sigma = epsilon;
}
sigma -= ep;
sigma += delta_epsilon;
sigma = elasticity(sigma);
},
range(stress), range(strain),
range(strain_increment),
range(*this->plastic_strain));
}
} // namespace tamaas
/* -------------------------------------------------------------------------- */
#endif // ISOTROPIC_HARDENING_HH
diff --git a/src/model/elasto_plastic/residual.hh b/src/model/elasto_plastic/residual.hh
index 75aca15..4da4727 100644
--- a/src/model/elasto_plastic/residual.hh
+++ b/src/model/elasto_plastic/residual.hh
@@ -1,138 +1,156 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 RESIDUAL_HH
#define RESIDUAL_HH
/* -------------------------------------------------------------------------- */
#include "boussinesq.hh"
#include "isotropic_hardening.hh"
#include "mindlin.hh"
#include "model_type.hh"
/* -------------------------------------------------------------------------- */
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
/**
* @brief Residual manager
*/
class Residual {
public:
/// Constructor
Residual(Model* model) : model(model) {}
/// Destructor
virtual ~Residual() = default;
// Pure virtual methods
public:
/// Compute the residual vector for a given strain increment
virtual void computeResidual(GridBase& strain_increment) = 0;
/// Apply tangent
virtual void applyTangent(GridBase& output, GridBase& input,
GridBase& current_strain_increment) = 0;
/// Compute the stresses for a given strain increment
virtual void computeStress(GridBase& strain_increment) = 0;
/// Update the plastic state
virtual void updateState(GridBase& converged_strain_increment) = 0;
/// Get residual vector
virtual const GridBase& getVector() const = 0;
/// Get plastic strain
virtual const GridBase& getPlasticStrain() const = 0;
/// Get stresses
virtual const GridBase& getStress() const = 0;
/// Compute displacement
virtual void
computeResidualDisplacement(GridBase& strain_increment) = 0;
/// Set integration method (cutoff parameter ignored if linear integration)
virtual void setIntegrationMethod(integration_method method,
Real cutoff = 0) = 0;
// Accessors
public:
Model& getModel() { return *model; }
+ virtual Real getYieldStress() const = 0;
+ virtual void setYieldStress(Real sigma_y) = 0;
+ virtual Real getHardeningModulus() const = 0;
+ virtual void setHardeningModulus(Real h) = 0;
protected:
Model* model;
};
/**
* @brief Templated residual manager
*/
template
class ResidualTemplate : public Residual {
using trait = model_type_traits;
static constexpr UInt dim = trait::dimension;
public:
ResidualTemplate(Model* model, Real sigma_0, Real h);
// Implementation
public:
/// Compute the residual vector for a given strain increment
void computeResidual(GridBase& strain_increment) override;
/// Apply tangent
void applyTangent(GridBase& output, GridBase& input,
GridBase& current_strain_increment) override;
/// Compute the stresses for a given strain increment
void computeStress(GridBase& strain_increment) override;
/// Update the plastic state
void updateState(GridBase& converged_strain_increment) override;
/// Compute displacement
void computeResidualDisplacement(GridBase& strain_increment) override;
/// Set integration method (cutoff parameter ignored if linear integration)
void setIntegrationMethod(integration_method method,
Real cutoff = 0) override;
// Accessors
public:
/// Get residual vector
const GridBase& getVector() const override { return *residual; }
/// Get plastic strain
const GridBase& getPlasticStrain() const override {
return hardening.getPlasticStrain();
}
/// Get stresses
const GridBase& getStress() const override { return *stress; }
+ Real getYieldStress() const override { return hardening.getYieldStress(); }
+
+ void setYieldStress(Real sigma_y) override {
+ hardening.setYieldStress(sigma_y);
+ }
+
+ Real getHardeningModulus() const override {
+ return hardening.getHardeningModulus();
+ }
+
+ void setHardeningModulus(Real h) override {
+ hardening.setHardeningModulus(h);
+ }
+
private:
/// Convenience function
const IntegralOperator& integralOperator(const std::string& name) const {
return *this->model->getIntegralOperator(name);
}
/// Add non-zero layers of plastic strain into the filter
void updateFilter(Grid& plastic_strain_increment);
protected:
IsotropicHardening hardening;
std::shared_ptr> strain, stress, residual, tmp;
std::unordered_set plastic_layers;
std::function plastic_filter;
};
/* -------------------------------------------------------------------------- */
} // namespace tamaas
#endif // RESIDUAL_HH
diff --git a/src/model/model.cpp b/src/model/model.cpp
index 15788a3..87eaeff 100644
--- a/src/model/model.cpp
+++ b/src/model/model.cpp
@@ -1,163 +1,163 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 "model.hh"
#include "be_engine.hh"
#include "logger.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
Model::Model(std::vector system_size, std::vector discretization)
: system_size(std::move(system_size)),
discretization(std::move(discretization)) {}
/* -------------------------------------------------------------------------- */
Model::~Model() = default;
/* -------------------------------------------------------------------------- */
-void Model::setElasticity(Real E, Real nu) {
- this->E = E;
- this->nu = nu;
+void Model::setElasticity(Real E_, Real nu_) {
+ setYoungModulus(E_);
+ setPoissonRatio(nu_);
updateOperators();
}
/* -------------------------------------------------------------------------- */
Real Model::getHertzModulus() const { return E / (1 - nu * nu); }
/* -------------------------------------------------------------------------- */
GridBase& Model::getTraction() { return getField("traction"); }
const GridBase& Model::getTraction() const {
return getField("traction");
}
/* -------------------------------------------------------------------------- */
GridBase& Model::getDisplacement() { return getField("displacement"); }
const GridBase& Model::getDisplacement() const {
return getField("displacement");
}
/* -------------------------------------------------------------------------- */
const std::vector& Model::getSystemSize() const { return system_size; }
/* -------------------------------------------------------------------------- */
const std::vector& Model::getDiscretization() const {
return discretization;
}
/* -------------------------------------------------------------------------- */
void Model::solveNeumann() {
engine->registerNeumann();
engine->solveNeumann(getField("traction"), getField("displacement"));
}
void Model::solveDirichlet() {
engine->registerDirichlet();
engine->solveDirichlet(getField("displacement"), getField("traction"));
}
/* -------------------------------------------------------------------------- */
IntegralOperator* Model::getIntegralOperator(const std::string& name) {
return operators[name].get();
}
/* -------------------------------------------------------------------------- */
void Model::updateOperators() {
for (auto& op : operators)
op.second->updateFromModel();
}
/* -------------------------------------------------------------------------- */
void Model::addDumper(std::shared_ptr dumper) {
this->dumpers.push_back(std::move(dumper));
}
void Model::dump() const {
for (auto& dumper : dumpers)
if (dumper)
dumper->dump(*this);
}
/* -------------------------------------------------------------------------- */
void Model::registerField(const std::string& name,
std::shared_ptr> field) {
fields[name] = std::move(field);
}
const GridBase& Model::getField(const std::string& name) const try {
return *fields.at(name);
} catch (std::out_of_range& e) {
Logger().get(LogLevel::warning)
<< "Field " << name << " not registered in model\n";
throw e;
}
GridBase& Model::getField(const std::string& name) try {
return *fields.at(name);
} catch (std::out_of_range& e) {
Logger().get(LogLevel::warning)
<< "Field " << name << " not registered in model\n";
throw e;
}
std::vector Model::getFields() const {
std::vector keys;
keys.reserve(fields.size());
std::transform(fields.begin(), fields.end(), std::back_inserter(keys),
[](auto&& pair) { return pair.first; });
return keys;
}
/* -------------------------------------------------------------------------- */
std::ostream& operator<<(std::ostream& o, const Model& _this) {
o << "Model<" << _this.getType() << "> (E = " << _this.getYoungModulus()
<< ", nu = " << _this.getPoissonRatio() << ")\n";
// Printing domain size
o << " - domain = [";
std::for_each(_this.getSystemSize().begin(), _this.getSystemSize().end() - 1,
[&o](const Real& x) { o << x << ", "; });
o << _this.getSystemSize().back() << "]\n";
// Printing discretization
o << " - discretization = [";
std::for_each(_this.getDiscretization().begin(),
_this.getDiscretization().end() - 1,
[&o](const UInt& x) { o << x << ", "; });
o << _this.getDiscretization().back() << "]";
return o;
}
/* -------------------------------------------------------------------------- */
} // namespace tamaas
diff --git a/src/model/model.hh b/src/model/model.hh
index 6044801..cbfdd1c 100644
--- a/src/model/model.hh
+++ b/src/model/model.hh
@@ -1,230 +1,241 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 MODEL_HH
#define MODEL_HH
/* -------------------------------------------------------------------------- */
#include "be_engine.hh"
#include "grid_base.hh"
#include "integral_operator.hh"
#include "model_dumper.hh"
#include "model_type.hh"
#include "tamaas.hh"
#include
#include
#include
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
/**
* @brief Class for linear isotropic elasticity
*/
namespace influence {
template
struct ElasticHelper {
ElasticHelper(Real mu, Real nu)
: mu(mu), nu(nu), lambda(2 * mu * nu / (1 - 2 * nu)) {}
template
Matrix, dim, dim>
operator()(const StaticMatrix& gradu) const {
auto trace = gradu.trace();
Matrix, dim, dim> sigma;
for (UInt i = 0; i < dim; ++i)
for (UInt j = 0; j < dim; ++j)
sigma(i, j) =
(i == j) * lambda * trace + mu * (gradu(i, j) + gradu(j, i));
return sigma;
}
template
SymMatrix, dim>
operator()(const StaticSymMatrix& eps) const {
SymMatrix, dim> sigma;
auto trace = eps.trace();
for (UInt i = 0; i < dim; ++i)
sigma(i) = lambda * trace + 2 * mu * eps(i);
for (UInt i = dim; i < voigt_size::value; ++i)
sigma(i) = 2 * mu * eps(i);
return sigma;
}
template
Matrix, dim, dim>
inverse(const StaticMatrix& sigma) const {
Matrix, dim, dim> epsilon;
auto trace = sigma.trace();
auto c1 = 1. / (2 * mu);
auto c2 = -lambda / (3 * lambda + 2 * mu);
for (UInt i = 0; i < dim; ++i)
for (UInt j = 0; j < dim; ++j)
epsilon(i, j) = c1 * (sigma(i, j) + c2 * trace * (i == j));
}
const Real mu, nu, lambda;
};
} // namespace influence
/* -------------------------------------------------------------------------- */
/**
* @brief Model containing pressure and displacement
* This class is a container for the model fields. It is supposed to be
* dimension agnostic, hence the GridBase members.
*/
class Model {
protected:
/// Constructor
Model(std::vector system_size, std::vector discretization);
public:
/// Destructor
virtual ~Model();
public:
/// Set elasticity parameters
void setElasticity(Real E, Real nu);
/// Get Hertz contact modulus
Real getHertzModulus() const;
/// Get Young's modulus
Real getYoungModulus() const { return E; }
/// Get Poisson's ratio
Real getPoissonRatio() const { return nu; }
/// Get shear modulus
Real getShearModulus() const { return E / (2 * (1 + nu)); }
/// Set Young's modulus
- void setYoungModulus(Real E_) { this->E = E_; }
+ void setYoungModulus(Real E_) {
+ if (E_ < 0)
+ TAMAAS_EXCEPTION("Elastic modulus should be positive");
+ this->E = E_;
+ updateOperators();
+ }
+
/// Set Poisson's ratio
- void setPoissonRatio(Real nu_) { this->nu = nu_; }
+ void setPoissonRatio(Real nu_) {
+ if (nu_ > 0.5 or nu_ <= -1)
+ TAMAAS_EXCEPTION("Poisson ratio should be in ]-1, 0.5]");
+ this->nu = nu_;
+ updateOperators();
+ }
public:
virtual void applyElasticity(GridBase& stress,
const GridBase& strain) const = 0;
public:
/// Get pressure
GridBase& getTraction();
/// Get pressure
const GridBase& getTraction() const;
/// Get displacement
GridBase& getDisplacement();
/// Get displacement
const GridBase& getDisplacement() const;
/// Get model type
virtual model_type getType() const = 0;
/// Get system physical size
const std::vector& getSystemSize() const;
/// Get boundary system physical size
virtual std::vector getBoundarySystemSize() const = 0;
/// Get discretization
const std::vector& getDiscretization() const;
/// Get boundary discretization
virtual std::vector getBoundaryDiscretization() const = 0;
/// Get boundary element engine
BEEngine& getBEEngine() {
TAMAAS_ASSERT(engine, "BEEngine was not initialized");
return *engine;
}
/// Solve Neumann problem using default neumann operator
void solveNeumann();
/// Solve Dirichlet problem using default dirichlet operator
void solveDirichlet();
public:
/// Register a new integral operator
template
IntegralOperator* registerIntegralOperator(const std::string& name) {
operators[name] = std::make_unique(this);
return operators[name].get();
}
/// Get a registerd integral operator
IntegralOperator* getIntegralOperator(const std::string& name);
/// Tell operators to update their cache
void updateOperators();
public:
/// Register a field
void registerField(const std::string& name,
std::shared_ptr> field);
/// Get a field
const GridBase& getField(const std::string& name) const;
/// Get a non-const field
GridBase& getField(const std::string& name);
/// Get fields
std::vector getFields() const;
/// Get fields map
const auto& getFieldsMap() const { return fields; }
public:
/// Set the dumper object
void addDumper(std::shared_ptr dumper);
/// Dump the model
void dump() const;
protected:
Real E = 1, nu = 0;
std::vector system_size;
std::vector discretization;
std::unique_ptr engine = nullptr;
std::unordered_map> operators;
std::unordered_map>> fields;
std::vector> dumpers;
};
/* --------------------------------------------------------------------------
*/
/* Output model to stream */
/* --------------------------------------------------------------------------
*/
std::ostream& operator<<(std::ostream& o, const Model& _this);
/* --------------------------------------------------------------------------
*/
/* Simpler grid allocation */
/* --------------------------------------------------------------------------
*/
template
std::unique_ptr> allocateGrid(Model& model) {
return allocateGrid(
model.getType(), (boundary) ? model.getBoundaryDiscretization()
: model.getDiscretization());
}
template
std::unique_ptr> allocateGrid(Model& model, UInt nc) {
return allocateGrid(model.getType(),
(boundary)
? model.getBoundaryDiscretization()
: model.getDiscretization(),
nc);
}
} // namespace tamaas
#endif // MODEL_HH