diff --git a/python/wrap/model.cpp b/python/wrap/model.cpp
index 2afa2f6..e8bd58e 100644
--- a/python/wrap/model.cpp
+++ b/python/wrap/model.cpp
@@ -1,320 +1,360 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#include "model.hh"
#include "adhesion_functional.hh"
#include "functional.hh"
#include "integral_operator.hh"
#include "model_dumper.hh"
#include "model_extensions.hh"
#include "model_factory.hh"
#include "numpy.hh"
#include "residual.hh"
#include "wrap.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
namespace wrap {
using namespace py::literals;
+struct model_operator_accessor {
+ Model& m;
+
+ decltype(auto) get(std::string name) {
+ return m.getIntegralOperator(std::move(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)
.def("computeGradF", &functional::Functional::computeGradF);
py::class_ adh(mod, "AdhesionFunctional",
func);
adh.def_property("parameters", &functional::AdhesionFunctional::getParameters,
&functional::AdhesionFunctional::setParameters)
// legacy wrapper code
.def("setParameters", &functional::AdhesionFunctional::setParameters);
py::class_(
mod, "ExponentialAdhesionFunctional", adh)
.def(py::init&>(), "surface"_a);
py::class_(
mod, "MaugisAdhesionFunctional", adh)
.def(py::init&>(), "surface"_a);
py::class_(
mod, "SquaredExponentialAdhesionFunctional", adh)
.def(py::init&>(), "surface"_a);
}
template
std::unique_ptr> instanciateFromNumpy(numpy& num) {
std::unique_ptr> result = nullptr;
switch (num.ndim()) {
case 2:
result = std::make_unique>>(num);
return result;
case 3:
result = std::make_unique>>(num);
return result;
case 4:
result = std::make_unique>>(num);
return result;
default:
TAMAAS_EXCEPTION("instanciateFromNumpy expects the last dimension of numpy "
"array to be the number of components");
}
}
/// Wrap IntegralOperator
void wrapIntegralOperator(py::module& mod) {
py::class_(mod, "IntegralOperator")
.def("apply",
[](IntegralOperator& op, numpy input, numpy output) {
auto in = instanciateFromNumpy(input);
auto out = instanciateFromNumpy(output);
op.apply(*in, *out);
})
.def("updateFromModel", &IntegralOperator::updateFromModel)
.def("getModel", &IntegralOperator::getModel,
py::return_value_policy::reference)
.def("getKind", &IntegralOperator::getKind)
.def("getType", &IntegralOperator::getType);
py::enum_(mod, "integration_method")
.value("linear", integration_method::linear)
.value("cutoff", integration_method::cutoff);
}
/// Wrap BEEngine classes
void wrapBEEngine(py::module& mod) {
py::class_(mod, "BEEngine")
.def("solveNeumann", &BEEngine::solveNeumann)
.def("solveDirichlet", &BEEngine::solveDirichlet)
.def("getModel", &BEEngine::getModel, py::return_value_policy::reference)
.def("registerNeumann", &BEEngine::registerNeumann)
.def("registerDirichlet", &BEEngine::registerDirichlet);
}
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; })
.def_property_readonly_static(
"components", [](py::object) { return trait::components; })
.def_property_readonly_static(
"boundary_dimension",
[](py::object) { return trait::boundary_dimension; })
.def_property_readonly_static("voigt",
[](py::object) { return trait::voigt; })
.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)
.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);
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_property_readonly("type", &Model::getType)
.def("setElasticity", &Model::setElasticity, "E"_a, "nu"_a)
.def_property("E", &Model::getYoungModulus, &Model::setYoungModulus)
.def_property("nu", &Model::getPoissonRatio, &Model::setPoissonRatio)
.def("getHertzModulus", &Model::getHertzModulus)
.def("getYoungModulus", &Model::getYoungModulus)
.def("getShearModulus", &Model::getShearModulus)
.def("getPoissonRatio", &Model::getPoissonRatio)
.def("getTraction", (GridBase & (Model::*)()) & Model::getTraction,
py::return_value_policy::reference_internal)
.def("getDisplacement",
(GridBase & (Model::*)()) & Model::getDisplacement,
py::return_value_policy::reference_internal)
.def("getSystemSize", &Model::getSystemSize)
.def("getDiscretization", &Model::getDiscretization)
.def("getBoundarySystemSize", &Model::getBoundarySystemSize)
.def("getBoundaryDiscretization", &Model::getBoundaryDiscretization)
.def("solveNeumann", &Model::solveNeumann)
.def("solveDirichlet", &Model::solveDirichlet)
.def("dump", &Model::dump)
.def("addDumper", &Model::addDumper, "dumper"_a, py::keep_alive<1, 2>())
.def("setDumper",
[](Model&, const ModelDumper&) {
TAMAAS_EXCEPTION("setDumper() is not a member of Model; use "
"addDumper() instead");
})
.def("getBEEngine", &Model::getBEEngine,
py::return_value_policy::reference_internal)
.def("getIntegralOperator", &Model::getIntegralOperator,
"operator_name"_a, py::return_value_policy::reference_internal)
.def("registerField",
[](Model& m, std::string name, numpy field) {
auto f = instanciateFromNumpy(field);
m.registerField(name, std::move(f));
},
"field_name"_a, "field"_a, py::keep_alive<1, 3>())
.def("getField",
(GridBase & (Model::*)(const std::string&)) & Model::getField,
"field_name"_a, py::return_value_policy::reference_internal)
.def("getFields", &Model::getFields)
.def("applyElasticity",
[](Model& model, numpy stress, numpy strain) {
auto out = instanciateFromNumpy(stress);
auto in = instanciateFromNumpy(strain);
model.applyElasticity(*out, *in);
})
- // Python magic functions
+ // Python magic functions
.def("__repr__",
[](const Model& m) {
std::stringstream ss;
ss << m;
return ss.str();
})
.def("__getitem__",
[](const Model& m, std::string key) -> const GridBase& {
try {
return m.getField(key);
} catch (std::out_of_range&) {
- throw py::key_error();
+ throw py::key_error(key);
}
},
py::return_value_policy::reference_internal)
+ .def("__setitem__",
+ [](Model& m, std::string name, numpy field) {
+ auto f = instanciateFromNumpy(field);
+ m.registerField(name, std::move(f));
+ },
+ py::keep_alive<1, 3>())
.def("__contains__",
[](const Model& m, std::string key) {
const auto fields = m.getFields();
return std::find(fields.begin(), fields.end(), key) !=
fields.end();
},
py::keep_alive<0, 1>())
.def("__iter__",
[](const Model& m) {
const auto& fields = m.getFieldsMap();
return py::make_key_iterator(fields.cbegin(), fields.cend());
},
py::keep_alive<0, 1>())
-
- // More python-like access to model properties
+ .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)
.def_property_readonly("global_shape", &Model::getGlobalDiscretization)
.def_property_readonly("boundary_shape",
&Model::getBoundaryDiscretization)
.def_property_readonly("system_size", &Model::getSystemSize);
py::class_>(
mod, "ModelDumper")
.def(py::init<>())
.def("dump", &ModelDumper::dump, "model"_a)
.def("__lshift__",
[](ModelDumper& dumper, Model& model) { dumper << model; });
}
/// Wrap factory for models
void wrapModelFactory(py::module& mod) {
py::class_(mod, "ModelFactory")
.def_static("createModel", &ModelFactory::createModel, "model_type"_a,
"system_size"_a, "global_discretization"_a)
- .def_static("createResidual", &ModelFactory::createResidual,
- "model"_a, "sigma_y"_a, "hardening"_a = 0.)
+ .def_static("createResidual", &ModelFactory::createResidual, "model"_a,
+ "sigma_y"_a, "hardening"_a = 0.)
.def_static("registerVolumeOperators",
&ModelFactory::registerVolumeOperators, "model"_a);
}
/// Wrap residual class
void wrapResidual(py::module& mod) {
// TODO adapt to n-dim
py::class_(mod, "Residual")
.def(py::init())
.def("computeResidual",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.computeResidual(*in);
})
.def("computeStress",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.computeStress(*in);
})
.def("updateState",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.updateState(*in);
})
.def("computeResidualDisplacement",
[](Residual& res, numpy& x) {
auto in = instanciateFromNumpy(x);
res.computeResidualDisplacement(*in);
})
.def("applyTangent",
[](Residual& res, numpy& output, numpy& input,
numpy& current_strain_inc) {
auto out = instanciateFromNumpy(output);
auto in = instanciateFromNumpy(input);
auto inc = instanciateFromNumpy(current_strain_inc);
res.applyTangent(*out, *in, *inc);
},
"output"_a, "input"_a, "current_strain_increment"_a)
.def("getVector", &Residual::getVector,
py::return_value_policy::reference_internal)
.def("getPlasticStrain", &Residual::getPlasticStrain,
py::return_value_policy::reference_internal)
.def("getStress", &Residual::getStress,
py::return_value_policy::reference_internal)
.def("setIntegrationMethod", &Residual::setIntegrationMethod, "method"_a,
"cutoff"_a = 1e-12)
.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/SConscript b/src/SConscript
index 4ec5151..853ca6c 100644
--- a/src/SConscript
+++ b/src/SConscript
@@ -1,188 +1,189 @@
# -*- mode:python; coding: utf-8 -*-
# vim: set ft=python:
# @file
# @section LICENSE
#
# Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see .
import os
Import('main_env')
def prepend(path, list):
return [os.path.join(path, x) for x in list]
env = main_env.Clone()
# Core
core_list = """
fft_engine.cpp
fftw_engine.cpp
grid.cpp
grid_hermitian.cpp
statistics.cpp
surface.cpp
tamaas.cpp
legacy_types.cpp
loop.cpp
computes.cpp
logger.cpp
""".split()
if env['legacy_bem']:
core_list.append('surface_statistics.cpp')
core_list = prepend('core', core_list)
info_file = main_env.Substfile('tamaas_info.cpp', 'tamaas_info.cpp.in')
core_list.append(info_file)
# Lib roughcontact
generator_list = """
surface_generator.cpp
surface_generator_filter.cpp
surface_generator_filter_fft.cpp
surface_generator_random_phase.cpp
isopowerlaw.cpp
regularized_powerlaw.cpp
""".split()
generator_list = prepend('surface', generator_list)
# Lib PERCOLATION
percolation_list = """
flood_fill.cpp
""".split()
percolation_list = prepend('percolation', percolation_list)
# BEM PERCOLATION
bem_list = """
bem_kato.cpp
bem_polonski.cpp
bem_gigi.cpp
bem_gigipol.cpp
bem_penalty.cpp
bem_uzawa.cpp
bem_fft_base.cpp
bem_functional.cpp
bem_meta_functional.cpp
elastic_energy_functional.cpp
exponential_adhesion_functional.cpp
squared_exponential_adhesion_functional.cpp
maugis_adhesion_functional.cpp
complimentary_term_functional.cpp
bem_grid.cpp
bem_grid_polonski.cpp
bem_grid_kato.cpp
bem_grid_teboulle.cpp
bem_grid_condat.cpp
""".split()
bem_list = prepend('bem', bem_list)
# Model
model_list = """
model.cpp
model_factory.cpp
model_type.cpp
model_template.cpp
integral_operator.cpp
be_engine.cpp
westergaard.cpp
elastic_functional.cpp
meta_functional.cpp
adhesion_functional.cpp
volume_potential.cpp
kelvin.cpp
mindlin.cpp
boussinesq.cpp
+hooke.cpp
elasto_plastic/isotropic_hardening.cpp
elasto_plastic/residual.cpp
integration/element.cpp
""".split()
model_list = prepend('model', model_list)
# Solvers
solvers_list = """
contact_solver.cpp
polonsky_keer_rey.cpp
kato_saturated.cpp
kato.cpp
beck_teboulle.cpp
condat.cpp
polonsky_keer_tan.cpp
ep_solver.cpp
dfsane_solver.cpp
epic.cpp
""".split()
solvers_list = prepend('solvers', solvers_list)
# GPU API
gpu_list = """
fftransform_cufft.cpp
""".split()
gpu_list = prepend('gpu', gpu_list)
# MPI API
mpi_list = """
fftw_mpi_engine.cpp
mpi_interface.cpp
""".split()
mpi_list = prepend('mpi', mpi_list)
# Assembling total list
rough_contact_list = \
core_list + generator_list + percolation_list + model_list + solvers_list
# Adding legacy code
if env['legacy_bem']:
rough_contact_list += bem_list
# Adding GPU if needed
if env['backend'] == 'cuda':
rough_contact_list += gpu_list
# Adding MPI if needed
if env['use_mpi']:
rough_contact_list += mpi_list
# Adding extra warnings for Tamaas base lib
env.AppendUnique(CXXFLAGS=['-Wextra'])
# Allowing libTamaas.so to find libs in the same directory
env.AppendUnique(RPATH=["'$$$$ORIGIN'"])
# Build static library for packaging
if env['build_static_lib']:
env.AppendUnique(CXXFLAGS='-fPIC')
libTamaas = env.StaticLibrary('Tamaas', rough_contact_list)
# Build shared library (default)
else:
libTamaas = env.SharedLibrary('Tamaas', rough_contact_list,
SHLIBVERSION=env['version'])
# Defining alias for cpp builds
main_env.Alias('build-cpp', libTamaas)
# Specify install target to install lib
lib_prefix = os.path.join(env['prefix'], 'lib')
lib_install = env.InstallVersionedLib(target=lib_prefix,
source=libTamaas)
main_env.Alias('install-lib', lib_install)
# Export target for use in python builds
Export('libTamaas')
diff --git a/src/model/hooke.cpp b/src/model/hooke.cpp
new file mode 100644
index 0000000..c48607d
--- /dev/null
+++ b/src/model/hooke.cpp
@@ -0,0 +1,70 @@
+/**
+ * @file
+ * @section LICENSE
+ *
+ * Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
+ * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Affero General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public License
+ * along with this program. If not, see .
+ *
+ */
+/* -------------------------------------------------------------------------- */
+#include "hooke.hh"
+#include "influence.hh"
+#include "loop.hh"
+#include "model.hh"
+/* -------------------------------------------------------------------------- */
+
+namespace tamaas {
+
+template
+void Hooke::apply(GridBase& strain, GridBase& stress) const {
+ constexpr auto dim = model_type_traits::dimension;
+ const influence::ElasticHelper elasticity(
+ this->model->getShearModulus(), this->model->getPoissonRatio());
+
+ // Checking incompressibility
+ if (elasticity.nu == 0.5)
+ TAMAAS_EXCEPTION("Incompressibility error");
+
+ // Checking number of components
+ if (strain.getNbComponents() == dim * dim) {
+
+ Loop::loop([&elasticity](auto sigma,
+ auto epsilon) { sigma = elasticity(epsilon); },
+ range>(stress),
+ range>(strain));
+ }
+
+ // in case we have voigt notation
+ else if (strain.getNbComponents() == voigt_size::value) {
+ Loop::loop([&elasticity](auto sigma,
+ auto epsilon) { sigma = elasticity(epsilon); },
+ range>(stress),
+ range>(strain));
+ }
+
+ else
+ TAMAAS_EXCEPTION("Strain components do not match dimension");
+}
+
+/* -------------------------------------------------------------------------- */
+template class Hooke;
+template class Hooke;
+template class Hooke;
+template class Hooke;
+template class Hooke;
+template class Hooke;
+
+} // namespace tamaas
diff --git a/src/model/model_template.hh b/src/model/hooke.hh
similarity index 52%
copy from src/model/model_template.hh
copy to src/model/hooke.hh
index 0d4092d..a10e0b2 100644
--- a/src/model/model_template.hh
+++ b/src/model/hooke.hh
@@ -1,68 +1,56 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
-#ifndef __MODEL_TEMPLATE_HH__
-#define __MODEL_TEMPLATE_HH__
+#ifndef HOOKE_HH
+#define HOOKE_HH
/* -------------------------------------------------------------------------- */
-#include "grid_view.hh"
-#include "model.hh"
+#include "integral_operator.hh"
#include "model_type.hh"
/* -------------------------------------------------------------------------- */
+
namespace tamaas {
/**
- * @brief Model class templated with model type
- * Specializations of this class should take care of dimension specific code
+ * @brief Applies Hooke's law of elasticity
*/
template
-class ModelTemplate : public Model {
- using trait = model_type_traits;
- static constexpr UInt dim = trait::dimension;
- using ViewType =
- GridView;
-
+class Hooke : public IntegralOperator {
public:
- /// Constructor
- ModelTemplate(std::vector system_size,
- std::vector discretization);
-
- // Get model type
+ Hooke(Model* m) : IntegralOperator(m) {}
+ /// Type of underlying model
model_type getType() const override { return type; }
- std::vector getGlobalDiscretization() const override;
- std::vector getBoundaryDiscretization() const override;
- std::vector getBoundarySystemSize() const override;
+ /// Hooke's law computes stresses ~ Neumann
+ IntegralOperator::kind getKind() const override {
+ return IntegralOperator::dirac;
+ }
- void applyElasticity(GridBase& stress,
- const GridBase& strain) const override;
+ /// Does not update
+ void updateFromModel() override {}
-protected:
- virtual void initializeBEEngine();
-
-protected:
- std::unique_ptr displacement_view = nullptr;
- std::unique_ptr traction_view = nullptr;
+ /// Apply Hooke's tensor
+ void apply(GridBase& input, GridBase& output) const override;
};
} // namespace tamaas
-/* -------------------------------------------------------------------------- */
-#endif // __MODEL_TEMPLATE_HH__
+
+#endif // HOOKE_HH
diff --git a/src/model/influence.hh b/src/model/influence.hh
index a7fd31b..39fb980 100644
--- a/src/model/influence.hh
+++ b/src/model/influence.hh
@@ -1,280 +1,327 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#ifndef INFLUENCE_HH
#define INFLUENCE_HH
/* -------------------------------------------------------------------------- */
#include "loop.hh"
#include "static_types.hh"
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
namespace influence {
namespace {
const Complex I(0, 1); /// < imaginary unit
constexpr inline Real sign(bool upper) { return (upper) ? 1. : -1.; }
} // namespace
template // dim + 1 to help dim template deduction
inline Vector
computeD(const VectorProxy& q) {
constexpr Real s = -sign(conjugate);
const auto q_norm = q.l2norm();
return {{{Complex{0, s * q(0) / q_norm}, Complex{0, s * q(1) / q_norm},
Complex{1, 0}}}};
}
template // dim + 1 to help dim template deduction
inline Vector
computeD2(const VectorProxy& q) {
const auto q_norm = (q.l2norm() < 1e-16) ? 1 : q.l2norm();
return {
{{Complex{0, q(0) / q_norm}, Complex{0, q(1) / q_norm}, Complex{0, 0}}}};
}
+/* -------------------------------------------------------------------------- */
+/// Functor to apply Hooke's tensor
+template
+struct ElasticHelper {
+ ElasticHelper(Real mu, Real nu)
+ : mu(mu), nu(nu), lambda(2 * mu * nu / (1 - 2 * nu)) {}
+
+ template
+ Matrix, dim, dim>
+ operator()(const StaticMatrix& gradu) const {
+ auto trace = gradu.trace();
+ Matrix, dim, dim> sigma;
+ for (UInt i = 0; i < dim; ++i)
+ for (UInt j = 0; j < dim; ++j)
+ sigma(i, j) =
+ (i == j) * lambda * trace + mu * (gradu(i, j) + gradu(j, i));
+ return sigma;
+ }
+
+ template
+ SymMatrix, dim>
+ operator()(const StaticSymMatrix& eps) const {
+ SymMatrix, dim> sigma;
+ auto trace = eps.trace();
+ for (UInt i = 0; i < dim; ++i)
+ sigma(i) = lambda * trace + 2 * mu * eps(i);
+ for (UInt i = dim; i < voigt_size::value; ++i)
+ sigma(i) = 2 * mu * eps(i);
+ return sigma;
+ }
+
+ template
+ Matrix, dim, dim>
+ inverse(const StaticMatrix& sigma) const {
+ Matrix, dim, dim> epsilon;
+ auto trace = sigma.trace();
+ auto c1 = 1. / (2 * mu);
+ auto c2 = -lambda / (3 * lambda + 2 * mu);
+
+ for (UInt i = 0; i < dim; ++i)
+ for (UInt j = 0; j < dim; ++j)
+ epsilon(i, j) = c1 * (sigma(i, j) + c2 * trace * (i == j));
+ }
+
+ const Real mu, nu, lambda;
+};
+
/* -------------------------------------------------------------------------- */
/// Class for the Kelvin tensor
template
class Kelvin;
/**
* @brief Kelvin base tensor
* See arXiv:1811.11558 for details.
*/
template <>
class Kelvin<3, 0> {
protected:
static constexpr UInt dim = 3;
static constexpr UInt order = 0;
public:
Kelvin(Real mu, Real nu) : mu(mu), b(4 * (1 - nu)) {}
template
inline Vector applyU0(const VectorProxy& q,
const StaticVector& f) const {
// upper is ignored here
auto d2 = computeD2(q);
auto q_power = (apply_q_power) ? 1. / q.l2norm() : 1.;
d2 *= d2.dot(f);
d2 += b * f;
d2(2) -= f(2);
d2 *= q_power / (2. * mu * b);
return d2;
}
template
inline Vector applyU1(const VectorProxy& q,
const StaticVector& f) const {
auto d = computeD(q);
auto q_power = (apply_q_power) ? 1. / q.l2norm() : 1.;
d *= d.dot(f);
d *= q_power * sign(upper) / (2. * mu * b);
return d;
}
protected:
const Real mu, b;
};
/* -------------------------------------------------------------------------- */
/// Kelvin first derivative
template <>
class Kelvin<3, 1> : protected Kelvin<3, 0> {
using parent = Kelvin<3, 0>;
protected:
static constexpr UInt dim = parent::dim;
static constexpr UInt order = parent::order + 1;
public:
using parent::parent;
template
inline Vector
applyU0(const VectorProxy& q,
const StaticMatrix& f) const {
// apply_q_power is ignored here
auto tmp = f * computeD(q);
auto res = Kelvin<3, 0>::applyU0(q, tmp);
res *= -sign(upper);
Vector e3{{{0, 0, 1}}};
tmp.template mul(f, e3);
res += Kelvin<3, 0>::applyU1(q, tmp);
return res;
}
template
inline Vector
applyU1(const VectorProxy& q,
const StaticMatrix& f) const {
// apply_q_power is ignored here
auto tmp = f * computeD(q);
auto res = Kelvin<3, 0>::applyU1(q, tmp);
res *= -sign(upper);
return res;
}
};
/* -------------------------------------------------------------------------- */
/// Kelvin second derivative
template <>
class Kelvin<3, 2> : protected Kelvin<3, 1> {
using parent = Kelvin<3, 1>;
static constexpr UInt dim = parent::dim;
static constexpr UInt order = parent::order + 1;
public:
using parent::parent;
template
inline Matrix
applyU0(const VectorProxy& q,
const StaticMatrix& f) const {
auto tmp = Kelvin<3, 1>::applyU0(q, f);
auto q_power = (apply_q_power) ? q.l2norm() : 1.;
tmp *= -sign(upper);
auto res = outer(computeD(q), tmp);
Vector e3{{{0, 0, 1}}};
res += outer(e3, Kelvin<3, 1>::applyU1(q, f));
res *= -1. * q_power; // << -1 for grad_x
return res;
}
template
inline Matrix
applyU1(const VectorProxy& q,
const StaticMatrix& f) const {
auto tmp = Kelvin<3, 1>::applyU1(q, f);
auto q_power = (apply_q_power) ? q.l2norm() : 1.;
auto res = outer(computeD(q), tmp);
res *= -1. * -sign(upper) * q_power; // << -1 for grad_x
return res;
}
template
inline Matrix
applyDiscontinuityTerm(const StaticMatrix& f) const {
Vector e3{{{0, 0, 1}}};
auto res = (f * e3);
res *= -this->b;
res(2) += 2. * f(2, 2);
e3(2) *= -1 / (this->mu * this->b); // << -1 for grad_x
return outer(e3, res);
}
};
/* -------------------------------------------------------------------------- */
/// Class for the Boussinesq tensor
template
class Boussinesq;
template <>
class Boussinesq<3, 0> {
protected:
static constexpr UInt dim = 3;
static constexpr UInt order = 0;
public:
/// Constructor
Boussinesq(Real mu, Real nu)
: mu(mu), nu(nu), lambda(2 * mu * nu / (1 - 2 * nu)) {}
template
inline Vector
applyU0(const StaticVector& t,
const VectorProxy& q) const {
// This guy has the property Bt(q) = B(-q). Since we know we recieve -q as
// an argument, we can just apply B without transpose
auto dplus = computeD(q), dminus = computeD(q),
d2 = computeD2(q);
auto q_power = (apply_q_power) ? 1. / q.l2norm() : 1.;
auto res = 2. * t;
res += ((1. - 2. * nu) * dminus.dot(t)) * dplus;
res += d2.dot(t) * d2;
res(2) -= t(2);
res *= q_power / (2. * mu);
return res;
}
template
inline Vector
applyU1(const StaticVector& t,
const VectorProxy& q) const {
// This guy does not have the property that Bt(q) = B(-q), so we have to
// apply Bt(q) for real. This makes us correct (in the line below) for the
// fact that we recieve -q as an argument.
auto dplus = computeD(q);
auto q_power = (apply_q_power) ? 1. / q.l2norm() : 1.;
auto res = q_power * (dplus.dot(t) / (2. * mu)) * dplus;
return res;
}
protected:
const Real mu, nu, lambda;
};
/* -------------------------------------------------------------------------- */
/// Boussinesq first gradient
template <>
class Boussinesq<3, 1> : protected Boussinesq<3, 0> {
protected:
using parent = Boussinesq<3, 0>;
static constexpr UInt dim = parent::dim;
static constexpr UInt order = parent::order + 1;
public:
using parent::parent;
template
inline Matrix
applyU0(const StaticVector& t,
const VectorProxy& q) const {
// taking into account K(-q) convolution
auto dbar = computeD(q);
Vector e3{{{0, 0, 1}}};
auto res = outer(dbar, Boussinesq<3, 0>::applyU0(t, q));
res *= -1.;
res += outer(e3, Boussinesq<3, 0>::applyU1(t, q));
return res;
}
template
inline Matrix
applyU1(const StaticVector& t,
const VectorProxy& q) const {
// taking into account K(-q) convolution
auto dbar = computeD(q);
auto res = outer(dbar, Boussinesq<3, 0>::applyU1(t, q));
res *= -1;
return res;
}
};
} // namespace influence
/* -------------------------------------------------------------------------- */
} // namespace tamaas
#endif // INFLUENCE_HH
diff --git a/src/model/integral_operator.hh b/src/model/integral_operator.hh
index ca1fbf1..f499167 100644
--- a/src/model/integral_operator.hh
+++ b/src/model/integral_operator.hh
@@ -1,84 +1,84 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#ifndef 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 };
+ 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,
std::function) 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");
}
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/model.cpp b/src/model/model.cpp
index 87eaeff..65f6ccc 100644
--- a/src/model/model.cpp
+++ b/src/model/model.cpp
@@ -1,163 +1,191 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#include "model.hh"
#include "be_engine.hh"
#include "logger.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
Model::Model(std::vector system_size, std::vector discretization)
: system_size(std::move(system_size)),
discretization(std::move(discretization)) {}
/* -------------------------------------------------------------------------- */
Model::~Model() = default;
/* -------------------------------------------------------------------------- */
void Model::setElasticity(Real E_, Real nu_) {
setYoungModulus(E_);
setPoissonRatio(nu_);
updateOperators();
}
+/* -------------------------------------------------------------------------- */
+void Model::applyElasticity(GridBase& stress,
+ const GridBase& strain) const {
+ operators.at("hooke")->apply(const_cast&>(strain), stress);
+}
+
/* -------------------------------------------------------------------------- */
Real Model::getHertzModulus() const { return E / (1 - nu * nu); }
/* -------------------------------------------------------------------------- */
GridBase& Model::getTraction() { return getField("traction"); }
const GridBase& Model::getTraction() const {
return getField("traction");
}
/* -------------------------------------------------------------------------- */
GridBase& Model::getDisplacement() { return getField("displacement"); }
const GridBase& Model::getDisplacement() const {
return getField("displacement");
}
/* -------------------------------------------------------------------------- */
const std::vector& Model::getSystemSize() const { return system_size; }
/* -------------------------------------------------------------------------- */
const std::vector& Model::getDiscretization() const {
return discretization;
}
/* -------------------------------------------------------------------------- */
void Model::solveNeumann() {
engine->registerNeumann();
engine->solveNeumann(getField("traction"), getField("displacement"));
}
void Model::solveDirichlet() {
engine->registerDirichlet();
engine->solveDirichlet(getField("displacement"), getField("traction"));
}
/* -------------------------------------------------------------------------- */
IntegralOperator* Model::getIntegralOperator(const std::string& name) {
- return operators[name].get();
+ return operators.at(name).get();
+}
+
+/* -------------------------------------------------------------------------- */
+std::vector Model::getIntegralOperators() const {
+ std::vector keys;
+ keys.reserve(operators.size());
+
+ std::transform(operators.begin(), operators.end(), std::back_inserter(keys),
+ [](auto&& pair) { return pair.first; });
+ return keys;
}
/* -------------------------------------------------------------------------- */
void Model::updateOperators() {
for (auto& op : operators)
op.second->updateFromModel();
}
/* -------------------------------------------------------------------------- */
void Model::addDumper(std::shared_ptr dumper) {
this->dumpers.push_back(std::move(dumper));
}
void Model::dump() const {
for (auto& dumper : dumpers)
if (dumper)
dumper->dump(*this);
}
/* -------------------------------------------------------------------------- */
void Model::registerField(const std::string& name,
std::shared_ptr> field) {
fields[name] = std::move(field);
}
const GridBase& Model::getField(const std::string& name) const try {
return *fields.at(name);
} catch (std::out_of_range& e) {
Logger().get(LogLevel::warning)
<< "Field " << name << " not registered in model\n";
throw e;
}
GridBase& Model::getField(const std::string& name) try {
return *fields.at(name);
} catch (std::out_of_range& e) {
Logger().get(LogLevel::warning)
<< "Field " << name << " not registered in model\n";
throw e;
}
std::vector Model::getFields() const {
std::vector keys;
keys.reserve(fields.size());
std::transform(fields.begin(), fields.end(), std::back_inserter(keys),
[](auto&& pair) { return pair.first; });
return keys;
}
/* -------------------------------------------------------------------------- */
std::ostream& operator<<(std::ostream& o, const Model& _this) {
o << "Model<" << _this.getType() << "> (E = " << _this.getYoungModulus()
<< ", nu = " << _this.getPoissonRatio() << ")\n";
+ auto out_collec = [&o](auto&& collec) {
+ std::for_each(collec.begin(), collec.end() - 1,
+ [&o](const auto& x) { o << x << ", "; });
+ o << collec.back();
+ };
+
// Printing domain size
o << " - domain = [";
- std::for_each(_this.getSystemSize().begin(), _this.getSystemSize().end() - 1,
- [&o](const Real& x) { o << x << ", "; });
- o << _this.getSystemSize().back() << "]\n";
+ out_collec(_this.getSystemSize());
+ o << "]\n";
// Printing discretization
o << " - discretization = [";
- std::for_each(_this.getDiscretization().begin(),
- _this.getDiscretization().end() - 1,
- [&o](const UInt& x) { o << x << ", "; });
- o << _this.getDiscretization().back() << "]";
+ out_collec(_this.getDiscretization());
+ o << "]\n";
+
+ // Print fields
+ o << " - registered fields = [";
+ out_collec(_this.getFields());
+ o << "]\n";
+
+ o << " - registered operators = [";
+ out_collec(_this.getIntegralOperators());
+ o << "]\n";
return o;
}
/* -------------------------------------------------------------------------- */
} // namespace tamaas
diff --git a/src/model/model.hh b/src/model/model.hh
index 3c6a9a8..36930ad 100644
--- a/src/model/model.hh
+++ b/src/model/model.hh
@@ -1,241 +1,195 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#ifndef MODEL_HH
#define MODEL_HH
/* -------------------------------------------------------------------------- */
#include "be_engine.hh"
#include "grid_base.hh"
#include "integral_operator.hh"
#include "model_dumper.hh"
#include "model_type.hh"
#include "tamaas.hh"
#include
#include
#include
#include
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
-/**
- * @brief Class for linear isotropic elasticity
- */
-namespace influence {
-template
-struct ElasticHelper {
- ElasticHelper(Real mu, Real nu)
- : mu(mu), nu(nu), lambda(2 * mu * nu / (1 - 2 * nu)) {}
-
- template
- Matrix, dim, dim>
- operator()(const StaticMatrix& gradu) const {
- auto trace = gradu.trace();
- Matrix, dim, dim> sigma;
- for (UInt i = 0; i < dim; ++i)
- for (UInt j = 0; j < dim; ++j)
- sigma(i, j) =
- (i == j) * lambda * trace + mu * (gradu(i, j) + gradu(j, i));
- return sigma;
- }
-
- template
- SymMatrix, dim>
- operator()(const StaticSymMatrix& eps) const {
- SymMatrix, dim> sigma;
- auto trace = eps.trace();
- for (UInt i = 0; i < dim; ++i)
- sigma(i) = lambda * trace + 2 * mu * eps(i);
- for (UInt i = dim; i < voigt_size::value; ++i)
- sigma(i) = 2 * mu * eps(i);
- return sigma;
- }
-
- template
- Matrix, dim, dim>
- inverse(const StaticMatrix& sigma) const {
- Matrix, dim, dim> epsilon;
- auto trace = sigma.trace();
- auto c1 = 1. / (2 * mu);
- auto c2 = -lambda / (3 * lambda + 2 * mu);
-
- for (UInt i = 0; i < dim; ++i)
- for (UInt j = 0; j < dim; ++j)
- epsilon(i, j) = c1 * (sigma(i, j) + c2 * trace * (i == j));
- }
-
- const Real mu, nu, lambda;
-};
-} // namespace influence
-
/* -------------------------------------------------------------------------- */
/**
* @brief Model containing pressure and displacement
* This class is a container for the model fields. It is supposed to be
* dimension agnostic, hence the GridBase members.
*/
class Model {
protected:
/// Constructor
Model(std::vector system_size, std::vector discretization);
public:
/// Destructor
virtual ~Model();
public:
/// Set elasticity parameters
void setElasticity(Real E, Real nu);
/// Get Hertz contact modulus
Real getHertzModulus() const;
/// Get Young's modulus
Real getYoungModulus() const { return E; }
/// Get Poisson's ratio
Real getPoissonRatio() const { return nu; }
/// Get shear modulus
Real getShearModulus() const { return E / (2 * (1 + nu)); }
/// Set Young's modulus
void setYoungModulus(Real E_) {
if (E_ < 0)
TAMAAS_EXCEPTION("Elastic modulus should be positive");
this->E = E_;
updateOperators();
}
/// Set Poisson's ratio
void setPoissonRatio(Real nu_) {
if (nu_ > 0.5 or nu_ <= -1)
TAMAAS_EXCEPTION("Poisson ratio should be in ]-1, 0.5]");
this->nu = nu_;
updateOperators();
}
public:
- virtual void applyElasticity(GridBase& stress,
- const GridBase& strain) const = 0;
+ void applyElasticity(GridBase& stress,
+ const GridBase& strain) const;
public:
/// Get pressure
GridBase& getTraction();
/// Get pressure
const GridBase& getTraction() const;
/// Get displacement
GridBase& getDisplacement();
/// Get displacement
const GridBase& getDisplacement() const;
/// Get model type
virtual model_type getType() const = 0;
/// Get system physical size
const std::vector& getSystemSize() const;
/// Get boundary system physical size
virtual std::vector getBoundarySystemSize() const = 0;
/// Get discretization
const std::vector& getDiscretization() const;
/// Get discretization of global MPI system
virtual std::vector getGlobalDiscretization() const = 0;
/// Get boundary discretization
virtual std::vector getBoundaryDiscretization() const = 0;
/// Get boundary element engine
BEEngine& getBEEngine() {
TAMAAS_ASSERT(engine, "BEEngine was not initialized");
return *engine;
}
/// Solve Neumann problem using default neumann operator
void solveNeumann();
/// Solve Dirichlet problem using default dirichlet operator
void solveDirichlet();
public:
/// Register a new integral operator
template
IntegralOperator* registerIntegralOperator(const std::string& name) {
Logger().get(LogLevel::debug)
<< TAMAAS_DEBUG_MSG("registering operator " + name);
operators[name] = std::make_unique(this);
return operators[name].get();
}
/// Get a registerd integral operator
IntegralOperator* getIntegralOperator(const std::string& name);
+ /// Get list of integral operators
+ std::vector getIntegralOperators() const;
+ /// Get operators mapcar
+ const auto& getIntegralOperatorsMap() const { return operators; }
/// Tell operators to update their cache
void updateOperators();
public:
/// Register a field
void registerField(const std::string& name,
std::shared_ptr> field);
/// Get a field
const GridBase& getField(const std::string& name) const;
/// Get a non-const field
GridBase& getField(const std::string& name);
/// Get fields
std::vector getFields() const;
/// Get fields map
const auto& getFieldsMap() const { return fields; }
public:
/// Set the dumper object
void addDumper(std::shared_ptr dumper);
/// Dump the model
void dump() const;
protected:
Real E = 1, nu = 0;
std::vector system_size;
std::vector discretization;
std::unique_ptr engine = nullptr;
std::unordered_map> operators;
std::unordered_map>> fields;
std::vector> dumpers;
};
/* -------------------------------------------------------------------------- */
/* Output model to stream */
/* -------------------------------------------------------------------------- */
std::ostream& operator<<(std::ostream& o, const Model& _this);
/* -------------------------------------------------------------------------- */
/* Simpler grid allocation */
/* -------------------------------------------------------------------------- */
template
std::unique_ptr> allocateGrid(Model& model) {
return allocateGrid(
model.getType(), (boundary) ? model.getBoundaryDiscretization()
: model.getDiscretization());
}
template
std::unique_ptr> allocateGrid(Model& model, UInt nc) {
return allocateGrid(model.getType(),
(boundary)
? model.getBoundaryDiscretization()
: model.getDiscretization(),
nc);
}
} // namespace tamaas
#endif // MODEL_HH
diff --git a/src/model/model_template.cpp b/src/model/model_template.cpp
index 0ddde97..1e264e5 100644
--- a/src/model/model_template.cpp
+++ b/src/model/model_template.cpp
@@ -1,160 +1,162 @@
/**
* @file
* @section LICENSE
*
* Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see .
*
*/
/* -------------------------------------------------------------------------- */
#include "model_template.hh"
+#include "computes.hh"
+#include "hooke.hh"
#include "influence.hh"
#include "partitioner.hh"
#include "westergaard.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
+namespace detail {
+template
+class ComputeOperator : public IntegralOperator {
+public:
+ /// Constructor
+ ComputeOperator(Model* model) : IntegralOperator(model) {}
+
+ IntegralOperator::kind getKind() const override {
+ return IntegralOperator::dirac;
+ }
+ model_type getType() const override { return model->getType(); }
+ void updateFromModel() override {}
+
+ /// Apply functor
+ void apply(GridBase& in, GridBase& out) const override {
+ applyCompute(model->getType(), out, in);
+ }
+};
+} // namespace detail
+
template
ModelTemplate