diff --git a/python/cast.hh b/python/cast.hh
index 05562597..2975a4d4 100644
--- a/python/cast.hh
+++ b/python/cast.hh
@@ -1,208 +1,208 @@
/**
* @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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __CAST_HH__
#define __CAST_HH__
/* -------------------------------------------------------------------------- */
#include "grid_base.hh"
#include "surface.hh"
#include "numpy.hh"
/* -------------------------------------------------------------------------- */
#include <pybind11/cast.h>
#include <boost/preprocessor/seq.hpp>
/* -------------------------------------------------------------------------- */
namespace pybind11 {
// Format descriptor necessary for correct wrap of tamaas complex type
template <typename T>
struct format_descriptor<
tamaas::complex<T>, detail::enable_if_t<std::is_floating_point<T>::value>> {
static constexpr const char c = format_descriptor<T>::c;
static constexpr const char value[3] = {'Z', c, '\0'};
static std::string format() { return std::string(value); }
};
#ifndef PYBIND11_CPP17
template <typename T>
constexpr const char format_descriptor<
tamaas::complex<T>,
detail::enable_if_t<std::is_floating_point<T>::value>>::value[3];
#endif
namespace detail {
// declare tamaas complex as a complex type for pybind11
template <typename T>
struct is_complex<tamaas::complex<T>> : std::true_type {};
template <typename T>
struct is_fmt_numeric<tamaas::complex<T>,
- detail::enable_if_t<std::is_floating_point<T>::value>> {
- static constexpr bool value = true;
+ detail::enable_if_t<std::is_floating_point<T>::value>>
+ : std::true_type {
static constexpr int index = is_fmt_numeric<T>::index + 3;
};
/**
* Type caster for grid classes
* inspired by https://tinyurl.com/y8m47qh3 from T. De Geus
* and pybind11/eigen.h
*/
template <template <typename, unsigned int> class G, typename T,
unsigned int dim>
struct type_caster<G<T, dim>> {
using type = G<T, dim>;
template <typename U>
using array_type = array_t<U, array::c_style | array::forcecast>;
public:
PYBIND11_TYPE_CASTER(type, _("GridWrap<T, dim>"));
/**
* Conversion part 1 (Python->C++): convert a PyObject into a grid
* instance or return false upon failure. The second argument
* indicates whether implicit conversions should be applied.
*/
bool load(handle src, bool convert) {
if (!convert && !array_type<typename type::value_type>::check_(src))
return false;
auto buf = array_type<typename type::value_type>::ensure(src);
value.move(tamaas::wrap::GridNumpy<G<T, dim>>(buf));
return true;
}
/**
* Conversion part 2 (C++ -> Python): convert a grid instance into
* a Python object. The second and third arguments are used to
* indicate the return value policy and parent object (for
* ``return_value_policy::reference_internal``) and are generally
* ignored by implicit casters.
*
* TODO: do not ignore policy (see pybind11/eigen.h)
*/
static handle cast(const type& src, return_value_policy /* policy */,
handle /*parent*/) {
if (src.getNbComponents() == 1) {
// none() passed as parent to get a correct nocopy
auto a = array_type<typename type::value_type>(
src.sizes(), src.getInternalData(), none());
return a.release();
} else {
std::array<unsigned int, dim+1> sizes;
std::copy_n(src.sizes().begin(), dim, sizes.begin());
sizes.back() = src.getNbComponents();
// none() passed as parent to get a correct nocopy
auto a = array_type<typename type::value_type>(
sizes, src.getInternalData(), none());
return a.release();
}
}
};
/**
* Type caster for GridBase classes
*/
template <typename T>
struct type_caster<tamaas::GridBase<T>> {
using type = tamaas::GridBase<T>;
template <typename U>
using array_type = array_t<U, array::c_style | array::forcecast>;
public:
PYBIND11_TYPE_CASTER(type, _("GridBaseWrap<T>"));
bool load(handle src, bool convert) {
if (!convert && !array_type<T>::check_(src))
return false;
auto buf = array_type<T>::ensure(src);
value.move(tamaas::wrap::GridBaseNumpy<T>(buf));
return true;
}
static handle cast(const type& src, return_value_policy /* policy */,
handle /*parent*/) {
#define GRID_BASE_CASE(unused1, unused2, dim) \
case dim: { \
const tamaas::Grid<T, dim>& conv = static_cast<decltype(conv)>(src); \
if (src.getNbComponents() != 1) { \
std::array<unsigned int, dim + 1> sizes; \
std::copy_n(conv.sizes().begin(), dim, sizes.begin()); \
sizes.back() = src.getNbComponents(); \
return array_type<T>(sizes, src.getInternalData(), none()).release(); \
} else { \
return array_type<T>(conv.sizes(), src.getInternalData(), none()) \
.release(); \
} \
}
switch (src.getDimension()) {
BOOST_PP_SEQ_FOR_EACH(GRID_BASE_CASE, ~, (1)(2)(3));
default:
return nullptr;
}
#undef GRID_BASE_CASE
}
};
/**
* Type caster for surface class
*/
template <typename T>
struct type_caster<tamaas::Surface<T>> {
using type = tamaas::Surface<T>;
template <typename U>
using array_type = array_t<U, array::c_style | array::forcecast>;
public:
PYBIND11_TYPE_CASTER(type, _("SurfaceWrap<T, dim>"));
bool load(handle src, bool convert) {
if (!convert && !array_type<T>::check_(src))
return false;
auto buf = array_type<T>::ensure(src);
value.move(tamaas::wrap::SurfaceNumpy<type>(buf));
return true;
}
static handle cast(const type& src, return_value_policy /* policy */,
handle /*parent*/) {
auto a = array_type<T>(src.sizes(), src.getInternalData(), none());
return a.release();
}
};
} // namespace detail
} // namespace pybind11
#endif // __CAST_HH__
diff --git a/python/wrap/model.cpp b/python/wrap/model.cpp
index 40890e58..a7269f4f 100644
--- a/python/wrap/model.cpp
+++ b/python/wrap/model.cpp
@@ -1,150 +1,147 @@
/**
* @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 "wrap.hh"
#include "model.hh"
#include "model_factory.hh"
#include "functional.hh"
#include "adhesion_functional.hh"
#include <pybind11/stl.h>
/* -------------------------------------------------------------------------- */
__BEGIN_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<Real>& variable,
GridBase<Real>& dual) const override {
PYBIND11_OVERLOAD_PURE(Real, Functional, computeF, variable, dual);
}
void computeGradF(GridBase<Real>& variable,
GridBase<Real>& gradient) const override {
PYBIND11_OVERLOAD_PURE(void, Functional, computeGradF, variable, gradient);
}
};
}
}
namespace wrap {
using namespace py::literals;
/// Wrap functional classes
void wrapFunctionals(py::module& mod) {
py::class_<functional::Functional, functional::wrap::PyFunctional> func(
mod, "Functional");
- func.def(py::init<const BEEngine&>())
+ func.def(py::init<>())
.def("computeF", &functional::Functional::computeF)
.def("computeGradF", &functional::Functional::computeGradF);
py::class_<functional::AdhesionFunctional> adh(mod, "AdhesionFunctional",
func);
adh.def_property("parameters", &functional::AdhesionFunctional::getParameters,
&functional::AdhesionFunctional::setParameters)
// legacy wrapper code
.def("setParameters", &functional::AdhesionFunctional::setParameters);
py::class_<functional::ExponentialAdhesionFunctional>(
mod, "ExponentialAdhesionFunctional", adh)
- .def(py::init<const BEEngine&, const GridBase<Real>&>(), "be_engine"_a,
- "surface"_a);
+ .def(py::init<const GridBase<Real>&>(), "surface"_a);
py::class_<functional::MaugisAdhesionFunctional>(
mod, "MaugisAdhesionFunctional", adh)
- .def(py::init<const BEEngine&, const GridBase<Real>&>(), "be_engine"_a,
- "surface"_a);
+ .def(py::init<const GridBase<Real>&>(), "surface"_a);
py::class_<functional::SquaredExponentialAdhesionFunctional>(
mod, "SquaredExponentialAdhesionFunctional", adh)
- .def(py::init<const BEEngine&, const GridBase<Real>&>(), "be_engine"_a,
- "surface"_a);
+ .def(py::init<const GridBase<Real>&>(), "surface"_a);
}
/// Wrap BEEngine classes
void wrapBEEngine(py::module& mod) {
py::class_<BEEngine>(mod, "BEEngine")
.def("solveNeumann", &BEEngine::solveNeumann)
.def("solveDirichlet", &BEEngine::solveDirichlet)
.def("getModel", &BEEngine::getModel, py::return_value_policy::reference);
}
/// Wrap Models
void wrapModelClass(py::module& mod) {
py::enum_<model_type>(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_<Model>(mod, "Model")
.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<Real> & (Model::*)()) & Model::getTraction)
.def("getDisplacement",
(GridBase<Real> & (Model::*)()) & Model::getDisplacement)
.def("getSystemSize", &Model::getSystemSize)
.def("getDiscretization", &Model::getDiscretization)
.def("solveNeumann", &Model::solveNeumann)
.def("solveDirichlet", &Model::solveDirichlet)
.def("getBEEngine", &Model::getBEEngine,
py::return_value_policy::reference)
.def("__repr__", [](const Model& m) {
std::stringstream ss;
ss << m;
return ss.str();
});
}
/// Wrap factor for models
void wrapModelFactory(py::module& mod) {
py::class_<ModelFactory>(mod, "ModelFactory")
.def_static("createModel", &ModelFactory::createModel, "model_type"_a,
"system_size"_a, "discretization"_a);
}
void wrapModel(py::module& mod) {
wrapBEEngine(mod);
wrapModelClass(mod);
wrapModelFactory(mod);
wrapFunctionals(mod);
}
} // namespace wrap
__END_TAMAAS__
diff --git a/src/model/adhesion_functional.hh b/src/model/adhesion_functional.hh
index 88cf66a4..7aad8d6b 100644
--- a/src/model/adhesion_functional.hh
+++ b/src/model/adhesion_functional.hh
@@ -1,110 +1,107 @@
/**
* @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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __ADHESION_FUNCTIONAL_HH__
#define __ADHESION_FUNCTIONAL_HH__
/* -------------------------------------------------------------------------- */
#include "functional.hh"
#include <map>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
namespace functional {
/// Functional class for adhesion energy
class AdhesionFunctional : public Functional {
protected:
/// Constructor
- AdhesionFunctional(const BEEngine& engine, const GridBase<Real>& surface)
- : Functional(engine), surface() {
+ AdhesionFunctional(const GridBase<Real>& surface)
+ : Functional(), surface() {
this->surface.wrap(surface);
}
public:
/// Get parameters
const std::map<std::string, Real>& getParameters() const {
return parameters;
}
/// Set parameters
void setParameters(std::map<std::string, Real> other) {
parameters = other;
}
protected:
GridBase<Real> surface;
std::map<std::string, Real> parameters;
};
/// Exponential adhesion functional
class ExponentialAdhesionFunctional : public AdhesionFunctional {
public:
/// Explicit declaration of constructor for swig
- ExponentialAdhesionFunctional(const BEEngine& engine,
- const GridBase<Real>& surface)
- : AdhesionFunctional(engine, surface) {}
+ ExponentialAdhesionFunctional(const GridBase<Real>& surface)
+ : AdhesionFunctional(surface) {}
/// Compute the total adhesion energy
Real computeF(GridBase<Real>& gap, GridBase<Real>& pressure) const override;
/// Compute the gradient of the adhesion functional
void computeGradF(GridBase<Real>& gap,
GridBase<Real>& gradient) const override;
};
/// Constant adhesion functional
class MaugisAdhesionFunctional : public AdhesionFunctional {
public:
/// Explicit declaration of constructor for swig
- MaugisAdhesionFunctional(const BEEngine& engine,
- const GridBase<Real>& surface)
- : AdhesionFunctional(engine, surface) {}
+ MaugisAdhesionFunctional(const GridBase<Real>& surface)
+ : AdhesionFunctional(surface) {}
/// Compute the total adhesion energy
Real computeF(GridBase<Real>& gap, GridBase<Real>& pressure) const override;
/// Compute the gradient of the adhesion functional
void computeGradF(GridBase<Real>& gap,
GridBase<Real>& gradient) const override;
};
/// Squared exponential adhesion functional
class SquaredExponentialAdhesionFunctional : public AdhesionFunctional {
public:
/// Explicit declaration of constructor for swig
- SquaredExponentialAdhesionFunctional(const BEEngine& engine,
- const GridBase<Real>& surface)
- : AdhesionFunctional(engine, surface) {}
+ SquaredExponentialAdhesionFunctional(const GridBase<Real>& surface)
+ : AdhesionFunctional(surface) {}
/// Compute the total adhesion energy
Real computeF(GridBase<Real>& gap, GridBase<Real>& pressure) const override;
/// Compute the gradient of the adhesion functional
void computeGradF(GridBase<Real>& gap,
GridBase<Real>& gradient) const override;
};
} // namespace functional
__END_TAMAAS__
#endif // __ADHESION_FUNCTIONAL_HH__
diff --git a/src/model/elastic_functional.cpp b/src/model/elastic_functional.cpp
index f79ab853..e88a4e34 100644
--- a/src/model/elastic_functional.cpp
+++ b/src/model/elastic_functional.cpp
@@ -1,77 +1,77 @@
/**
* @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 "elastic_functional.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
namespace functional {
Real ElasticFunctionalPressure::computeF(GridBase<Real>& pressure,
GridBase<Real>& dual) const {
// This doesn't work if dual is multi-component
*buffer = dual;
// this can get large
unsigned long long n = dual.getNbPoints();
n *= n;
return (0.5 * pressure.dot(*buffer) - pressure.dot(surface)) / n;
}
/* -------------------------------------------------------------------------- */
void ElasticFunctionalPressure::computeGradF(GridBase<Real>& pressure,
GridBase<Real>& gradient) const {
- this->engine.solveNeumann(pressure, *buffer);
+ this->op.apply(pressure, *buffer);
// This doesn't work if dual is multi-component
*buffer -= surface;
gradient += *buffer;
}
/* -------------------------------------------------------------------------- */
Real ElasticFunctionalGap::computeF(GridBase<Real>& gap,
GridBase<Real>& dual) const {
// This doesn't work if dual is multi-component
*buffer = gap;
*buffer += surface;
UInt n = dual.getNbPoints();
Real F = 0.5 * buffer->dot(dual);
return F / n / n; // avoid integer overflow
}
/* -------------------------------------------------------------------------- */
void ElasticFunctionalGap::computeGradF(GridBase<Real>& gap,
GridBase<Real>& gradient) const {
*buffer = gap;
*buffer += surface; // Now we have displacements
- this->engine.solveDirichlet(*buffer, *buffer);
+ this->op.apply(*buffer, *buffer);
gradient += *buffer;
}
} // namespace functional
__END_TAMAAS__
diff --git a/src/model/elastic_functional.hh b/src/model/elastic_functional.hh
index 42fc94b6..04484e3b 100644
--- a/src/model/elastic_functional.hh
+++ b/src/model/elastic_functional.hh
@@ -1,84 +1,84 @@
/**
* @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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __ELASTIC_FUNCTIONAL_HH__
#define __ELASTIC_FUNCTIONAL_HH__
/* -------------------------------------------------------------------------- */
#include "functional.hh"
#include "model.hh"
#include "tamaas.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
namespace functional {
/// Generic functional for elastic energy
class ElasticFunctional : public Functional {
public:
- ElasticFunctional(BEEngine& engine, const GridBase<Real>& surface)
- : Functional(engine), surface(surface) {
- const Model& model = engine.getModel();
- buffer =
- allocateGrid<true, Real>(model.getType(), model.getDiscretization(),
- model.getTraction().getNbComponents());
+ ElasticFunctional(const IntegralOperator& op, const GridBase<Real>& surface)
+ : Functional(), op(op), surface(surface) {
+ const Model& model = op.getModel();
+ buffer = allocateGrid<true, Real>(op.getType(),
+ model.getBoundaryDiscretization());
}
protected:
+ const IntegralOperator& op;
const GridBase<Real>& surface;
std::unique_ptr<GridBase<Real>> mutable buffer;
};
/* -------------------------------------------------------------------------- */
/// Functional with pressure as primal field
class ElasticFunctionalPressure : public ElasticFunctional {
public:
using ElasticFunctional::ElasticFunctional;
/// Compute functional with input pressure
Real computeF(GridBase<Real>& pressure, GridBase<Real>& dual) const override;
/// Compute functional gradient with input pressure
void computeGradF(GridBase<Real>& pressure,
GridBase<Real>& gradient) const override;
};
/* -------------------------------------------------------------------------- */
/// Functional with gap as primal field
class ElasticFunctionalGap : public ElasticFunctional {
public:
using ElasticFunctional::ElasticFunctional;
/// Compute functional with input gap
Real computeF(GridBase<Real>& gap, GridBase<Real>& dual) const override;
/// Compute functional gradient with input gap
void computeGradF(GridBase<Real>& gap,
GridBase<Real>& gradient) const override;
};
} // namespace functional
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
#endif // __ELASTIC_FUNCTIONAL_HH__
diff --git a/src/model/functional.hh b/src/model/functional.hh
index 0b48d1ef..0250689b 100644
--- a/src/model/functional.hh
+++ b/src/model/functional.hh
@@ -1,61 +1,56 @@
/**
* @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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __FUNCTIONAL_HH__
#define __FUNCTIONAL_HH__
/* -------------------------------------------------------------------------- */
#include "be_engine.hh"
#include "tamaas.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
namespace functional {
/// Generic functional class for the cost function of the optimization problem
class Functional {
public:
- /// Constructor
- explicit Functional(const BEEngine& engine) : engine(engine) {}
/// Destructor
virtual ~Functional() = default;
public:
/// Compute functional
virtual Real computeF(GridBase<Real>& variable,
GridBase<Real>& dual) const = 0;
/// Compute functional gradient
virtual void computeGradF(GridBase<Real>& variable,
GridBase<Real>& gradient) const = 0;
-
-protected:
- const BEEngine& engine;
};
} // namespace functional
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
#endif // ___FUNCTIONAL_HH__
diff --git a/src/model/integral_operator.hh b/src/model/integral_operator.hh
index 4a2a9c95..b872d88f 100644
--- a/src/model/integral_operator.hh
+++ b/src/model/integral_operator.hh
@@ -1,91 +1,98 @@
/**
* @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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef INTEGRAL_OPERATOR_HH
#define INTEGRAL_OPERATOR_HH
/* -------------------------------------------------------------------------- */
#include "grid_base.hh"
+#include "model_type.hh"
#include <boost/preprocessor/seq.hpp>
#include <boost/preprocessor/cat.hpp>
#include <boost/preprocessor/stringize.hpp>
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
class Model;
/**
* @brief Generic class for integral operators
*/
class IntegralOperator {
public:
/// Kind of operator
enum kind { neumann, dirichlet };
public:
/// Constructor
IntegralOperator(Model* model) : model(model) {}
/// Virtual destructor for subclasses
virtual ~IntegralOperator() = default;
/// Apply operator on input
virtual void apply(GridBase<Real>& input, GridBase<Real>& output) const = 0;
/// Update any data dependent on model parameters
virtual void updateFromModel() = 0;
/// Get model
const Model& getModel() const { return *model; }
+ /// Kind
+ virtual kind getKind() const = 0;
+
+ /// Type
+ virtual model_type getType() const = 0;
+
protected:
Model* model = nullptr;
};
/* -------------------------------------------------------------------------- */
/* Printing IntegralOperator::kind to string */
/* -------------------------------------------------------------------------- */
#define INTEGRAL_KINDS (neumann)(dirichlet)
inline std::ostream& operator<<(std::ostream& o,
const IntegralOperator::kind& val) {
switch (val) {
#define PRINT_INTEGRAL_KIND(r, data, kind) \
case data::kind: \
o << BOOST_PP_STRINGIZE(kind); \
break;
BOOST_PP_SEQ_FOR_EACH(PRINT_INTEGRAL_KIND, IntegralOperator,
INTEGRAL_KINDS);
#undef PRINT_INTEGRAL_KIND
}
return o;
}
#undef INTEGRAL_KINDS
/* -------------------------------------------------------------------------- */
} // namespace tamaas
#endif // INTEGRAL_OPERATOR_HH
diff --git a/src/model/meta_functional.hh b/src/model/meta_functional.hh
index 65dd5aad..57543457 100644
--- a/src/model/meta_functional.hh
+++ b/src/model/meta_functional.hh
@@ -1,70 +1,68 @@
/**
* @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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __META_FUNCTIONAL_HH__
#define __META_FUNCTIONAL_HH__
#include "functional.hh"
#include "tamaas.hh"
#include <list>
#include <memory>
__BEGIN_TAMAAS__
namespace functional {
/// Meta functional that contains list of functionals
class MetaFunctional : public Functional {
using FunctionalList = std::list<std::pair<Functional*, bool>>;
public:
- using Functional::Functional;
-
~MetaFunctional() override;
public:
/// Compute functional
Real computeF(GridBase<Real>& variable, GridBase<Real>& dual) const override;
/// Compute functional gradient
void computeGradF(GridBase<Real>& variable,
GridBase<Real>& gradient) const override;
/// Add functional to the list
void addFunctionalTerm(Functional* term, bool own);
/// Add functional to the list
void addFunctionalTerm(std::pair<Functional*, bool> element);
protected:
FunctionalList functionals;
};
} // namespace functional
__END_TAMAAS__
#endif // __META_FUNCTIONAL_HH__
diff --git a/src/model/model.hh b/src/model/model.hh
index 9abfcb36..b9e8db4c 100644
--- a/src/model/model.hh
+++ b/src/model/model.hh
@@ -1,135 +1,135 @@
/**
* @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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef MODEL_HH
#define MODEL_HH
/* -------------------------------------------------------------------------- */
#include "grid_base.hh"
#include "model_type.hh"
#include "integral_operator.hh"
#include "be_engine.hh"
#include "tamaas.hh"
#include <map>
#include <memory>
#include <vector>
#include <algorithm>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/**
* @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(const std::vector<Real>& system_size,
const std::vector<UInt>& 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_; }
/// Set Poisson's ratio
void setPoissonRatio(Real nu_) { this->nu = nu_; }
public:
/// Get pressure
GridBase<Real>& getTraction();
/// Get pressure
const GridBase<Real>& getTraction() const;
/// Get displacement
GridBase<Real>& getDisplacement();
/// Get displacement
const GridBase<Real>& getDisplacement() const;
- /// Get displacement on boundary
- virtual GridBase<Real>& getBoundaryDisplacement() = 0;
/// Get model type
virtual ::tamaas::model_type getType() const = 0;
/// Get system physical size
const std::vector<Real>& getSystemSize() const;
/// Get discretization
const std::vector<UInt>& getDiscretization() const;
+ /// Get boundary discretization
+ virtual std::vector<UInt> 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 <typename Operator>
std::shared_ptr<IntegralOperator>
registerIntegralOperator(const std::string& name) {
operators[name] = std::make_shared<Operator>(this);
return operators[name];
}
/// Get a registerd integral operator
std::shared_ptr<IntegralOperator>
getIntegralOperator(const std::string& name);
/// Tell operators to update their cache
void updateOperators();
protected:
Real E = 1, nu = 0;
std::vector<Real> system_size;
std::vector<UInt> discretization;
std::unique_ptr<GridBase<Real>> traction = nullptr;
std::unique_ptr<GridBase<Real>> displacement = nullptr;
std::unique_ptr<BEEngine> engine = nullptr;
std::map<std::string, std::shared_ptr<IntegralOperator>> operators;
};
/* -------------------------------------------------------------------------- */
/* Output model to stream */
/* -------------------------------------------------------------------------- */
std::ostream& operator<<(std::ostream& o, const Model& _this);
__END_TAMAAS__
#endif // MODEL_HH
diff --git a/src/model/model_factory.cpp b/src/model/model_factory.cpp
index b85dd5b0..a6702d8e 100644
--- a/src/model/model_factory.cpp
+++ b/src/model/model_factory.cpp
@@ -1,58 +1,56 @@
/**
* @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 "model_factory.hh"
#include "model_template.hh"
#include <boost/preprocessor/seq.hpp>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
#define MODEL_TYPE_CASE_MACRO(truc, unused, type) \
case type: { \
model = new ModelTemplate<type>(system_size, discretization); \
break; \
}
-#define IMPLEMENTED_MODEL_TYPES \
- (model_type::basic_2d)(model_type::basic_1d)(model_type::volume_2d)( \
- model_type::surface_2d)
+#define IMPLEMENTED_MODEL_TYPES TAMAAS_MODEL_TYPES
Model* ModelFactory::createModel(model_type type,
const std::vector<Real>& system_size,
const std::vector<UInt>& discretization) {
Model* model = nullptr;
switch (type) {
BOOST_PP_SEQ_FOR_EACH(MODEL_TYPE_CASE_MACRO, ~, IMPLEMENTED_MODEL_TYPES);
default:
TAMAAS_EXCEPTION("Model type not implemented");
}
return model;
}
__END_TAMAAS__
diff --git a/src/model/model_template.cpp b/src/model/model_template.cpp
index d4de8079..b6d40728 100644
--- a/src/model/model_template.cpp
+++ b/src/model/model_template.cpp
@@ -1,85 +1,99 @@
/**
* @file
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016-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 "model_template.hh"
#include "westergaard.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
template <model_type type>
ModelTemplate<type>::ModelTemplate(const std::vector<Real>& system_size,
const std::vector<UInt>& discretization)
: Model(system_size, discretization) {
constexpr UInt dim = trait::dimension;
constexpr UInt dim_b = trait::boundary_dimension;
constexpr UInt nb_components = trait::components;
if (system_size.size() != dim)
TAMAAS_EXCEPTION("System size does not match model type");
if (discretization.size() != dim)
TAMAAS_EXCEPTION("Discretization size does not match model type");
// Copying sizes for traction grid
std::array<UInt, dim_b> traction_size;
auto disc_it = discretization.begin();
if (dim > dim_b)
++disc_it;
std::copy(disc_it, discretization.end(), traction_size.begin());
// Allocating
traction = std::make_unique<Grid<Real, dim_b>>(traction_size, nb_components);
displacement =
std::make_unique<Grid<Real, dim>>(discretization, nb_components);
- boundary_displacement =
- std::make_unique<ViewType>(*displacement, trait::indices);
this->initializeBEEngine();
}
/* -------------------------------------------------------------------------- */
template <model_type type>
void ModelTemplate<type>::initializeBEEngine() {
engine = std::make_unique<BEEngineTmpl<type>>(this);
}
-/* -------------------------------------------------------------------------- */
-
template <model_type type>
-GridBase<Real>& ModelTemplate<type>::getBoundaryDisplacement() {
- return *boundary_displacement;
+std::vector<UInt> ModelTemplate<type>::getBoundaryDiscretization() const {
+ return getDiscretization();
+}
+
+template <>
+std::vector<UInt>
+ModelTemplate<model_type::volume_1d>::getBoundaryDiscretization() const {
+ auto begin = getDiscretization().begin() + 1, end = getDiscretization().end();
+ std::vector<UInt> size(getDiscretization().size() - 1);
+ std::copy(begin, end, size.begin());
+ return size;
+}
+
+template <>
+std::vector<UInt>
+ModelTemplate<model_type::volume_2d>::getBoundaryDiscretization() const {
+ auto begin = getDiscretization().begin() + 1, end = getDiscretization().end();
+ std::vector<UInt> size(getDiscretization().size() - 1);
+ std::copy(begin, end, size.begin());
+ return size;
}
/* -------------------------------------------------------------------------- */
/* Template instanciation */
/* -------------------------------------------------------------------------- */
template class ModelTemplate<model_type::basic_1d>;
template class ModelTemplate<model_type::basic_2d>;
template class ModelTemplate<model_type::surface_1d>;
template class ModelTemplate<model_type::surface_2d>;
template class ModelTemplate<model_type::volume_1d>;
template class ModelTemplate<model_type::volume_2d>;
__END_TAMAAS__
diff --git a/src/model/model_template.hh b/src/model/model_template.hh
index bf59eb47..51fcfbbf 100644
--- a/src/model/model_template.hh
+++ b/src/model/model_template.hh
@@ -1,64 +1,65 @@
/**
* @file
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016-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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __MODEL_TEMPLATE_HH__
#define __MODEL_TEMPLATE_HH__
/* -------------------------------------------------------------------------- */
#include "grid_view.hh"
#include "model.hh"
#include "model_type.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/**
* @brief Model class templated with model type
* Specializations of this class should take care of dimension specific code
*/
template <model_type type>
class ModelTemplate : public Model {
using trait = model_type_traits<type>;
using ViewType =
GridView<Grid, Real, trait::dimension, trait::boundary_dimension>;
public:
/// Constructor
ModelTemplate(const std::vector<Real>& system_size,
const std::vector<UInt>& discretization);
- GridBase<Real>& getBoundaryDisplacement() override;
-
// Get model type
model_type getType() const override { return type; }
+ std::vector<UInt> getBoundaryDiscretization() const override;
+
protected:
virtual void initializeBEEngine();
protected:
- std::unique_ptr<ViewType> boundary_displacement = nullptr;
+ std::unique_ptr<ViewType> displacement_view = nullptr;
+ std::unique_ptr<ViewType> traction_view = nullptr;
};
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
#endif // __MODEL_TEMPLATE_HH__
diff --git a/src/model/model_type.hh b/src/model/model_type.hh
index 50fc3685..23aff783 100644
--- a/src/model/model_type.hh
+++ b/src/model/model_type.hh
@@ -1,122 +1,140 @@
/**
* @file
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016-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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __MODEL_TYPE_HH__
#define __MODEL_TYPE_HH__
/* -------------------------------------------------------------------------- */
#include "grid.hh"
#include "grid_base.hh"
#include "tamaas.hh"
#include <boost/preprocessor/seq.hpp>
#include <boost/preprocessor/cat.hpp>
#include <boost/preprocessor/stringize.hpp>
#include <memory>
+#include <algorithm>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
// clang-format off
+/// Models types **for printing only**
#define TAMAAS_MODEL_TYPES \
(basic_1d) \
(basic_2d) \
(surface_1d) \
(surface_2d) \
(volume_1d) \
(volume_2d)
// clang-format on
/// Types for grid dimensions and number of components
enum class model_type {
basic_1d, ///< one component line
basic_2d, ///< one component surface
surface_1d, ///< two components line
surface_2d, ///< three components surface
volume_1d, ///< two components volume
volume_2d ///< three components volume
};
inline std::ostream& operator<<(std::ostream& o, const model_type& val) {
switch (val) {
#define PRINT_MODEL_TYPE(r, data, type) \
case data::type: \
o << BOOST_PP_STRINGIZE(type); \
break;
BOOST_PP_SEQ_FOR_EACH(PRINT_MODEL_TYPE, model_type, TAMAAS_MODEL_TYPES);
#undef PRINT_MODEL_TYPE
}
return o;
}
/// Trait class to store physical dimension of domain, of boundary and
/// number of components
template <model_type type>
struct model_type_traits {};
#define MODEL_TYPE_TRAITS_MACRO(type, dim, comp, bdim) \
template <> \
struct model_type_traits<model_type::type> { \
static constexpr UInt dimension = dim; \
static constexpr UInt components = comp; \
static constexpr UInt boundary_dimension = bdim; \
static std::vector<UInt> indices; \
}
MODEL_TYPE_TRAITS_MACRO(basic_1d, 1, 1, 1);
MODEL_TYPE_TRAITS_MACRO(basic_2d, 2, 1, 2);
MODEL_TYPE_TRAITS_MACRO(surface_1d, 1, 2, 1);
MODEL_TYPE_TRAITS_MACRO(surface_2d, 2, 3, 2);
MODEL_TYPE_TRAITS_MACRO(volume_1d, 2, 2, 1);
MODEL_TYPE_TRAITS_MACRO(volume_2d, 3, 3, 2);
#undef MODEL_TYPE_TRAITS_MACRO
#define ALLOC_GRID_CASE_MACRO(a, data, type) \
case data::type: { \
- constexpr UInt dim = \
- (boundary) ? model_type_traits<data::type>::boundary_dimension \
- : model_type_traits<data::type>::dimension; \
- ptr = new Grid<T, dim>(n, nc); \
+ constexpr UInt bdim = model_type_traits<data::type>::boundary_dimension; \
+ constexpr UInt dim = model_type_traits<data::type>::dimension; \
+ ptr = new Grid < T, \
+ (boundary) ? bdim \
+ : dim > (std::begin(n), std::end(n), \
+ model_type_traits<data::type>::components); \
break; \
}
/// Helper function for grid allocation
-template <bool boundary, typename T>
+template <bool boundary, typename T, typename Container>
std::unique_ptr<GridBase<T>> allocateGrid(model_type type,
- const std::vector<UInt> n, UInt nc) {
+ Container&& n) {
GridBase<T>* ptr = nullptr;
+
switch (type) {
BOOST_PP_SEQ_FOR_EACH(ALLOC_GRID_CASE_MACRO, model_type,
TAMAAS_MODEL_TYPES);
}
return std::unique_ptr<GridBase<T>>{ptr};
}
#undef ALLOC_GRID_CASE_MACRO
+// clang-format off
+#ifdef TAMAAS_MODEL_TYPES
+#undef TAMAAS_MODEL_TYPES
+#endif
+#define TAMAAS_MODEL_TYPES \
+ (model_type::basic_1d) \
+ (model_type::basic_2d) \
+ (model_type::surface_1d) \
+ (model_type::surface_2d) \
+ (model_type::volume_1d) \
+ (model_type::volume_2d)
+// clang-format on
+
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
#endif // __MODEL_TYPE_HH__
diff --git a/src/model/volume_potential.hh b/src/model/volume_potential.hh
index a1f0c4ee..6c72af0b 100644
--- a/src/model/volume_potential.hh
+++ b/src/model/volume_potential.hh
@@ -1,97 +1,105 @@
/**
* @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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef VOLUME_POTENTIAL_HH
#define VOLUME_POTENTIAL_HH
/* -------------------------------------------------------------------------- */
#include "integral_operator.hh"
#include "model_type.hh"
#include "grid_view.hh"
#include "grid_hermitian.hh"
#include "fft_plan_manager.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/**
* @brief Volume potential operator class. Applies the operators for computation
* of displacements and strains due to residual/eigen strains
*/
template <model_type type>
class VolumePotential : public IntegralOperator {
using trait = model_type_traits<type>;
public:
VolumePotential(Model * model);
/// Update from model (does nothing)
void updateFromModel() override {}
+ /// Kind
+ IntegralOperator::kind getKind() const override {
+ return IntegralOperator::neumann;
+ }
+
+ /// Type
+ model_type getType() const override { return this->model->getType(); }
+
protected:
/// Function to handle layer-by-layer Fourier treatment of operators
template <typename Func>
void fourierApply(Func func, GridBase<Real>& in, GridBase<Real>& out) const;
protected:
Grid<Real, trait::boundary_dimension> wavevectors;
using BufferType = GridHermitian<Real, trait::boundary_dimension>;
mutable std::vector<BufferType> source_buffers;
mutable BufferType disp_buffer;
};
/* -------------------------------------------------------------------------- */
/* Template implementation */
/* -------------------------------------------------------------------------- */
template <model_type type>
template <typename Func>
void VolumePotential<type>::fourierApply(Func func, GridBase<Real>& in,
GridBase<Real>& out) const {
constexpr UInt dim = trait::dimension;
Grid<Real, dim>& i = dynamic_cast<decltype(i)>(in);
Grid<Real, dim>& o = dynamic_cast<decltype(o)>(out);
// TAMAAS_ASSERT(i.sizes().front() == o.sizes().front(),
// "Number of layers does not match");
// ^^^^ not applicable for Boussinesq
for (UInt layer : Loop::range(i.sizes().front())) {
auto in_layer = make_view(i, layer);
FFTPlanManager::get()
.createPlan(in_layer, source_buffers[layer])
.forwardTransform();
}
for (UInt layer : Loop::range(o.sizes().front())) {
func(layer, source_buffers, disp_buffer);
auto out_layer = make_view(o, layer);
FFTPlanManager::get()
.createPlan(out_layer, disp_buffer)
.backwardTransform();
}
}
__END_TAMAAS__
#endif // VOLUME_POTENTIAL_HH
diff --git a/src/model/westergaard.cpp b/src/model/westergaard.cpp
index 8d58549e..7e2fd877 100644
--- a/src/model/westergaard.cpp
+++ b/src/model/westergaard.cpp
@@ -1,215 +1,209 @@
/**
* @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 "fft_plan_manager.hh"
#include "influence.hh"
#include "loop.hh"
#include "model.hh"
#include "grid_view.hh"
#include "static_types.hh"
#include "westergaard.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/* -------------------------------------------------------------------------- */
template <model_type mtype, IntegralOperator::kind otype>
Westergaard<mtype, otype>::Westergaard(Model* model) : IntegralOperator(model) {
// Copy sizes
- std::array<UInt, trait::dimension> sizes;
- std::copy(model->getDiscretization().begin(),
- model->getDiscretization().end(), sizes.begin());
+ std::array<UInt, trait::boundary_dimension> sizes;
+ auto bdisc = model->getBoundaryDiscretization();
- // Copy boundary sizes
- auto sizes_it = sizes.begin();
- std::array<UInt, trait::boundary_dimension> boundary_hermitian_sizes;
// Ignoring first dimension if model is volumetric
- if (trait::dimension > trait::boundary_dimension)
- ++sizes_it;
- std::copy(sizes_it, sizes.end(), boundary_hermitian_sizes.begin());
- boundary_hermitian_sizes =
+ std::copy(bdisc.begin(),
+ bdisc.end(), sizes.begin());
+
+ auto boundary_hermitian_sizes =
GridHermitian<Real, trait::boundary_dimension>::hermitianDimensions(
- boundary_hermitian_sizes);
+ sizes);
constexpr UInt nb_components = trait::components;
buffer.setNbComponents(nb_components);
buffer.resize(boundary_hermitian_sizes);
influence.setNbComponents(nb_components * nb_components);
influence.resize(boundary_hermitian_sizes);
initInfluence();
}
/* -------------------------------------------------------------------------- */
template <model_type mtype, IntegralOperator::kind otype>
template <typename Functor>
void Westergaard<mtype, otype>::fourierApply(Functor func, GridBase<Real>& in,
GridBase<Real>& out) const try {
Grid<Real, bdim>& i = dynamic_cast<decltype(i)>(in);
Grid<Real, dim>& full_out = dynamic_cast<decltype(full_out)>(out);
FFTPlanManager::get().createPlan(i, buffer).forwardTransform();
/// Applying influence
func(buffer, influence);
/// Backward fourier transform on boundary only
if constexpr (bdim == dim)
FFTPlanManager::get().createPlan(full_out, buffer).backwardTransform();
else {
auto view = make_view(full_out, 0);
FFTPlanManager::get()
.createPlan(view, buffer)
.backwardTransform();
}
} catch (const std::bad_cast& e) {
TAMAAS_EXCEPTION("Neumann and dirichlet types do not match model type");
}
/* -------------------------------------------------------------------------- */
template <model_type mtype, IntegralOperator::kind otype>
template <typename Functor>
void Westergaard<mtype, otype>::initFromFunctor(Functor func) {
// Compute wavevectors for influence
auto wavevectors = FFTransform<Real, bdim>::template computeFrequencies<true>(
influence.sizes());
// Get boundary physical size
VectorProxy<const Real, bdim> domain(
this->model->getSystemSize()[(bdim == dim) ? 0 : 1]);
// Normalize wavevectors
wavevectors *= 2 * M_PI;
wavevectors /= domain;
// Apply functor
Loop::stridedLoop(func, wavevectors, influence);
// Set fundamental mode to zero
MatrixProxy<Complex, trait::components, trait::components> mat(influence(0));
mat = 0;
}
/* -------------------------------------------------------------------------- */
#define NEUMANN_BASIC(type) \
template <> \
void Westergaard<type, IntegralOperator::neumann>::initInfluence() { \
auto E_star = model->getHertzModulus(); \
auto basic = [E_star] CUDA_LAMBDA(VectorProxy<Real, bdim> && q, \
Complex & k) { \
k = 2. / (E_star * q.l2norm()); \
}; \
initFromFunctor(basic); \
}
#define DIRICHLET_BASIC(type) \
template <> \
void Westergaard<type, IntegralOperator::dirichlet>::initInfluence() { \
auto E_star = model->getHertzModulus(); \
auto basic = [E_star] CUDA_LAMBDA(VectorProxy<Real, bdim> && q, \
Complex & k) { \
k = E_star * q.l2norm() / 2; \
}; \
initFromFunctor(basic); \
}
NEUMANN_BASIC(model_type::basic_1d);
NEUMANN_BASIC(model_type::basic_2d);
DIRICHLET_BASIC(model_type::basic_1d);
DIRICHLET_BASIC(model_type::basic_2d);
#undef NEUMANN_BASIC
#undef DIRICHLET_BASIC
/* -------------------------------------------------------------------------- */
template <>
void Westergaard<model_type::surface_2d,
IntegralOperator::neumann>::initInfluence() {
auto E = model->getYoungModulus();
auto nu = model->getPoissonRatio();
const Complex I(0, 1);
auto surface = [E, nu, I] CUDA_LAMBDA(VectorProxy<Real, bdim> && q,
MatrixProxy<Complex, comp, comp>&& F) {
Real q_norm = q.l2norm();
q /= q_norm;
F(0, 0) = 2 * (1 + nu) * (1 - nu * q(0) * q(0));
F(1, 1) = 2 * (1 + nu) * (1 - nu * q(1) * q(1));
F(2, 2) = 2 * (1 - nu * nu);
F(0, 1) = F(1, 0) = -q(0) * q(1) * 2 * nu * (1 + nu);
F(0, 2) = I * q(0) * (1 + nu) * (1. - 2. * nu);
F(1, 2) = I * q(1) * (1 + nu) * (1 - 2 * nu);
F(2, 0) = -F(0, 2);
F(2, 1) = -F(1, 2);
F *= 1. / (E * q_norm);
};
initFromFunctor(surface);
}
/* -------------------------------------------------------------------------- */
template <model_type mtype, IntegralOperator::kind otype>
-void Westergaard<mtype, otype>::initInfluence() {
- TAMAAS_EXCEPTION("the requested operator has not been implemented");
-}
+void Westergaard<mtype, otype>::initInfluence() {}
/* ------------------------------------------------------------------------ */
template <model_type mtype, IntegralOperator::kind otype>
void Westergaard<mtype, otype>::apply(GridBase<Real>& input,
GridBase<Real>& output) const {
auto apply = [](decltype(buffer)& buffer,
const decltype(influence)& influence) {
Loop::stridedLoop([] CUDA_LAMBDA(VectorProxy<Complex, comp> && i,
MatrixProxy<const Complex, comp, comp> &&
F) { i = F * i; },
buffer, influence);
};
fourierApply(apply, input, output);
}
/* -------------------------------------------------------------------------- */
/* Template instanciation */
/* -------------------------------------------------------------------------- */
template class Westergaard<model_type::basic_1d, IntegralOperator::neumann>;
template class Westergaard<model_type::basic_2d, IntegralOperator::neumann>;
template class Westergaard<model_type::basic_1d, IntegralOperator::dirichlet>;
template class Westergaard<model_type::basic_2d, IntegralOperator::dirichlet>;
template class Westergaard<model_type::surface_1d, IntegralOperator::neumann>;
template class Westergaard<model_type::surface_2d, IntegralOperator::neumann>;
template class Westergaard<model_type::surface_1d, IntegralOperator::dirichlet>;
template class Westergaard<model_type::surface_2d, IntegralOperator::dirichlet>;
template class Westergaard<model_type::volume_1d, IntegralOperator::neumann>;
template class Westergaard<model_type::volume_2d, IntegralOperator::neumann>;
template class Westergaard<model_type::volume_1d, IntegralOperator::dirichlet>;
template class Westergaard<model_type::volume_2d, IntegralOperator::dirichlet>;
/* -------------------------------------------------------------------------- */
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
diff --git a/src/model/westergaard.hh b/src/model/westergaard.hh
index acd2f85c..2272861b 100644
--- a/src/model/westergaard.hh
+++ b/src/model/westergaard.hh
@@ -1,87 +1,93 @@
/**
* @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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __WESTERGAARD_HH__
#define __WESTERGAARD_HH__
/* -------------------------------------------------------------------------- */
#include "be_engine.hh"
#include "fft_plan_manager.hh"
#include "grid_hermitian.hh"
#include "integral_operator.hh"
#include "model_type.hh"
#include "tamaas.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
/**
* @brief Operator based on Westergaard solution and the Dicrete Fourier
* Transform.
* This class is templated with model type to allow efficient storage of the
* influence coefficients.
* The integral operator is only applied to surface pressure/displacements,
* even for volume models.
*/
template <model_type mtype, IntegralOperator::kind otype>
class Westergaard : public IntegralOperator {
using trait = model_type_traits<mtype>;
static constexpr UInt dim = trait::dimension;
static constexpr UInt bdim = trait::boundary_dimension;
static constexpr UInt comp = trait::components;
public:
/// Constuctor: initalizes influence coefficients and allocates buffer
Westergaard(Model* model);
/// Get influence coefficients
const GridHermitian<Real, bdim>& getInfluence() const { return influence; }
/// Apply influence coefficients to input
void apply(GridBase<Real>& input, GridBase<Real>& ouput) const override;
/// Update the influence coefficients
void updateFromModel() override { initInfluence(); }
+ /// Kind
+ IntegralOperator::kind getKind() const override { return otype; }
+
+ /// Type
+ model_type getType() const override { return mtype; }
+
private:
/// Initialize influence coefficients
void initInfluence();
template <typename Functor>
void initFromFunctor(Functor func);
/// Apply a functor in Fourier space
template <typename Functor>
void fourierApply(Functor func, GridBase<Real>& in,
GridBase<Real>& out) const;
public:
GridHermitian<Real, bdim> influence;
mutable GridHermitian<Real, bdim> buffer;
};
__END_TAMAAS__
#endif // __WESTERGAARD_HH__
diff --git a/src/solvers/contact_solver.cpp b/src/solvers/contact_solver.cpp
index d2d313c0..40c77c5e 100644
--- a/src/solvers/contact_solver.cpp
+++ b/src/solvers/contact_solver.cpp
@@ -1,67 +1,67 @@
/**
* @file
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016-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 "contact_solver.hh"
#include "fft_plan_manager.hh"
#include <iomanip>
#include <iostream>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
ContactSolver::ContactSolver(Model& model, const GridBase<Real>& surface,
Real tolerance)
- : model(model), surface(), functional(model.getBEEngine()),
+ : model(model), surface(), functional(),
tolerance(tolerance), surface_stddev(std::sqrt(surface.var())) {
this->surface.wrap(surface);
- _gap = allocateGrid<true, Real>(model.getType(), model.getDiscretization(),
- model.getTraction().getNbComponents());
+ _gap = allocateGrid<true, Real>(model.getType(),
+ model.getBoundaryDiscretization());
}
/* -------------------------------------------------------------------------- */
ContactSolver::~ContactSolver() {
FFTPlanManager::get().clean(); /// TODO find a smarter way to clean up
}
/* -------------------------------------------------------------------------- */
void ContactSolver::printState(UInt iter, Real cost_f, Real error) {
// UInt precision = std::cout.precision();
if (iter % dump_frequency)
return;
std::cout << std::setw(5) << iter << " ";
std::cout << std::setw(15) << std::scientific << cost_f << " ";
std::cout << std::setw(15) << error << std::endl;
std::cout << std::fixed;
}
/* -------------------------------------------------------------------------- */
void ContactSolver::addFunctionalTerm(functional::Functional* func) {
functional.addFunctionalTerm(func, false);
}
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
diff --git a/src/solvers/polonsky_keer_rey.cpp b/src/solvers/polonsky_keer_rey.cpp
index 77a9216d..9c773293 100644
--- a/src/solvers/polonsky_keer_rey.cpp
+++ b/src/solvers/polonsky_keer_rey.cpp
@@ -1,295 +1,289 @@
/**
* @file
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016-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 "polonsky_keer_rey.hh"
#include "elastic_functional.hh"
#include "fft_plan_manager.hh"
#include "loop.hh"
#include "model_type.hh"
+#include <boost/preprocessor/seq.hpp>
#include <iomanip>
#include <iterator>
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
PolonskyKeerRey::PolonskyKeerRey(Model& model, const GridBase<Real>& surface,
Real tolerance, type variable_type,
type constraint_type)
: ContactSolver(model, surface, tolerance), variable_type(variable_type),
constraint_type(constraint_type) {
- if (model.getType() != model_type::basic_1d &&
- model.getType() != model_type::basic_2d) {
- TAMAAS_EXCEPTION(
- "Model type is not compatible with PolonskyKeerRey solver");
+
+#define SET_VIEWS(_, __, type) \
+ case type: { \
+ setViews<type>(); \
+ break; \
+ }
+
+ switch (model.getType()) {
+ BOOST_PP_SEQ_FOR_EACH(SET_VIEWS, ~, TAMAAS_MODEL_TYPES);
}
- auto pressure_ = &model.getTraction();
+#undef SET_VIEWS
+
+ search_direction = allocateGrid<true, Real>(
+ operation_type, model.getBoundaryDiscretization());
+
+ projected_search_direction = allocateGrid<true, Real>(
+ operation_type, model.getBoundaryDiscretization());
switch (variable_type) {
case gap: {
model.getBEEngine().registerDirichlet();
- primal = this->_gap.get();
- dual = pressure_;
- this->functional.addFunctionalTerm(new functional::ElasticFunctionalGap(
- model.getBEEngine(), this->surface),
- true);
+ primal = gap_view.get();
+ dual = pressure_view.get();
+ this->functional.addFunctionalTerm(
+ new functional::ElasticFunctionalGap(*integral_op, this->surface),
+ true);
break;
}
case pressure: {
model.getBEEngine().registerNeumann();
- primal = pressure_;
- dual = this->_gap.get();
+ primal = pressure_view.get();
+ dual = gap_view.get();
this->functional.addFunctionalTerm(
- new functional::ElasticFunctionalPressure(model.getBEEngine(),
- this->surface),
+ new functional::ElasticFunctionalPressure(*integral_op, this->surface),
true);
break;
}
}
-
- search_direction =
- allocateGrid<true, Real>(model.getType(), model.getDiscretization(),
- model.getTraction().getNbComponents());
-
- projected_search_direction =
- allocateGrid<true, Real>(model.getType(), model.getDiscretization(),
- model.getTraction().getNbComponents());
}
/* -------------------------------------------------------------------------- */
Real PolonskyKeerRey::solve(Real target) {
Real G = 0, Gold = 1, error = 0, error_norm = 1;
UInt n = 0;
bool delta = false;
*search_direction = 0;
//*dual = 0;
// Printing column headers
std::cout << std::setw(5) << "Iter"
<< " " << std::setw(15) << "Cost_f"
<< " " << std::setw(15) << "Error" << '\n'
<< std::fixed;
// Setting uniform value if constraint
if (constraint_type == variable_type && std::abs(primal->sum()) <= 0)
*primal = target;
else if (constraint_type == variable_type)
*primal *= target / primal->mean();
else if (constraint_type != variable_type)
*primal = this->surface_stddev;
do {
// Computing functional gradient
functional.computeGradF(*primal, *dual);
Real dbar = meanOnUnsaturated(*dual);
// Enforcing dual constraint via gradient
if (constraint_type != variable_type) {
*dual += 2 * target + dbar;
} else {
// Centering dual on primal > 0
*dual -= dbar;
}
// Computing gradient norm
G = computeSquaredNorm(*dual);
// Updating search direction (conjugate gradient)
updateSearchDirection(delta * G / Gold);
Gold = G;
// Computing critical step
Real tau = computeCriticalStep(target);
// Update primal variable
delta = updatePrimal(tau);
// We should scale to constraint
if (constraint_type == variable_type)
*primal *= target / primal->mean();
error = computeError() / error_norm;
Real cost_f = functional.computeF(*primal, *dual);
printState(n, cost_f, error);
} while (error > this->tolerance && n++ < this->max_iterations);
// Final update of dual variable
functional.computeGradF(*primal, *dual);
Real dual_min = dual->min();
*dual -= dual_min;
// Primal is pressure: need to make sure gap is positive
if (variable_type == pressure) {
- model.getDisplacement() = *dual;
- model.getDisplacement() += this->surface;
+ *displacement_view = *dual;
+ *displacement_view += this->surface;
}
// Primal is gap: need to make sure dual is positive (pressure + adhesion)
else {
- model.getDisplacement() = *primal;
- model.getDisplacement() += this->surface;
- this->model.solveDirichlet();
- model.getTraction() -= dual_min;
+ *displacement_view = *primal;
+ *displacement_view += this->surface;
+ integral_op->apply(*displacement_view, *pressure_view);
+ *pressure_view -= dual_min;
}
return error;
}
/* -------------------------------------------------------------------------- */
/**
* Computes \f$ \frac{1}{\mathrm{card}(\{p > 0\})} \sum_{\{p > 0\}}{f_i} \f$
*/
Real PolonskyKeerRey::meanOnUnsaturated(const GridBase<Real>& field) const {
// Sum on unsaturated
Real sum = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(Real & p, const Real& f) { return (p > 0) ? f : 0; },
*primal, field);
// Number of unsaturated points
UInt n_unsat = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(Real & p) -> UInt { return (p > 0); }, *primal);
return sum / n_unsat;
}
/* -------------------------------------------------------------------------- */
/**
* Computes \f$ \sum_{\{p > 0\}}{f_i^2} \f$
*/
Real PolonskyKeerRey::computeSquaredNorm(const GridBase<Real>& field) const {
return Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(Real & p, const Real& f) { return (p > 0) ? f * f : 0; },
*primal, field);
}
/* -------------------------------------------------------------------------- */
/**
* Computes \f$ \tau = \frac{ \sum_{\{p > 0\}}{q_i't_i} }{ \sum_{\{p > 0\}}{r_i'
* t_i} } \f$
*/
Real PolonskyKeerRey::computeCriticalStep(Real target) {
- switch (variable_type) {
- case gap:
- model.getBEEngine().solveDirichlet(*search_direction,
- *projected_search_direction);
- break;
- case pressure:
- model.getBEEngine().solveNeumann(*search_direction,
- *projected_search_direction);
- break;
- }
+ integral_op->apply(*search_direction, *projected_search_direction);
Real rbar = meanOnUnsaturated(*projected_search_direction);
if (constraint_type == variable_type)
*projected_search_direction -= rbar;
else {
*projected_search_direction += 2 * target + rbar;
}
Real num = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(Real & p, Real & q, Real & t) {
return (p > 0) ? q * t : 0;
},
*primal, *dual, *search_direction);
Real denum = Loop::reduce<operation::plus>(
[] CUDA_LAMBDA(Real & p, Real & r, Real & t) {
return (p > 0) ? r * t : 0;
},
*primal, *projected_search_direction, *search_direction);
return num / denum;
}
/* -------------------------------------------------------------------------- */
/**
* Update steps:
* 1. \f$\mathbf{p} = \mathbf{p} - \tau \mathbf{t} \f$
* 2. Truncate all \f$p\f$ negative
* 3. For all points in \f$I_\mathrm{na} = \{p = 0 \land q < 0 \}\f$ do \f$p_i =
* p_i - \tau q_i\f$
*/
bool PolonskyKeerRey::updatePrimal(Real step) {
UInt na_num = Loop::reduce<operation::plus>(
[step] CUDA_LAMBDA(Real & p, Real & q, Real & t) -> UInt {
p -= step * t; // Updating primal
if (p < 0)
p = 0; // Truncating neg values
if (p == 0 && q < 0) { // If non-admissible state
p -= step * q;
return 1;
} else
return 0;
},
*primal, *dual, *search_direction);
return na_num == 0;
}
/* -------------------------------------------------------------------------- */
/**
* Error is based on \f$ \sum{p_i q_i} \f$
*/
Real PolonskyKeerRey::computeError() {
/// Making sure dual respects constraint
*dual -= dual->min();
Real error = primal->dot(*dual);
Real norm = 1;
if (variable_type == pressure)
norm = std::abs(primal->sum() * this->surface_stddev);
else
norm = std::abs(dual->sum() * this->surface_stddev);
norm *= primal->getNbPoints();
return std::abs(error) / norm;
}
/* -------------------------------------------------------------------------- */
/**
* Do \f$\mathbf{t} = \mathbf{q}' + \delta \frac{R}{R_\mathrm{old}}\mathbf{t}
* \f$
*/
void PolonskyKeerRey::updateSearchDirection(Real factor) {
Loop::loop(
[factor] CUDA_LAMBDA(Real & p, Real & q, Real & t) {
t = (p > 0) ? q + factor * t : 0;
},
*primal, *dual, *search_direction);
}
/* -------------------------------------------------------------------------- */
void PolonskyKeerRey::addFunctionalTerm(functional::Functional* func) {
functional.addFunctionalTerm(func, false);
}
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
diff --git a/src/solvers/polonsky_keer_rey.hh b/src/solvers/polonsky_keer_rey.hh
index dd30bc65..759f4315 100644
--- a/src/solvers/polonsky_keer_rey.hh
+++ b/src/solvers/polonsky_keer_rey.hh
@@ -1,74 +1,130 @@
/**
* @file
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2016-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/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef __POLONSKY_KEER_REY_HH__
#define __POLONSKY_KEER_REY_HH__
/* -------------------------------------------------------------------------- */
#include "contact_solver.hh"
+#include "grid_view.hh"
#include "meta_functional.hh"
-#include "tamaas.hh"
+#include "westergaard.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_TAMAAS__
class PolonskyKeerRey : public ContactSolver {
public:
/// Types of algorithm (primal/dual) or constraint
enum type { gap, pressure };
public:
/// Constructor
PolonskyKeerRey(Model& model, const GridBase<Real>& surface, Real tolerance,
type variable_type, type constraint_type);
public:
/// Solve
Real solve(Real target);
/// Mean on unsaturated constraint zone
Real meanOnUnsaturated(const GridBase<Real>& field) const;
/// Compute squared norm
Real computeSquaredNorm(const GridBase<Real>& field) const;
/// Update search direction
void updateSearchDirection(Real factor);
/// Compute critical step
Real computeCriticalStep(Real target = 0);
/// Update primal and check non-admissible state
bool updatePrimal(Real step);
/// Compute error/stopping criterion
Real computeError();
/// Add term to functional
- void addFunctionalTerm(functional::Functional * func);
+ void addFunctionalTerm(functional::Functional* func);
+
+private:
+ /// Set correct views on normal traction and gap
+ template <model_type type>
+ void setViews();
protected:
type variable_type, constraint_type;
- GridBase<Real>* primal = nullptr;
- GridBase<Real>* dual = nullptr;
+ model_type operation_type;
+ GridBase<Real>* primal = nullptr; ///< non-owning pointer for primal varaible
+ GridBase<Real>* dual = nullptr; ///< non-owning pointer for dual variable
+
+ /// CG search direction
std::unique_ptr<GridBase<Real>> search_direction = nullptr;
+ /// Projected CG search direction
std::unique_ptr<GridBase<Real>> projected_search_direction = nullptr;
+ /// View on normal pressure, gap, displacement
+ std::unique_ptr<GridBase<Real>> pressure_view, gap_view,
+ displacement_view = nullptr;
+
+ /// Integral operator for gradient computation
+ std::unique_ptr<IntegralOperator> integral_op = nullptr;
};
+/* -------------------------------------------------------------------------- */
+template <model_type mtype>
+void PolonskyKeerRey::setViews() {
+ constexpr UInt dim = model_type_traits<mtype>::dimension;
+ constexpr UInt bdim = model_type_traits<mtype>::boundary_dimension;
+ constexpr UInt comp = model_type_traits<mtype>::components;
+
+ pressure_view = std::unique_ptr<GridBase<Real>>{
+ new GridView<Grid, Real, bdim, bdim>(model.getTraction(), {}, comp - 1)};
+ gap_view = std::unique_ptr<GridBase<Real>>{
+ new GridView<Grid, Real, bdim, bdim>(*this->_gap, {}, comp - 1)};
+ displacement_view =
+ std::unique_ptr<GridBase<Real>>{new GridView<Grid, Real, dim, bdim>(
+ model.getDisplacement(), model_type_traits<mtype>::indices,
+ comp - 1)};
+
+ // I noz is fugly - TODO factory?
+ if (bdim == 1) {
+ operation_type = model_type::basic_1d;
+ if (variable_type == gap)
+ integral_op = std::unique_ptr<IntegralOperator>{
+ new Westergaard<model_type::basic_1d, IntegralOperator::dirichlet>(
+ &this->model)};
+ else
+ integral_op = std::unique_ptr<IntegralOperator>{
+ new Westergaard<model_type::basic_1d, IntegralOperator::neumann>(
+ &this->model)};
+ } else if (bdim == 2) {
+ operation_type = model_type::basic_2d;
+ if (variable_type == gap)
+ integral_op = std::unique_ptr<IntegralOperator>{
+ new Westergaard<model_type::basic_2d, IntegralOperator::dirichlet>(
+ &this->model)};
+ else
+ integral_op = std::unique_ptr<IntegralOperator>{
+ new Westergaard<model_type::basic_2d, IntegralOperator::neumann>(
+ &this->model)};
+ }
+}
+
__END_TAMAAS__
#endif // __POLONSKY_KEER_REY_HH__
diff --git a/tests/test_westergaard.py b/tests/test_westergaard.py
index a9bc094a..383097be 100644
--- a/tests/test_westergaard.py
+++ b/tests/test_westergaard.py
@@ -1,62 +1,81 @@
#!/usr/bin/env python
# coding: utf-8
# -----------------------------------------------------------------------------
# @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/>.
# -----------------------------------------------------------------------------
from __future__ import division, print_function
import sys
import numpy as np
from numpy.linalg import norm
import tamaas as tm
import pytest
@pytest.fixture(scope="module",
params=[tm.PolonskyKeerRey.pressure,
tm.PolonskyKeerRey.gap])
def pkr(westergaard, request):
loads = {
tm.PolonskyKeerRey.pressure: westergaard.load,
tm.PolonskyKeerRey.gap: 0.06697415181446396
}
model = tm.ModelFactory.createModel(tm.model_type.basic_1d, [1.], [westergaard.n])
model.setElasticity(westergaard.e_star, 0)
solver = tm.PolonskyKeerRey(model, westergaard.surface, 1e-12,
request.param,
request.param)
solver.solve(loads[request.param])
return model, westergaard
def test_energy(pkr):
model, sol = pkr
traction, displacement = model.getTraction(), model.getDisplacement()
error = norm((traction - sol.pressure)*(displacement - sol.displacement))
error /= norm(sol.pressure * sol.displacement)
assert error < 4e-6
+def test_pkr_multi(westergaard):
+ model_basic = tm.ModelFactory.createModel(tm.model_type.basic_1d, [1.], [westergaard.n])
+ model_volume = tm.ModelFactory.createModel(tm.model_type.volume_1d, [1., 1.],
+ [20, westergaard.n])
+
+ pressures = {}
+ for model in [model_basic, model_volume]:
+ print(model)
+ solver = tm.PolonskyKeerRey(model, westergaard.surface, 1e-12,
+ tm.PolonskyKeerRey.pressure,
+ tm.PolonskyKeerRey.pressure)
+ solver.solve(westergaard.load)
+ pressures[model] = model.getTraction()
+
+ error = norm(pressures[model_basic] - pressures[model_volume][:, 1]) / westergaard.load
+ assert error < 1e-16
+
+
+
if __name__ == "__main__":
print("Test is meant to run with pytest")