diff --git a/CHANGELOG.md b/CHANGELOG.md
index 41c18e6..252d002 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,314 +1,317 @@
# Changelog
All notable changes to this project will be documented in this file.
The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html) for final versions and [PEP440](https://www.python.org/dev/peps/pep-0440/) in case intermediate versions need to be released (e.g. development version `2.2.3.dev1` or release candidates `2.2.3rc1`), or individual commits are packaged.
## Unrelease
### Added
- Added `tamaas.utils.publications` to print a list of publications relevant to
the parts of Tamaas used in a script. Call at the end of a script for an
exhaustive list.
- Added `tamaas.mpi.gather` and `tamaas.mpi.scatter` to help handling of 2D data
in MPI contexts.
- Added `getBoundaryFields()/boundary_fields` function/property to `Model`.
- Added constructor to `Model`, which is an alias to `ModelFactory.createModel`.
- Added convenience `tamaas.dumpers._helper.netCDFtoParaview` function to
convert model sequence to Paraview data file.
- Added dump of surface fields as `FieldData` in VTK files
+- Added `scipy.sparse.linalg.LinearOperator` interface to `Westergaard`. A
+ `LinearOperator` can be created by calling
+ `scipy.sparse.linalg.aslinearoperator` on the operator object.
## Changed
- Filter for boundary fields in dumper should be smarter.
- `read` functions in dumpers are now `@classmethod`
- Changed sign of plastic correction in KatoSaturated solver to match the
EPICSolver sign.
- Plastic correction of KatoSaturated is registered in the model as
"KatoSaturated::residual_displacement".
- Changed to modern Python packaging (PEP517) with `setup.cfg` and
`pyproject.toml`. A `setup.py` script is still provided for editable
installations, see [gh#7953](https://github.com/pypa/pip/issues/7953). Once
setuptools has stable support for
[PEP621](https://peps.python.org/pep-0621/), `setup.cfg` will be merged into
`pyproject.toml`. Note that source distributions that can be built with
[`build`](https://pypi.org/project/build/) are not true source distributions,
since they contain the compiled Python module of Tamaas.
## Fixed
- Fixed `tamaas.dumpers.NetCDFDumper` to dump in MPI context.
## v2.4.0 -- 2022-03-22
### Added
- Added a `tamaas.utils` module.
- Added `tamaas.utils.load_path` generator, which yields model objects for a
sequence of applied loads.
- Added `tamaas.utils.seeded_surfaces` generator, which yields surfaces for a
sequence of random seeds.
- Added `tamaas.utils.hertz_surface` which generates a parabolic (Hertzian)
surface.
- Added `tamaas.utils.log_context` context manager which locally sets the log
level for Tamaas' logger.
- Added `tamaas.compute.from_voigt` to convert Voigt/Mendel fields to dense
tensor fields.
- Added a deep-copy function to `ModelFactory` to copy `Model` objects. Use
`model_copy = copy.deepcopy(model)` in Python to make a copy. Currently only
copies registered fields, same as dumpers/readers.
- Added a `read_sequence` method for dumpers to read model frames.
- Automatic draft release on Zenodo.
### Changed
- `*_ROOT` variables for vendored dependencies GoogleTest and pybind11 now
default to empty strings, so that the vendored trees in `third-party/` are not
selected by default. This is so system packages are given priority. Vendored
dependency submodules will eventually be deprecated.
### Fixed
- Fixed an issue with `scons -c` when GTest was missing.
## v2.3.1 -- 2021-11-08
### Added
- Now using `clang-format`, `clang-tidy` and `flake8` for linting
- Pre-built Tamaas images can be pulled with `docker pull
registry.gitlab.com/tamaas/tamaas`
- Added a `--version` flag to the `tamaas` command
### Changed
- The root `Dockerfile` now compiles Tamaas, so using Tamaas in Docker is easier.
- The model attributes dumped to Numpy files are now written in a JSON-formatted
string to avoid unsafe loading/unpickling of objects.
- Removed the `build_doc` build option: now the doc targets are automatically
added if the dependencies are met, and built if `scons doc` is called.
- Removed the `use_googletest` build option: if tests are built and gtest is
present, the corresponding tests will be built
### Fixed
- The command `tamaas plot` gracefully fails if Matplotlib is missing
- Better heuristic to guess which fields are defined on boundary in dumpers
(still not perfect and may give false negatives)
## v2.3.0 -- 2021-06-15
### Added
- Added `read()` method to dumpers to create a model from a dump file
- `getClusters()` can be called in MPI contact with partial contact maps
- Added a JSON encoder class for models and a JSON dumper
- CUDA compatibility is re-established, but has not been tested
- Docstrings in the Python bindings for many classes/methods
### Changed
- Tamaas version numbers are now managed by
[versioneer](https://github.com/python-versioneer/python-versioneer). This
means that Git tags prefixed with `v` (e.g. `v2.2.3`) carry meaning and
determine the version. When no tag is set, versioneer uses the last tag,
specifies the commit short hash and the distance to the last tag (e.g.
`2.2.2+33.ge314b0e`). This version string is used in the compiled library, the
`setup.py` script and the `__version__` variable in the python module.
- Tamaas migrated to [GitLab](https://gitlab.com/tamaas/tamaas)
- Continuous delivery has been implemented:
- the `master` branch will now automatically build and publish Python wheels
to `https://gitlab.com/api/v4/projects/19913787/packages/pypi/simple`. These
"nightly" builds can be installed with:
pip install \
--extra-index-url https://gitlab.com/api/v4/projects/19913787/packages/pypi/simple \
tamaas
- version tags pushed to `master` will automatically publish the wheels to
[PyPI](https://pypi.org/project/tamaas/)
### Deprecated
- The `finalize()` function is now deprecated, since it is automatically called
when the process terminates
- Python versions 3.5 and below are not supported anymore
### Fixed
- Fixed a host of dump read/write issues when model type was not `volume_*d`.
Dumper tests are now streamlined and systematic.
- Fixed a bug where `Model::solveDirichlet` would not compute correctly
- Fixed a bug where `Statistics::contact` would not normalize by the global
number of surface points
## v2.2.2 -- 2021-04-02
### Added
- Entry-point `tamaas` defines a grouped CLI for `examples/pipe_tools`. Try
executing `tamaas surface -h` from the command-line!
### Changed
- `CXXFLAGS` are now passed to the linker
- Added this changelog
- Using absolute paths for environmental variables when running `scons test`
- Reorganized documentation layout
- Gave the build system a facelift (docs are now generated directly with SCons
instead of a Makefile)
### Deprecated
- Python 2 support is discontinued. Version `v2.2.1` is the last PyPi build with
a Python 2 wheel.
- The scripts in `examples/pipe_tools` have been replaced by the `tamaas` command
### Fixed
- `UVWDumper` no longer imports `mpi4py` in sequential
- Compiling with different Thrust/FFTW backends
## v2.2.1 -- 2021-03-02
### Added
- Output registered fields and dumpers in `print(model)`
- Added `operator[]` to the C++ model class (for fields)
- Added `traction` and `displacement` properties to Python model bindings
- Added `operators` property to Python model bindings, which provides a
dict-like access to registered operators
- Added `shape` and `spectrum` to properties to Python surface generator
bindings
- Surface generator constructor accepts surface global shape as argument
- Choice of FFTW thread model
### Changed
- Tests use `/tmp` for temporary files
- Updated dependency versions (Thrust, Pybind11)
### Deprecated
- Most `get___()` and `set___()` in Python bindings have been deprecated. They
will generate a `DeprecationWarning`.
### Removed
- All legacy code
## v2.2.0 -- 2020-12-31
### Added
- More accurate function for computation of contact area
- Function to compute deviatoric of tensor fields
- MPI implementation
- Convenience `hdf5toVTK` function
- Readonly properties `shape`, `global_shape`, `boundary_shape` on model to give
shape information
### Changed
- Preprocessor defined macros are prefixed with `TAMAAS_`
- Moved `tamaas.to_voigt` to `tamaas.compute.to_voigt`
### Fixed
- Warning about deprecated constructors with recent GCC versions
- Wrong computation of grid strides
- Wrong computation of grid sizes in views
## v2.1.4 -- 2020-08-07
### Added
- Possibility to generate a static `libTamaas`
- C++ implementation of DFSANE solver
- Allowing compilation without OpenMP
### Changed
- NetCDF dumper writes frames to a single file
### Fixed
- Compatibility with SCons+Python 3
## v2.1.3 -- 2020-07-27
### Added
- Version number to `TamaasInfo`
### Changed
- Prepending root directory when generating archive
## v2.1.2 -- 2020-07-24
This release changes some core internals related to discrete Fourier transforms
for future MPI support.
### Added
- Caching `CXXFLAGS` in SCons build
- SCons shortcut to create code archive
- Test of the elastic-plastic contact solver
- Paraview data dumper (`.pvd` files)
- Compression for UVW dumper
- `__contains__` and `__iter__` Python bindings of model
- Warning message of possible overflow in Kelvin
### Changed
- Simplified `tamaas_info.cpp`, particularly the diff part
- Using a new class `FFTEngine` to manage discrete Fourier transforms. Plans are
re-used as much as possible with different data with the same shape. This is
in view of future MPI developments
- Redirecting I/O streams in solve functions so they can be used from Python
(e.g. in Jupyter notebooks)
- Calling `initialize()` and `finalize()` is no longer necessary
### Fixed
- Convergence issue with non-linear solvers
- Memory error in volume potentials
## v2.1.1 -- 2020-04-22
### Added
- SCons shortcut to run tests
### Fixed
- Correct `RPATH` for shared libraries
- Issues with SCons commands introduced in v2.1.0
- Tests with Python 2.7
## v2.1.0 -- 2020-04-17
### Added
- SCons shortcuts to build/install Tamaas and its components
- Selection of integration method for Kelvin operator
- Compilation option to remove the legacy part of Tamaas
- NetCDF dumper
### Fixed
- Link bug with clang
- NaNs in Kato saturated solver
## v2.0.0 -- 2019-11-11
First public release. Contains relatively mature elastic-plastic contact code.
diff --git a/python/cast.hh b/python/cast.hh
index 07b203a..e5cc85c 100644
--- a/python/cast.hh
+++ b/python/cast.hh
@@ -1,183 +1,196 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 CAST_HH
#define CAST_HH
/* -------------------------------------------------------------------------- */
-#include "grid_base.hh"
#include "grid.hh"
+#include "grid_base.hh"
#include "numpy.hh"
/* -------------------------------------------------------------------------- */
#include
#include
/* -------------------------------------------------------------------------- */
namespace pybind11 {
// Format descriptor necessary for correct wrap of tamaas complex type
template
struct format_descriptor<
tamaas::complex, detail::enable_if_t::value>> {
static constexpr const char c = format_descriptor::c;
static constexpr const char value[3] = {'Z', c, '\0'};
static std::string format() { return std::string(value); }
};
#ifndef PYBIND11_CPP17
template
constexpr const char format_descriptor<
tamaas::complex,
detail::enable_if_t::value>>::value[3];
#endif
namespace detail {
// declare tamaas complex as a complex type for pybind11
template
struct is_complex> : std::true_type {};
template
struct is_fmt_numeric,
detail::enable_if_t::value>>
: std::true_type {
static constexpr int index = is_fmt_numeric::index + 3;
};
static inline handle policy_switch(return_value_policy policy, handle parent) {
switch (policy) {
case return_value_policy::copy:
case return_value_policy::move:
return handle();
case return_value_policy::automatic_reference: // happens in python-derived
// classes
case return_value_policy::reference:
return none();
case return_value_policy::reference_internal:
return parent;
default:
TAMAAS_EXCEPTION("Policy is not handled");
}
}
template
handle grid_to_python(const tamaas::Grid& grid,
return_value_policy policy, handle parent) {
parent = policy_switch(policy, parent); // reusing variable
std::vector sizes(dim);
std::copy(grid.sizes().begin(), grid.sizes().end(), sizes.begin());
if (grid.getNbComponents() != 1)
sizes.push_back(grid.getNbComponents());
return array(sizes, grid.getInternalData(), parent).release();
}
+template
+handle grid_to_python(const tamaas::GridBase& grid,
+ return_value_policy policy, handle parent) {
+ parent = policy_switch(policy, parent); // reusing variable
+ std::vector sizes = {grid.dataSize()};
+ return array(sizes, grid.getInternalData(), parent).release();
+}
+
/**
* Type caster for grid classes
* inspired by https://tinyurl.com/y8m47qh3 from T. De Geus
* and pybind11/eigen.h
*/
template class G, typename T,
tamaas::UInt dim>
struct type_caster> {
using type = G;
using array_type =
array_t;
public:
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_TYPE_CASTER(type, _("GridWrap"));
/**
* 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 (!array_type::check_(src) or !convert)
return false;
auto buf = array_type::ensure(src);
- if (!buf) return false;
+ if (!buf)
+ return false;
value.move(tamaas::wrap::GridNumpy>(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.
*/
static handle cast(const type& src, return_value_policy policy,
handle parent) {
return grid_to_python(src, policy, parent);
}
};
/**
* Type caster for GridBase classes
*/
template
struct type_caster> {
using type = tamaas::GridBase;
using array_type =
array_t;
public:
// NOLINTNEXTLINE(readability-else-after-return)
PYBIND11_TYPE_CASTER(type, _("GridBaseWrap"));
bool load(handle src, bool convert) {
if (!array_type::check_(src) or !convert)
return false;
auto buf = array_type::ensure(src);
- if (!buf) return false;
+ if (!buf)
+ return false;
value.move(tamaas::wrap::GridBaseNumpy(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 auto& conv = dynamic_cast&>(src); \
- return grid_to_python(conv, policy, parent); \
+ const auto* conv = dynamic_cast*>(&src); \
+ if (conv) \
+ return grid_to_python(*conv, policy, parent); \
+ else \
+ return grid_to_python(src, policy, parent); \
}
switch (src.getDimension()) {
BOOST_PP_SEQ_FOR_EACH(GRID_BASE_CASE, ~, (1)(2)(3));
default:
- return nullptr;
+ return grid_to_python(src, policy, parent);
}
#undef GRID_BASE_CASE
}
};
} // namespace detail
} // namespace pybind11
#endif // CAST_HH
diff --git a/python/wrap/model.cpp b/python/wrap/model.cpp
index b883fcb..3f4a636 100644
--- a/python/wrap/model.cpp
+++ b/python/wrap/model.cpp
@@ -1,499 +1,508 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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;
struct model_operator_accessor {
Model& m;
decltype(auto) get(const std::string& name) {
return m.getIntegralOperator(name);
}
};
/// Wrap functional classes
void wrapFunctionals(py::module& mod) {
py::class_,
functional::wrap::PyFunctional>
func(mod, "Functional");
func.def(py::init<>())
.def("computeF", &functional::Functional::computeF,
"Compute functional value")
.def("computeGradF", &functional::Functional::computeGradF,
"Compute functional gradient");
py::class_ adh(mod, "AdhesionFunctional",
func);
adh.def_property("parameters", &functional::AdhesionFunctional::getParameters,
&functional::AdhesionFunctional::setParameters,
"Parameters dictionary")
.def("setParameters", [](functional::AdhesionFunctional& f,
const std::map& m) {
TAMAAS_DEPRECATE("setParameters()", "the parameters property");
f.setParameters(m);
});
py::class_(
mod, "ExponentialAdhesionFunctional", adh,
"Potential of the form F = -γ·exp(-g/ρ)")
.def(py::init&>(), "surface"_a);
py::class_(
mod, "MaugisAdhesionFunctional", adh,
"Cohesive zone potential F = H(g - ρ)·γ/ρ")
.def(py::init&>(), "surface"_a);
py::class_(
mod, "SquaredExponentialAdhesionFunctional", adh,
"Potential of the form F = -γ·exp(-0.5·(g/ρ)²)")
.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) {
TAMAAS_DEPRECATE("apply()", "the () operator");
auto in = instanciateFromNumpy(input);
auto out = instanciateFromNumpy(output);
op.apply(*in, *out);
})
.def(TAMAAS_DEPRECATE_ACCESSOR(getModel, IntegralOperator, "model"),
py::return_value_policy::reference)
.def(TAMAAS_DEPRECATE_ACCESSOR(getKind, IntegralOperator, "kind"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getType, IntegralOperator, "type"))
.def(
"__call__",
[](IntegralOperator& op, numpy input, numpy output) {
auto in = instanciateFromNumpy(input);
auto out = instanciateFromNumpy(output);
op.apply(*in, *out);
},
"Apply the integral operator")
.def("updateFromModel", &IntegralOperator::updateFromModel,
"Resets internal persistent variables from the model")
.def_property_readonly("kind", &IntegralOperator::getKind)
.def_property_readonly("model", &IntegralOperator::getModel)
- .def_property_readonly("type", &IntegralOperator::getType);
+ .def_property_readonly("type", &IntegralOperator::getType)
+ .def_property_readonly("shape", &IntegralOperator::matvecShape)
+ .def(
+ "matvec",
+ [](const IntegralOperator& op, numpy X) -> GridBase {
+ GridBaseNumpy x(X);
+ auto y = op.matvec(x);
+ return y;
+ },
+ py::return_value_policy::move);
py::enum_(mod, "integration_method",
"Integration method used for the computation "
"of volumetric Fourier operators")
.value("linear", integration_method::linear,
"No approximation error, O(N₁·N₂·N₃) time complexity, may cause "
"float overflow/underflow")
.value("cutoff", integration_method::cutoff,
"Approximation, O(sqrt(N₁²+N₂²)·N₃²) time complexity, no "
"overflow/underflow risk");
}
/// Wrap BEEngine classes
void wrapBEEngine(py::module& mod) {
py::class_(mod, "BEEngine")
.def("solveNeumann", &BEEngine::solveNeumann)
.def("solveDirichlet", &BEEngine::solveDirichlet)
.def("registerNeumann", &BEEngine::registerNeumann)
.def("registerDirichlet", &BEEngine::registerDirichlet)
.def(TAMAAS_DEPRECATE_ACCESSOR(getModel, BEEngine, "model"),
py::return_value_policy::reference)
.def_property_readonly("model", &BEEngine::getModel);
}
template
void wrapModelTypeTrait(py::module& mod) {
using trait = model_type_traits;
py::class_(mod, trait::repr)
.def_property_readonly_static(
"dimension", [](py::object) { return trait::dimension; },
"Dimension of computational domain")
.def_property_readonly_static(
"components", [](py::object) { return trait::components; },
"Number of components of vector fields")
.def_property_readonly_static(
"boundary_dimension",
[](py::object) { return trait::boundary_dimension; },
"Dimension of boundary of computational domain")
.def_property_readonly_static(
"voigt", [](py::object) { return trait::voigt; },
"Number of components of symmetrical tensor fields")
.def_property_readonly_static("indices",
[](py::object) { return trait::indices; });
}
/// Wrap Models
void wrapModelClass(py::module& mod) {
py::enum_(mod, "model_type")
.value("basic_1d", model_type::basic_1d,
"Normal contact with 1D interface")
.value("basic_2d", model_type::basic_2d,
"Normal contact with 2D interface")
.value("surface_1d", model_type::surface_1d,
"Normal & tangential contact with 1D interface")
.value("surface_2d", model_type::surface_2d,
"Normal & tangential contact with 2D interface")
.value("volume_1d", model_type::volume_1d,
"Contact with volumetric representation and 1D interface")
.value("volume_2d", model_type::volume_2d,
"Contact with volumetric representation and 2D interface");
auto trait_mod = mod.def_submodule("_type_traits");
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
wrapModelTypeTrait(trait_mod);
py::class_(mod, "_model_operator_acessor")
.def(py::init())
.def(
"__getitem__",
[](model_operator_accessor& acc, std::string name) {
try {
return acc.get(name);
} catch (std::out_of_range&) {
throw py::key_error(name);
}
},
py::return_value_policy::reference_internal)
.def("__contains__",
[](model_operator_accessor& acc, std::string key) {
const auto ops = acc.m.getIntegralOperators();
return std::find(ops.begin(), ops.end(), key) != ops.end();
})
.def(
"__iter__",
[](const model_operator_accessor& acc) {
const auto& ops = acc.m.getIntegralOperatorsMap();
return py::make_key_iterator(ops.cbegin(), ops.cend());
},
py::keep_alive<0, 1>());
py::class_(mod, "Model")
.def(py::init(py::overload_cast&,
const std::vector&>(
&ModelFactory::createModel)))
.def_property_readonly("type", &Model::getType)
.def_property("E", &Model::getYoungModulus, &Model::setYoungModulus,
"Young's modulus")
.def_property("nu", &Model::getPoissonRatio, &Model::setPoissonRatio,
"Poisson's ratio")
.def_property_readonly("mu", &Model::getShearModulus, "Shear modulus")
.def_property_readonly("E_star", &Model::getHertzModulus,
"Contact (Hertz) modulus")
.def_property_readonly("be_engine", &Model::getBEEngine,
"Boundary element engine")
.def(
"setElasticity",
[](Model& m, Real E, Real nu) {
TAMAAS_DEPRECATE("setElasticity()", "the E and nu properties");
m.setElasticity(E, nu);
},
"E"_a, "nu"_a)
.def(TAMAAS_DEPRECATE_ACCESSOR(getHertzModulus, Model, "E_star"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getYoungModulus, Model, "E"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getShearModulus, Model, "mu"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getPoissonRatio, Model, "nu"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getTraction, Model, "traction"),
py::return_value_policy::reference_internal)
.def(TAMAAS_DEPRECATE_ACCESSOR(getDisplacement, Model, "displacement"),
py::return_value_policy::reference_internal)
.def(TAMAAS_DEPRECATE_ACCESSOR(getSystemSize, Model, "system_size"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getDiscretization, Model, "shape"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getBoundarySystemSize, Model,
"boundary_system_size"))
.def(TAMAAS_DEPRECATE_ACCESSOR(getBoundaryDiscretization, Model,
"boundary_shape"))
.def("solveNeumann", &Model::solveNeumann,
"Solve surface tractions -> displacements")
.def("solveDirichlet", &Model::solveDirichlet,
"Solve surface displacemnts -> tractions")
.def("dump", &Model::dump, "Write model data to registered dumpers")
.def("addDumper", &Model::addDumper, "dumper"_a, py::keep_alive<1, 2>(),
"Register a dumper")
.def(
"getBEEngine",
[](Model& m) -> decltype(m.getBEEngine()) {
TAMAAS_DEPRECATE("getBEEngine()", "the be_engine property");
return m.getBEEngine();
},
py::return_value_policy::reference_internal)
.def(
"getIntegralOperator",
[](const Model& m, std::string name) {
TAMAAS_DEPRECATE("getIntegralOperator()", "the operators property");
return m.getIntegralOperator(name);
},
"operator_name"_a, py::return_value_policy::reference_internal)
.def(
"registerField",
[](Model& m, std::string name, numpy field) {
TAMAAS_DEPRECATE("registerField()", "the [] operator");
auto f = instanciateFromNumpy(field);
m.registerField(name, std::move(f));
},
"field_name"_a, "field"_a, py::keep_alive<1, 3>())
.def(
"getField",
[](const Model& m, std::string name) -> decltype(m.getField(name)) {
TAMAAS_DEPRECATE("getField()", "the [] operator");
return m.getField(name);
},
"field_name"_a, py::return_value_policy::reference_internal)
.def(
"getFields",
[](const Model& m) {
TAMAAS_DEPRECATE("getFields()", "list(model)");
return m.getFields();
},
"Return fields list")
.def(
"applyElasticity",
[](Model& model, numpy stress, numpy strain) {
auto out = instanciateFromNumpy(stress);
auto in = instanciateFromNumpy(strain);
model.applyElasticity(*out, *in);
},
"Apply Hooke's law")
// Python magic functions
.def("__repr__",
[](const Model& m) {
std::stringstream ss;
ss << m;
return ss.str();
})
.def(
"__getitem__",
[](const Model& m, std::string key) -> decltype(m[key]) {
try {
return m[key];
} catch (std::out_of_range&) {
throw py::key_error(key);
}
},
py::return_value_policy::reference_internal, "Get field")
.def(
"__setitem__",
[](Model& m, std::string name, numpy field) {
auto f = instanciateFromNumpy(field);
m.registerField(name, std::move(f));
},
py::keep_alive<1, 3>(), "Register new field")
.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>(), "Test field existence")
.def(
"__iter__",
[](const Model& m) {
const auto& fields = m.getFieldsMap();
return py::make_key_iterator(fields.cbegin(), fields.cend());
},
py::keep_alive<0, 1>(), "Iterator on fields")
.def(
"__copy__",
[](const Model&) {
throw std::runtime_error("__copy__ not implemented");
},
"Shallow copy of model. Not implemented.")
.def(
"__deepcopy__",
[](const Model& m, py::dict) { return ModelFactory::createModel(m); },
"Deep copy of model.")
.def_property_readonly("boundary_fields", &Model::getBoundaryFields)
.def_property_readonly(
"operators", [](Model& m) { return model_operator_accessor{m}; },
"Returns a dict-like object allowing access to the model's "
"integral "
"operators")
// More python-like access to model properties
.def_property_readonly("shape", &Model::getDiscretization,
"Discretization (local in MPI environment)")
.def_property_readonly("global_shape", &Model::getGlobalDiscretization,
"Global discretization (in MPI environement)")
.def_property_readonly("boundary_shape",
&Model::getBoundaryDiscretization,
"Number of points on boundary")
.def_property_readonly("system_size", &Model::getSystemSize,
"Size of physical domain")
.def_property_readonly("boundary_system_size",
&Model::getBoundarySystemSize,
"Physical size of surface")
.def_property_readonly("traction",
(const GridBase& (Model::*)() const) &
Model::getTraction,
"Surface traction field")
.def_property_readonly("displacement",
(const GridBase& (Model::*)() const) &
Model::getDisplacement,
"Displacement field");
py::class_>(
mod, "ModelDumper")
.def(py::init<>())
.def("dump", &ModelDumper::dump, "model"_a, "Dump model")
.def(
"__lshift__",
[](ModelDumper& dumper, Model& model) { dumper << model; },
"Dump model");
}
/// Wrap factory for models
void wrapModelFactory(py::module& mod) {
py::class_(mod, "ModelFactory")
.def_static(
"createModel",
py::overload_cast&,
const std::vector&>(
&ModelFactory::createModel),
"model_type"_a, "system_size"_a, "global_discretization"_a,
R"-(Create a new model of a given type, physical size and *global*
discretization.
:param model_type: the type of desired model
:param system_size: the physical size of the domain in each direction
:param global_discretization: number of points in each direction)-")
.def_static("createModel",
py::overload_cast(&ModelFactory::createModel),
"model"_a, "Create a deep copy of a model.")
.def_static("createResidual", &ModelFactory::createResidual, "model"_a,
"sigma_y"_a, "hardening"_a = 0.,
R"-(Create an isotropic linear hardening residual.
:param model: the model on which to define the residual
:param sigma_y: the (von Mises) yield stress
:param hardening: the hardening modulus)-")
.def_static("registerVolumeOperators",
&ModelFactory::registerVolumeOperators, "model"_a,
"Register Boussinesq and Mindlin operators to model.");
}
/// 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)
.def_property("yield_stress", &Residual::getYieldStress,
&Residual::setYieldStress)
.def_property("hardening_modulus", &Residual::getHardeningModulus,
&Residual::setHardeningModulus)
.def_property_readonly("model", &Residual::getModel);
}
void wrapModel(py::module& mod) {
wrapBEEngine(mod);
wrapModelClass(mod);
wrapModelFactory(mod);
wrapFunctionals(mod);
wrapResidual(mod);
wrapIntegralOperator(mod);
}
} // namespace wrap
} // namespace tamaas
diff --git a/src/model/integral_operator.hh b/src/model/integral_operator.hh
index 043ad8d..597472a 100644
--- a/src/model/integral_operator.hh
+++ b/src/model/integral_operator.hh
@@ -1,83 +1,95 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 INTEGRAL_OPERATOR_HH
#define INTEGRAL_OPERATOR_HH
/* -------------------------------------------------------------------------- */
#include "grid_base.hh"
#include "model_type.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
class Model;
/**
* @brief Generic class for integral operators
*/
class IntegralOperator {
public:
/// Kind of operator
enum kind { neumann, dirichlet, dirac };
public:
/// Constructor
IntegralOperator(Model* model) : model(model) {}
/// Virtual destructor for subclasses
virtual ~IntegralOperator() = default;
/// Apply operator on input
virtual void apply(GridBase& input, GridBase& output) const = 0;
/// Apply operator on filtered input
virtual void applyIf(GridBase& input, GridBase& output,
const std::function& /*unused*/) const {
apply(input, output);
}
/// 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;
/// Norm
virtual Real getOperatorNorm() {
TAMAAS_EXCEPTION("operator does not implement norm");
}
+ /// Dense shape (for Scipy integration)
+ virtual std::pair matvecShape() const {
+ TAMAAS_EXCEPTION("operator does not implement "
+ "scipy.sparse.linalg.LinearOperator interface");
+ }
+
+ /// matvec definition
+ virtual GridBase matvec(GridBase& /*x*/) const {
+ TAMAAS_EXCEPTION("operator does not implement "
+ "scipy.sparse.linalg.LinearOperator interface");
+ }
+
protected:
Model* model = nullptr;
};
/* -------------------------------------------------------------------------- */
/* Printing IntegralOperator::kind to string */
/* -------------------------------------------------------------------------- */
std::ostream& operator<<(std::ostream& o, const IntegralOperator::kind& val);
} // namespace tamaas
#endif // INTEGRAL_OPERATOR_HH
diff --git a/src/model/westergaard.cpp b/src/model/westergaard.cpp
index 00463db..d2ffc7c 100644
--- a/src/model/westergaard.cpp
+++ b/src/model/westergaard.cpp
@@ -1,277 +1,305 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 "grid_view.hh"
#include "influence.hh"
#include "loop.hh"
#include "model.hh"
#include "static_types.hh"
#include "westergaard.hh"
+
+#include
+#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
template
Westergaard::Westergaard(Model* model)
: IntegralOperator(model), influence(), buffer(),
engine(FFTEngine::makeEngine()) {
// Copy sizes
auto bdisc = model->getBoundaryDiscretization();
auto boundary_hermitian_sizes =
GridHermitian::hermitianDimensions(
bdisc);
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();
}
/* -------------------------------------------------------------------------- */
namespace detail {
template
struct boundary_fft_helper {
template
static void backwardTransform(FFTEngine& e, Buffer&& buffer, Out&& out) {
auto view = make_view(out, 0);
e.backward(view, buffer);
}
};
template
struct boundary_fft_helper {
template
static void backwardTransform(FFTEngine& e, Buffer&& buffer, Out&& out) {
e.backward(out, buffer);
}
};
} // namespace detail
template
template
void Westergaard::fourierApply(Functor func, GridBase& in,
GridBase& out) const try {
auto& i = dynamic_cast&>(in);
auto& full_out = dynamic_cast&>(out);
engine->forward(i, buffer);
/// Applying influence
func(buffer, influence);
/// Backward fourier transform on boundary only
detail::boundary_fft_helper::backwardTransform(*engine, buffer,
full_out);
} catch (const std::bad_cast& e) {
TAMAAS_EXCEPTION("Neumann and dirichlet types do not match model type");
}
/* -------------------------------------------------------------------------- */
template
template
void Westergaard::initFromFunctor(Functor func) {
// Compute wavevectors for influence
auto wavevectors = FFTEngine::template computeFrequencies(
influence.sizes());
// Get boundary physical size
auto system_size = this->model->getBoundarySystemSize();
VectorProxy domain(system_size[0]);
// Normalize wavevectors
wavevectors *= 2 * M_PI;
wavevectors /= domain;
// Apply functor
Loop::loop(func, range>(wavevectors),
range>(influence));
if (mpi::rank() == 0) {
// Set fundamental mode to zero
MatrixProxy mat(influence(0));
mat = 0;
}
}
/* -------------------------------------------------------------------------- */
/// \cond DO_NOT_DOCUMENT
#define NEUMANN_BASIC(type) \
template <> \
void Westergaard::initInfluence() { \
auto E_star = model->getHertzModulus(); \
auto basic = [E_star] CUDA_LAMBDA(VectorProxy q, \
MatrixProxy k) { \
k(0, 0) = 2. / (E_star * q.l2norm()); \
}; \
initFromFunctor(basic); \
}
#define DIRICHLET_BASIC(type) \
template <> \
void Westergaard::initInfluence() { \
auto E_star = model->getHertzModulus(); \
auto basic = [E_star] CUDA_LAMBDA(VectorProxy q, \
MatrixProxy k) { \
k(0, 0) = 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::initInfluence() {
auto E = model->getYoungModulus();
auto nu = model->getPoissonRatio();
const Complex I(0, 1);
auto surface = [E, nu, I] CUDA_LAMBDA(VectorProxy q_,
MatrixProxy F) {
Real q = q_(0);
q *= 1. / q_.l2norm();
F(0, 0) = 2 * (1 + nu) * (1 - nu * q * q);
F(1, 1) = 2 * (1 - nu * nu);
F(0, 1) = I * q * (1 + nu) * (1 - 2 * nu);
F(1, 0) = -F(0, 1);
F *= 1. / (E * q_.l2norm());
};
initFromFunctor(surface);
}
/* -------------------------------------------------------------------------- */
template <>
void Westergaard::initInfluence() {
auto E = model->getYoungModulus();
auto nu = model->getPoissonRatio();
const Complex I(0, 1);
auto surface = [E, nu, I] CUDA_LAMBDA(VectorProxy q_,
MatrixProxy F) {
Vector q(q_);
q *= 1 / q_.l2norm();
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_.l2norm());
};
initFromFunctor(surface);
}
/* -------------------------------------------------------------------------- */
template
void Westergaard::initInfluence() {}
/// \endcond
/* ------------------------------------------------------------------------ */
template
void Westergaard::apply(GridBase& input,
GridBase& output) const {
auto apply = [](decltype(buffer)& buffer,
const decltype(influence)& influence) {
Loop::loop(
[] CUDA_LAMBDA(VectorProxy i,
MatrixProxy F) { i = F * i; },
range>(buffer),
range>(influence));
};
fourierApply(apply, input, output);
}
/* -------------------------------------------------------------------------- */
template
Real Westergaard::getOperatorNorm() {
constexpr UInt comp = model_type_traits::components;
Real _norm = 0.0;
_norm = Loop::reduce(
[] CUDA_LAMBDA(MatrixProxy m) {
Real n = thrust::norm(m.l2squared());
return std::isnan(n) ? 0 : n;
},
range>(influence));
auto size = model->getSystemSize();
auto disc = model->getDiscretization();
auto dim = model_type_traits::dimension;
/// TODO: why?
switch (dim) {
case 1:
_norm /= (size[0] / disc[0]) * (size[0] / disc[0]);
break;
default:
for (UInt i = 0; i < dim; i++) {
_norm /= size[i] / disc[i];
}
}
return std::sqrt(_norm);
}
+/* -------------------------------------------------------------------------- */
+template
+std::pair Westergaard::matvecShape() const {
+ auto n = model->getBoundaryDiscretization();
+ auto N = std::accumulate(n.begin(), n.end(), comp, std::multiplies<>());
+ return std::make_pair(N, N);
+}
+
+/* -------------------------------------------------------------------------- */
+template
+GridBase Westergaard::matvec(GridBase& X) const {
+ auto n = model->getBoundaryDiscretization();
+ // Creating correct input
+ Grid x(n, comp);
+ thrust::copy(X.begin(), X.end(), x.begin());
+
+ // Creating correct output
+ if (bdim != dim)
+ n.insert(n.begin(), 1);
+
+ Grid y(n, comp);
+ apply(x, y);
+ return std::move(y);
+}
+
/* -------------------------------------------------------------------------- */
/* Template instanciation */
/* -------------------------------------------------------------------------- */
template class Westergaard;
template class Westergaard;
template class Westergaard;
template class Westergaard;
template class Westergaard;
template class Westergaard;
template class Westergaard;
template class Westergaard;
template class Westergaard;
template class Westergaard;
template class Westergaard;
template class Westergaard;
/* -------------------------------------------------------------------------- */
} // namespace tamaas
/* -------------------------------------------------------------------------- */
diff --git a/src/model/westergaard.hh b/src/model/westergaard.hh
index 880c621..ebe9168 100644
--- a/src/model/westergaard.hh
+++ b/src/model/westergaard.hh
@@ -1,90 +1,96 @@
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* 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 WESTERGAARD_HH
#define WESTERGAARD_HH
/* -------------------------------------------------------------------------- */
+#include "fft_engine.hh"
#include "grid_hermitian.hh"
#include "integral_operator.hh"
#include "model_type.hh"
-#include "fft_engine.hh"
#include "tamaas.hh"
/* -------------------------------------------------------------------------- */
namespace 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
class Westergaard : public IntegralOperator {
using trait = model_type_traits;
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& getInfluence() const { return influence; }
/// Apply influence coefficients to input
void apply(GridBase& input, GridBase& output) 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; }
/// Initialize influence coefficients
void initInfluence();
template
void initFromFunctor(Functor func);
/// Apply a functor in Fourier space
template
void fourierApply(Functor func, GridBase& in,
GridBase& out) const;
/// Compute L_2 norm of influence functions
Real getOperatorNorm() override;
+ /// Dense shape
+ std::pair matvecShape() const override;
+
+ /// Dense matvec
+ GridBase matvec(GridBase& x) const override;
+
public:
GridHermitian influence;
mutable GridHermitian buffer;
mutable std::unique_ptr engine;
};
} // namespace tamaas
#endif // WESTERGAARD_HH
diff --git a/tests/SConscript b/tests/SConscript
index 7b78016..ad75d07 100644
--- a/tests/SConscript
+++ b/tests/SConscript
@@ -1,223 +1,224 @@
# -*- mode:python; coding: utf-8 -*-
# vim: set ft=python:
#
# Copyright (©) 2016-2022 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 .
from __future__ import print_function
from SCons.Script import Split, Copy, Dir, Import
from detect import FindGTest, FindPybind11
# ------------------------------------------------------------------------------
def copyComStr(env, main):
if 'SHCXXCOMSTR' in main:
env['CXXCOMSTR'] = main['SHCXXCOMSTR']
if 'SHLINKCOMSTR' in main:
env['LINKCOMSTR'] = main['SHLINKCOMSTR']
# ------------------------------------------------------------------------------
def make_python_tests(env):
"""Copy python tests to build directory"""
test_env = env.Clone()
test_files = Split("""
test_hertz.py
+ test_matvec.py
test_westergaard.py
test_patch_westergaard.py
test_patch_plasticity.py
test_surface.py
test_hertz_disp.py
test_hertz_kato.py
test_saturated_pressure.py
test_flood_fill.py
test_integral_operators.py
test_dumper.py
test_tangential.py
test_boussinesq_surface.py
test_voigt.py
test_memory.py
test_copy.py
test_epic.py
test_publications.py
fftfreq.py
conftest.py
pytest.ini
""")
if env['use_mpi']:
test_files += ['test_mpi_routines.py', 'mpi_routines.py']
src_dir = "#/tests"
targets = [
test_env.Command(file, test_env.File(file, src_dir),
Copy("$TARGET", "$SOURCE"))
for file in test_files
]
test_env = env.Clone(tools=[pybind11])
# Helper module for integral operators
test_env['SHLIBPREFIX'] = ''
test_env.PrependUnique(LIBS=['Tamaas'])
register = test_env.Pybind11Module(
target="register_integral_operators",
source=["register_integral_operators.cpp"])
Import('libTamaas')
test_env.Depends(register, libTamaas)
targets.append(register)
return targets
# ------------------------------------------------------------------------------
def compile_google_test(env, gtest_path):
gtest_obj = env.Object('gtest.o',
[env.File("src/gtest-all.cc", gtest_path)])
return env.StaticLibrary('gtest', gtest_obj)
# ------------------------------------------------------------------------------
def make_google_tests(env):
gtest_dir = Dir(env['GTEST_ROOT'])
gtest_env = env.Clone(CPPPATH=[gtest_dir],
CXXFLAGS=['-pthread', '-isystem',
env.Dir('include', gtest_dir).path])
FindGTest(gtest_env)
if not gtest_env.get('has_gtest'):
return []
libgtest = None
# Hugly hack to detect if we need to compile gtest submodule
if env['GTEST_ROOT'] == '#third-party/googletest/googletest':
gtest_path = str(gtest_dir)
libgtest = compile_google_test(gtest_env, gtest_path)
env.AppendUnique(CXXFLAGS=gtest_env['CXXFLAGS'])
env.PrependUnique(LIBS=['Tamaas', env.subst('python${py_version}')])
google_test_files = Split("""
test_grid.cpp
test_loop.cpp
test_model.cpp
test_static_types.cpp
test_integration.cpp
""")
if env['use_fftw']:
google_test_files += ['test_fftw.cpp']
if env['use_cuda']:
google_test_files += ['test_cufft.cpp']
# Necessary for the tests that use pybind11 calls to python
uses = []
if env['build_python']:
google_test_files.append('test_fftfreq.cpp')
uses = ['TAMAAS_USE_PYTHON']
if env['use_mpi']:
google_test_files.append('test_mpi.cpp')
defines = env['CPPDEFINES']
if type(defines) is not list:
defines = [defines]
gtest_main = env.Object("tamaas_gtest_main.o", 'tamaas_gtest_main.cc',
CPPDEFINES=defines + uses)
gtest_all = env.Program('test_gtest_all', google_test_files + [gtest_main],
LIBS=(env['LIBS'] + ['gtest']))
Import('libTamaas')
env.Depends(gtest_all, libTamaas)
env.Depends(gtest_all, libgtest)
copy_gtest = gtest_env.Command('test_gtest.py', '#tests/test_gtest.py',
Copy('$TARGET', '$SOURCE'))
return [gtest_all, copy_gtest]
# ------------------------------------------------------------------------------
def make_bare_tests(env):
rough = env.Program("test_rough_surface.cpp")
Import('libTamaas')
env.Depends(rough, libTamaas)
return [rough]
# ------------------------------------------------------------------------------
Import('main_env')
# Setup of test environment
test_env = main_env.Clone()
test_env.AppendUnique(
LIBPATH=['.', '../src'],
RPATH=["'$$$$ORIGIN/../src'"],
LINKFLAGS=['-pthread'],
)
test_env.PrependUnique(LIBS=['Tamaas'])
# Building tests that do not require any third party
targets = make_bare_tests(test_env)
# Build tests that required python bindings
if test_env['build_python']:
FindPybind11(test_env)
test_env.Tool(pybind11)
test_env.ParseConfig("${py_exec}-config --ldflags")
test_env['CCFLAGS'] = []
targets += make_python_tests(test_env)
# Building google tests
targets += make_google_tests(test_env)
# Target alias to build tests
main_env.Alias('build-tests', targets)
# Check if pytest is installed
conf = Configure(test_env,
custom_tests={'CheckPythonModule': CheckPythonModule})
has_pytest = conf.CheckPythonModule('pytest')
conf.Finish()
# Define a command to execute tests
if has_pytest:
pytest_env = test_env.Clone()
pytest_env['pythonpath'] = pytest_env.subst('${build_dir}/python')
pytest_env['ld_library_path'] = pytest_env.subst('${build_dir}/src')
pytest_env.PrependENVPath('PYTHONPATH', pytest_env['pythonpath'])
pytest_env.PrependENVPath('LD_LIBRARY_PATH', pytest_env['ld_library_path'])
# Setting a moderate thread number
pytest_env['ENV']['OMP_NUM_THREADS'] = "1"
pytest_env['PYTESTOPTS'] = ['${"-v" if verbose else "-q"}']
test_target = pytest_env.Command(
'.phony_test', targets,
'echo $$PYTHONPATH; echo $$LD_LIBRARY_PATH; '
'${py_exec} -m pytest $PYTESTOPTS ${TARGET.dir}')
main_env.Alias('test', test_target)
else:
# We still define a target here so that `scons test` still works
dummy_command(main_env, 'test',
'Cannot run tests: pytest is not installed')
diff --git a/tests/test_matvec.py b/tests/test_matvec.py
new file mode 100644
index 0000000..fc79569
--- /dev/null
+++ b/tests/test_matvec.py
@@ -0,0 +1,63 @@
+# -*- coding: utf-8 -*-
+#
+# Copyright (©) 2016-2022 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 .
+
+import pytest
+import numpy as np
+import tamaas as tm
+
+from numpy.testing import assert_allclose
+from test_dumper import model_fixture # noqa
+
+
+linalg = pytest.importorskip('scipy.sparse.linalg')
+
+
+def test_matvec_westergaard(model_fixture):
+ m = model_fixture
+ m.be_engine.registerNeumann()
+ op = m.operators['Westergaard::neumann']
+ linop = linalg.aslinearoperator(op)
+
+ x = [np.linspace(0, L, N, endpoint=False)
+ for (L, N) in zip(m.boundary_system_size, m.boundary_shape)]
+
+ coords = np.meshgrid(*x, indexing='ij')
+
+ t = m.traction
+ res = np.zeros_like(m.displacement)
+ volume = len(m.shape) != len(m.boundary_shape)
+ basic = m.type in (tm.model_type.basic_1d, tm.model_type.basic_2d)
+
+ # Adapting to special cases
+ if basic:
+ t = t[..., np.newaxis]
+ res = res[..., np.newaxis]
+
+ t[:] = 0
+ t[..., -1] = np.sin(2 * np.pi * coords[0])
+
+ # Computing using native API
+ op(t, res)
+
+ # Only look at surface displacement
+ if volume:
+ res = res[0]
+
+ # Computing using LinearOperator API
+ y = np.reshape(linop * np.ravel(t), res.shape)
+ assert_allclose(res, y)