diff --git a/python/wrap/test_features.cpp b/python/wrap/test_features.cpp index 06a1b63..7e864e5 100644 --- a/python/wrap/test_features.cpp +++ b/python/wrap/test_features.cpp @@ -1,75 +1,62 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "kelvin.hh" #include "model.hh" #include "model_type.hh" -#include "ve_engine.hh" +#include "volume_potential.hh" #include "wrap.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ namespace wrap { -void wrapVEEngine(py::module& mod) { - py::class_(mod, "VEEngine") - .def("applyVolumeForcePotential", &VEEngine::applyVolumeForcePotential) - .def("applyVolumeStressPotential", &VEEngine::applyVolumeStressPotential) - .def("applyHypersingularStressPotential", - &VEEngine::applyHypersingularStressPotential); -} - -template -void wrapVEEngineTmpl(py::module& mod) { - py::class_, VEEngine>(mod, "VEEngineTmpl"); -} -template +template void wrapKelvin(py::module& mod) { constexpr UInt dim = model_type_traits::dimension; - py::class_, VEEngine>(mod, "Kelvin") + py::class_>(mod, "Kelvin") .def(py::init()) - .def("applyVolumeForcePotential", - [](const Kelvin& engine, Grid& in, + .def("apply", + [](const Kelvin& engine, Grid& in, Grid& out) { - engine.applyVolumeForcePotential(in, out); + engine.apply(in, out); }); } /// Wrap temporary features for testing void wrapTestFeatures(py::module& mod) { auto test_module = mod.def_submodule("_test_features"); test_module.doc() = "Module for testing new features.\n" "DISCLAIMER: this API is subject to frequent and unannounced changes " "and should **not** be relied upon!"; - wrapVEEngine(test_module); - wrapKelvin(test_module); + wrapKelvin(test_module); } } __END_TAMAAS__ diff --git a/src/SConscript b/src/SConscript index d8d085c..6fc1464 100644 --- a/src/SConscript +++ b/src/SConscript @@ -1,115 +1,114 @@ import os Import('main_env') def prepend(path, list): return [os.path.join(path, x) for x in list] env = main_env # Core core_list = """ fft_plan_manager.cpp fftransform.cpp fftransform_fftw.cpp grid.cpp grid_hermitian.cpp statistics.cpp surface.cpp tamaas.cpp legacy_types.cpp loop.cpp """.split() core_list = prepend('core', core_list) core_list.append('tamaas_info.cpp') # Lib roughcontact generator_list = """ surface_generator.cpp surface_generator_filter.cpp surface_generator_filter_fft.cpp isopowerlaw.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 be_engine.cpp westergaard.cpp elastic_functional.cpp meta_functional.cpp adhesion_functional.cpp elasto_plastic_model.cpp -ve_engine.cpp +volume_potential.cpp kelvin.cpp -mindlin.cpp """.split() model_list = prepend('model', model_list) # Solvers solvers_list = """ contact_solver.cpp polonsky_keer_rey.cpp """.split() solvers_list = prepend('solvers', solvers_list) # GPU API gpu_list = """ fftransform_cufft.cpp """.split() gpu_list = prepend('gpu', gpu_list) # Assembling total list rough_contact_list = \ core_list + generator_list + percolation_list + model_list + solvers_list + bem_list # Adding GPU if needed if env['backend'] == 'cuda': rough_contact_list += gpu_list # Generating dependency files # env.AppendUnique(CXXFLAGS=['-MMD']) # for file_name in rough_contact_list: # obj = file_name.replace('.cpp', '.os') # dep_file_name = file_name.replace('.cpp', '.d') # env.SideEffect(dep_file_name, obj) # env.ParseDepends(dep_file_name) libTamaas = env.SharedLibrary('Tamaas', rough_contact_list) Export('libTamaas') diff --git a/src/model/be_engine.cpp b/src/model/be_engine.cpp index b6bd6e9..1c192a9 100644 --- a/src/model/be_engine.cpp +++ b/src/model/be_engine.cpp @@ -1,85 +1,85 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "be_engine.hh" #include "model.hh" /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ -void solve(IntegralOperator::type kind, - const std::map>& operators, GridBase& input, GridBase& output) { operators.at(kind)->apply(input, output); } -template +template void registerWestergaardOperator( - std::map>& + std::map>& operators, Model& model) { std::stringstream sstr; sstr << "Westergaard::" << otype; if (operators.find(otype) == operators.end()) operators[otype] = model.template registerIntegralOperator>( sstr.str()); } template void BEEngineTmpl::solveNeumann(GridBase& neumann, GridBase& dirichlet) const { solve(IntegralOperator::neumann, this->operators, neumann, dirichlet); } template void BEEngineTmpl::solveDirichlet(GridBase& dirichlet, GridBase& neumann) const { solve(IntegralOperator::dirichlet, this->operators, dirichlet, neumann); } template void BEEngineTmpl::registerNeumann() { registerWestergaardOperator(this->operators, *this->model); } /// Register dirichlet operator template void BEEngineTmpl::registerDirichlet() { registerWestergaardOperator( this->operators, *this->model); } template class BEEngineTmpl; template class BEEngineTmpl; template class BEEngineTmpl; template class BEEngineTmpl; template class BEEngineTmpl; template class BEEngineTmpl; /* -------------------------------------------------------------------------- */ } // namespace tamaas diff --git a/src/model/be_engine.hh b/src/model/be_engine.hh index 3c7c934..067a728 100644 --- a/src/model/be_engine.hh +++ b/src/model/be_engine.hh @@ -1,85 +1,85 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef BE_ENGINE_HH #define BE_ENGINE_HH /* -------------------------------------------------------------------------- */ #include "model_type.hh" #include "integral_operator.hh" #include "westergaard.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /** * @brief Boundary equation engine class. Solves the Neumann/Dirichlet problem * This class should be dimension and model-type agnostic. */ class BEEngine { public: BEEngine(Model* model) : model(model) {} /// Destructor virtual ~BEEngine() = default; public: /// Solve Neumann problem (expects boundary data) virtual void solveNeumann(GridBase& neumann, GridBase& dirichlet) const = 0; /// Solve Dirichlet problem (expects boundary data) virtual void solveDirichlet(GridBase& dirichlet, GridBase& neumann) const = 0; /// Register neumann operator virtual void registerNeumann() = 0; /// Register dirichlet operator virtual void registerDirichlet() = 0; /// Get model const Model& getModel() const { return *model; } protected: Model* model; - std::map> operators; + std::map> operators; }; template class BEEngineTmpl : public BEEngine { public: BEEngineTmpl(Model* model): BEEngine(model) {} void solveNeumann(GridBase& neumann, GridBase& dirichlet) const override; void solveDirichlet(GridBase& dirichlet, GridBase& neumann) const override; /// Register neumann operator void registerNeumann() override; /// Register dirichlet operator void registerDirichlet() override; }; __END_TAMAAS__ #endif // BE_ENGINE_HH diff --git a/src/model/integral_operator.hh b/src/model/integral_operator.hh index f4b9696..4a2a9c9 100644 --- a/src/model/integral_operator.hh +++ b/src/model/integral_operator.hh @@ -1,91 +1,91 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef INTEGRAL_OPERATOR_HH #define INTEGRAL_OPERATOR_HH /* -------------------------------------------------------------------------- */ #include "grid_base.hh" #include #include #include /* -------------------------------------------------------------------------- */ namespace tamaas { /* -------------------------------------------------------------------------- */ class Model; /** * @brief Generic class for integral operators */ class IntegralOperator { public: - /// Type of operator - enum type { neumann, dirichlet }; + /// Kind of operator + enum kind { neumann, dirichlet }; public: /// Constructor IntegralOperator(Model* model) : model(model) {} /// Virtual destructor for subclasses virtual ~IntegralOperator() = default; /// Apply operator on input virtual void apply(GridBase& input, GridBase& output) const = 0; /// Update any data dependent on model parameters virtual void updateFromModel() = 0; /// Get model const Model& getModel() const { return *model; } protected: Model* model = nullptr; }; /* -------------------------------------------------------------------------- */ -/* Printing IntegralOperator::type to string */ +/* Printing IntegralOperator::kind to string */ /* -------------------------------------------------------------------------- */ -#define INTEGRAL_TYPES (neumann)(dirichlet) +#define INTEGRAL_KINDS (neumann)(dirichlet) inline std::ostream& operator<<(std::ostream& o, - const IntegralOperator::type& val) { + const IntegralOperator::kind& val) { switch (val) { -#define PRINT_INTEGRAL_TYPE(r, data, type) \ - case data::type: \ - o << BOOST_PP_STRINGIZE(type); \ +#define PRINT_INTEGRAL_KIND(r, data, kind) \ + case data::kind: \ + o << BOOST_PP_STRINGIZE(kind); \ break; - BOOST_PP_SEQ_FOR_EACH(PRINT_INTEGRAL_TYPE, IntegralOperator, - INTEGRAL_TYPES); -#undef PRINT_INTEGRAL_TYPE + BOOST_PP_SEQ_FOR_EACH(PRINT_INTEGRAL_KIND, IntegralOperator, + INTEGRAL_KINDS); +#undef PRINT_INTEGRAL_KIND } return o; } -#undef INTEGRAL_TYPES +#undef INTEGRAL_KINDS /* -------------------------------------------------------------------------- */ } // namespace tamaas #endif // INTEGRAL_OPERATOR_HH diff --git a/src/model/kelvin.cpp b/src/model/kelvin.cpp index 1d19532..f8369f4 100644 --- a/src/model/kelvin.cpp +++ b/src/model/kelvin.cpp @@ -1,131 +1,103 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "kelvin.hh" #include "elasto_plastic_model.hh" #include "influence.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* Volume 2D implementation */ /* -------------------------------------------------------------------------- */ template <> -void Kelvin::applyVolumeForcePotential( +void Kelvin::apply( GridBase& source, GridBase& out) const { Real nu = model->getPoissonRatio(), mu = model->getShearModulus(); VectorProxy domain(model->getSystemSize()[0]); influence::Kelvin kelvin(mu, nu, domain); auto apply = [this, &kelvin](UInt i, decltype(source_buffers)& source_buffers, decltype(disp_buffer)& displacement) { constexpr UInt dim = trait::dimension; const Real L = this->model->getSystemSize().front(); const UInt N = this->model->getDiscretization().front(); const Real dl = L / N; const Real xi = i * dl; // Compute displacements u_i displacement = 0; for (UInt j : Loop::range(N)) { const Real yj = j * dl; auto& source = source_buffers[j]; #define VOLUME_FORCE_POTENTIAL(yj_xi) \ Loop::stridedLoop( \ [&kelvin, xi, yj, dl](VectorProxy&& u, \ VectorProxy&& f, \ VectorProxy&& q) { \ constexpr UInt dim = trait::dimension; \ Matrix F; \ influence::KelvinIntegrator::integrate(kelvin, yj, xi, \ dl, q, F); \ u += F * f; \ }, \ displacement, source, this->wavevectors) if (j > i) { VOLUME_FORCE_POTENTIAL(1); } else if (j == i) { VOLUME_FORCE_POTENTIAL(0); } else { VOLUME_FORCE_POTENTIAL(-1); } #undef VOLUME_FORCE_POTENTIAL } // Setting fundamental frequency to zero VectorProxy u_fundamental(displacement(0)); u_fundamental = 0; }; this->fourierApply(apply, source, out); } -template <> -void Kelvin::applyVolumeStressPotential( - GridBase& source, GridBase& out) const { - // -} - -template <> -void Kelvin::applyHypersingularStressPotential( - GridBase& source, GridBase& out) const { - // -} - -/* -------------------------------------------------------------------------- */ -/* Volume 1D implementation */ -/* -------------------------------------------------------------------------- */ - -template <> -void Kelvin::applyVolumeForcePotential( - GridBase& source, GridBase& out) const { - // -} - -template <> -void Kelvin::applyVolumeStressPotential( - GridBase& source, GridBase& out) const { - // -} - -template <> -void Kelvin::applyHypersingularStressPotential( - GridBase& source, GridBase& out) const { - // +template +void Kelvin::apply(GridBase& source, + GridBase& out) const { + TAMAAS_EXCEPTION("The requested operator has not been implemented"); } /* -------------------------------------------------------------------------- */ /* Template instanciation */ /* -------------------------------------------------------------------------- */ -template class Kelvin; -template class Kelvin; +template class Kelvin; +template class Kelvin; /* -------------------------------------------------------------------------- */ __END_TAMAAS__ diff --git a/src/model/kelvin.hh b/src/model/kelvin.hh index 1af52f0..3e5902d 100644 --- a/src/model/kelvin.hh +++ b/src/model/kelvin.hh @@ -1,64 +1,56 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __KELVIN_HH__ #define __KELVIN_HH__ /* -------------------------------------------------------------------------- */ -#include "ve_engine.hh" +#include "volume_potential.hh" #include "model_type.hh" #include "grid_hermitian.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /** * @brief Kelvin tensor */ -template -class Kelvin : public VEEngineTmpl { +template +class Kelvin : public VolumePotential { static_assert(type == model_type::volume_1d || type == model_type::volume_2d, "Only volume types are supported"); using trait = model_type_traits; public: /// Constructor - Kelvin(Model * model): VEEngineTmpl(model) {} + Kelvin(Model * model): VolumePotential(model) {} public: - /// Apply the Kelvin tensor to a volume force distribution - void applyVolumeForcePotential(GridBase& source, - GridBase& out) const override; - /// Apply the Kelvin first derivative tensor to an initial stress distribution - void applyVolumeStressPotential(GridBase& source, - GridBase& out) const override; - /// Apply the Kelvin second derivative tensor to an initial stress distribution - void applyHypersingularStressPotential(GridBase& source, - GridBase& out) const override; - + /// Apply the Kelvin-tensor_order operator + void apply(GridBase& source, GridBase& out) const override; }; __END_TAMAAS__ #endif // __KELVIN_HH__ diff --git a/src/model/mindlin.cpp b/src/model/mindlin.cpp deleted file mode 100644 index 5180630..0000000 --- a/src/model/mindlin.cpp +++ /dev/null @@ -1,83 +0,0 @@ -/** - * @file - * - * @author Lucas Frérot - * - * @section LICENSE - * - * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de - * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des - * Solides) - * - * Tamaas is free software: you can redistribute it and/or modify it under the - * terms of the GNU Lesser General Public License as published by the Free - * Software Foundation, either version 3 of the License, or (at your option) any - * later version. - * - * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY - * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR - * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more - * details. - * - * You should have received a copy of the GNU Lesser General Public License - * along with Tamaas. If not, see . - * - */ -/* -------------------------------------------------------------------------- */ -#include "mindlin.hh" -#include "elasto_plastic_model.hh" -/* -------------------------------------------------------------------------- */ -__BEGIN_TAMAAS__ -/* -------------------------------------------------------------------------- */ - -/* -------------------------------------------------------------------------- */ -/* Volume 2D implementation */ -/* -------------------------------------------------------------------------- */ -template <> -void Mindlin::applyVolumeForcePotential( - GridBase& source, GridBase& out) const { - // -} - -template <> -void Mindlin::applyVolumeStressPotential( - GridBase& source, GridBase& out) const { - // -} - -template <> -void Mindlin::applyHypersingularStressPotential( - GridBase& source, GridBase& out) const { - // -} - -/* -------------------------------------------------------------------------- */ -/* Volume 1D implementation */ -/* -------------------------------------------------------------------------- */ - -template <> -void Mindlin::applyVolumeForcePotential( - GridBase& source, GridBase& out) const { - // -} - -template <> -void Mindlin::applyVolumeStressPotential( - GridBase& source, GridBase& out) const { - // -} - -template <> -void Mindlin::applyHypersingularStressPotential( - GridBase& source, GridBase& out) const { - // -} - -/* -------------------------------------------------------------------------- */ -/* Template instanciation */ -/* -------------------------------------------------------------------------- */ -template class Mindlin; -template class Mindlin; - -/* -------------------------------------------------------------------------- */ -__END_TAMAAS__ diff --git a/src/model/mindlin.hh b/src/model/mindlin.hh deleted file mode 100644 index 6a86a94..0000000 --- a/src/model/mindlin.hh +++ /dev/null @@ -1,66 +0,0 @@ -/** - * @file - * - * @author Lucas Frérot - * - * @section LICENSE - * - * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de - * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des - * Solides) - * - * Tamaas is free software: you can redistribute it and/or modify it under the - * terms of the GNU Lesser General Public License as published by the Free - * Software Foundation, either version 3 of the License, or (at your option) any - * later version. - * - * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY - * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR - * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more - * details. - * - * You should have received a copy of the GNU Lesser General Public License - * along with Tamaas. If not, see . - * - */ -/* -------------------------------------------------------------------------- */ -#ifndef __MINDLIN_HH__ -#define __MINDLIN_HH__ -/* -------------------------------------------------------------------------- */ -#include "ve_engine.hh" -#include "model_type.hh" -#include "grid_hermitian.hh" -/* -------------------------------------------------------------------------- */ -__BEGIN_TAMAAS__ - -/** - * @brief Mindlin tensor - */ -template -class Mindlin : public VEEngine { - static_assert(type == model_type::volume_1d || type == model_type::volume_2d, - "Only volume types are supported"); - using trait = model_type_traits; - -public: - /// Constructor - Mindlin(Model * model): VEEngine(model) {} - -public: - /// Apply the Mindlin tensor to a volume force distribution - void applyVolumeForcePotential(GridBase& source, - GridBase& out) const override; - /// Apply the Mindlin first derivative tensor to an initial stress distribution - void applyVolumeStressPotential(GridBase& source, - GridBase& out) const override; - /// Apply the Mindlin second derivative tensor to an initial stress distribution - void applyHypersingularStressPotential(GridBase& source, - GridBase& out) const override; - -protected: - mutable GridHermitian buffer; -}; - - __END_TAMAAS__ - -#endif // __MINDLIN_HH__ diff --git a/src/model/ve_engine.cpp b/src/model/volume_potential.cpp similarity index 92% rename from src/model/ve_engine.cpp rename to src/model/volume_potential.cpp index 3ea1497..3364dc9 100644 --- a/src/model/ve_engine.cpp +++ b/src/model/volume_potential.cpp @@ -1,75 +1,75 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ -#include "ve_engine.hh" +#include "volume_potential.hh" #include "model.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ template -VEEngineTmpl::VEEngineTmpl(Model * model): VEEngine(model) { +VolumePotential::VolumePotential(Model * model): IntegralOperator(model) { // Copy horizontal sizes std::array sizes; std::copy(model->getDiscretization().begin() + 1, model->getDiscretization().end(), sizes.begin()); auto hermitian_sizes = GridHermitian::hermitianDimensions( sizes); wavevectors = FFTransform:: template computeFrequencies(hermitian_sizes); // Normalize wavevectors VectorProxy boundary_domain{ model->getSystemSize()[1]}; wavevectors *= 2 * M_PI; wavevectors /= boundary_domain; /// Initializing buffers auto initialize = [this, &hermitian_sizes](std::vector& buffers) { buffers.resize(this->model->getDiscretization()[0]); std::for_each(buffers.begin(), buffers.end(), [&hermitian_sizes](BufferType& buffer) { buffer.setNbComponents(trait::components); buffer.resize(hermitian_sizes); }); }; initialize(source_buffers); disp_buffer.setNbComponents(trait::components); disp_buffer.resize(hermitian_sizes); } /* ------------------------------------------------------------------------ */ /* Template instanciation */ /* ------------------------------------------------------------------------ */ -template class VEEngineTmpl; -template class VEEngineTmpl; +template class VolumePotential; +template class VolumePotential; __END_TAMAAS__ diff --git a/src/model/ve_engine.hh b/src/model/volume_potential.hh similarity index 67% rename from src/model/ve_engine.hh rename to src/model/volume_potential.hh index caa2ed1..3f85f1a 100644 --- a/src/model/ve_engine.hh +++ b/src/model/volume_potential.hh @@ -1,122 +1,95 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ -#ifndef __VE_ENGINE_HH__ -#define __VE_ENGINE_HH__ +#ifndef VOLUME_POTENTIAL_HH +#define VOLUME_POTENTIAL_HH /* -------------------------------------------------------------------------- */ -#include "tamaas.hh" -#include "grid_base.hh" -#include "grid_view.hh" -#include "grid_hermitian.hh" +#include "integral_operator.hh" #include "model_type.hh" +#include "grid_hermitian.hh" #include "fft_plan_manager.hh" -#include "model.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /** - * @brief Volume equation engine class. Applies the operators for computation + * @brief Volume potential operator class. Applies the operators for computation * of displacements and strains due to residual/eigen strains */ -class VEEngine { -public: - /// Constructor - VEEngine(Model* model) : model(model) {} - /// Destructor - virtual ~VEEngine() = default; - - // Methods to override -public: - /// Compute displacement field in volume to to a volume force distribution - virtual void applyVolumeForcePotential(GridBase& source, - GridBase& out) const = 0; - /// Compute displacement field in volume due to an initial stress distribution - virtual void applyVolumeStressPotential(GridBase& source, - GridBase& out) const = 0; - /// Compute strain field in volume due to an initial stress distribution - virtual void applyHypersingularStressPotential(GridBase& source, - GridBase& out) const = 0; - -protected: - Model* model; -}; - -/** - * @brief Templated class to manage application of functor on fields in Fourier - * space. - */ template -class VEEngineTmpl : public VEEngine { +class VolumePotential : public IntegralOperator { using trait = model_type_traits; public: - VEEngineTmpl(Model * model); + VolumePotential(Model * model); + + /// Update from model (does nothing) + void updateFromModel() override {} protected: /// Function to handle layer-by-layer Fourier treatment of operators template void fourierApply(Func func, GridBase& in, GridBase& out) const; protected: Grid wavevectors; using BufferType = GridHermitian; mutable std::vector source_buffers; mutable BufferType disp_buffer; }; /* -------------------------------------------------------------------------- */ /* Template implementation */ /* -------------------------------------------------------------------------- */ template template -void VEEngineTmpl::fourierApply(Func func, GridBase& in, - GridBase& out) const { +void VolumePotential::fourierApply(Func func, GridBase& in, + GridBase& out) const { constexpr UInt dim = trait::dimension; Grid& i = dynamic_cast(in); Grid& o = dynamic_cast(out); TAMAAS_ASSERT(i.sizes().front() == o.sizes().front(), "Number of layers does not match"); for (UInt layer : Loop::range(i.sizes().front())) { auto in_layer = make_view(i, layer); FFTPlanManager::get() .createPlan(in_layer, source_buffers[layer]) .forwardTransform(); } for (UInt layer : Loop::range(i.sizes().front())) { func(layer, source_buffers, disp_buffer); auto out_layer = make_view(o, layer); FFTPlanManager::get() .createPlan(out_layer, disp_buffer) .backwardTransform(); } } __END_TAMAAS__ -#endif // __VE_ENGINE_HH__ +#endif // VOLUME_POTENTIAL_HH diff --git a/src/model/westergaard.cpp b/src/model/westergaard.cpp index ec2cc26..8d58549 100644 --- a/src/model/westergaard.cpp +++ b/src/model/westergaard.cpp @@ -1,215 +1,215 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "fft_plan_manager.hh" #include "influence.hh" #include "loop.hh" #include "model.hh" #include "grid_view.hh" #include "static_types.hh" #include "westergaard.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /* -------------------------------------------------------------------------- */ -template +template Westergaard::Westergaard(Model* model) : IntegralOperator(model) { // Copy sizes std::array sizes; std::copy(model->getDiscretization().begin(), model->getDiscretization().end(), sizes.begin()); // Copy boundary sizes auto sizes_it = sizes.begin(); std::array boundary_hermitian_sizes; // Ignoring first dimension if model is volumetric if (trait::dimension > trait::boundary_dimension) ++sizes_it; std::copy(sizes_it, sizes.end(), boundary_hermitian_sizes.begin()); boundary_hermitian_sizes = GridHermitian::hermitianDimensions( boundary_hermitian_sizes); constexpr UInt nb_components = trait::components; buffer.setNbComponents(nb_components); buffer.resize(boundary_hermitian_sizes); influence.setNbComponents(nb_components * nb_components); influence.resize(boundary_hermitian_sizes); initInfluence(); } /* -------------------------------------------------------------------------- */ -template +template template void Westergaard::fourierApply(Functor func, GridBase& in, GridBase& out) const try { Grid& i = dynamic_cast(in); Grid& full_out = dynamic_cast(out); FFTPlanManager::get().createPlan(i, buffer).forwardTransform(); /// Applying influence func(buffer, influence); /// Backward fourier transform on boundary only if constexpr (bdim == dim) FFTPlanManager::get().createPlan(full_out, buffer).backwardTransform(); else { auto view = make_view(full_out, 0); FFTPlanManager::get() .createPlan(view, buffer) .backwardTransform(); } } catch (const std::bad_cast& e) { TAMAAS_EXCEPTION("Neumann and dirichlet types do not match model type"); } /* -------------------------------------------------------------------------- */ -template +template template void Westergaard::initFromFunctor(Functor func) { // Compute wavevectors for influence auto wavevectors = FFTransform::template computeFrequencies( influence.sizes()); // Get boundary physical size VectorProxy domain( this->model->getSystemSize()[(bdim == dim) ? 0 : 1]); // Normalize wavevectors wavevectors *= 2 * M_PI; wavevectors /= domain; // Apply functor Loop::stridedLoop(func, wavevectors, influence); // Set fundamental mode to zero MatrixProxy mat(influence(0)); mat = 0; } /* -------------------------------------------------------------------------- */ #define NEUMANN_BASIC(type) \ template <> \ void Westergaard::initInfluence() { \ auto E_star = model->getHertzModulus(); \ auto basic = [E_star] CUDA_LAMBDA(VectorProxy && q, \ Complex & k) { \ k = 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, \ Complex & k) { \ k = E_star * q.l2norm() / 2; \ }; \ initFromFunctor(basic); \ } NEUMANN_BASIC(model_type::basic_1d); NEUMANN_BASIC(model_type::basic_2d); DIRICHLET_BASIC(model_type::basic_1d); DIRICHLET_BASIC(model_type::basic_2d); #undef NEUMANN_BASIC #undef DIRICHLET_BASIC /* -------------------------------------------------------------------------- */ template <> void Westergaard::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_norm = q.l2norm(); q /= q_norm; F(0, 0) = 2 * (1 + nu) * (1 - nu * q(0) * q(0)); F(1, 1) = 2 * (1 + nu) * (1 - nu * q(1) * q(1)); F(2, 2) = 2 * (1 - nu * nu); F(0, 1) = F(1, 0) = -q(0) * q(1) * 2 * nu * (1 + nu); F(0, 2) = I * q(0) * (1 + nu) * (1. - 2. * nu); F(1, 2) = I * q(1) * (1 + nu) * (1 - 2 * nu); F(2, 0) = -F(0, 2); F(2, 1) = -F(1, 2); F *= 1. / (E * q_norm); }; initFromFunctor(surface); } /* -------------------------------------------------------------------------- */ -template +template void Westergaard::initInfluence() { TAMAAS_EXCEPTION("the requested operator has not been implemented"); } /* ------------------------------------------------------------------------ */ -template +template void Westergaard::apply(GridBase& input, GridBase& output) const { auto apply = [](decltype(buffer)& buffer, const decltype(influence)& influence) { Loop::stridedLoop([] CUDA_LAMBDA(VectorProxy && i, MatrixProxy && F) { i = F * i; }, buffer, influence); }; fourierApply(apply, input, output); } /* -------------------------------------------------------------------------- */ /* 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; /* -------------------------------------------------------------------------- */ __END_TAMAAS__ /* -------------------------------------------------------------------------- */ diff --git a/src/model/westergaard.hh b/src/model/westergaard.hh index b3b1905..acd2f85 100644 --- a/src/model/westergaard.hh +++ b/src/model/westergaard.hh @@ -1,87 +1,87 @@ /** * @file * * @author Lucas Frérot * * @section LICENSE * * Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de * Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des * Solides) * * Tamaas is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Tamaas. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __WESTERGAARD_HH__ #define __WESTERGAARD_HH__ /* -------------------------------------------------------------------------- */ #include "be_engine.hh" #include "fft_plan_manager.hh" #include "grid_hermitian.hh" #include "integral_operator.hh" #include "model_type.hh" #include "tamaas.hh" /* -------------------------------------------------------------------------- */ __BEGIN_TAMAAS__ /** * @brief Operator based on Westergaard solution and the Dicrete Fourier * Transform. * This class is templated with model type to allow efficient storage of the * influence coefficients. * The integral operator is only applied to surface pressure/displacements, * even for volume models. */ -template +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& ouput) const override; /// Update the influence coefficients void updateFromModel() override { initInfluence(); } private: /// 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; public: GridHermitian influence; mutable GridHermitian buffer; }; __END_TAMAAS__ #endif // __WESTERGAARD_HH__ diff --git a/tests/test_integral_operators.py b/tests/test_integral_operators.py index 613aac2..a64baee 100644 --- a/tests/test_integral_operators.py +++ b/tests/test_integral_operators.py @@ -1,68 +1,68 @@ #!/usr/bin/env python # coding: utf-8 # ----------------------------------------------------------------------------- # @author Lucas Frérot # # @section LICENSE # # Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de # Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des # Solides) # # Tamaas is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Tamaas is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Tamaas. If not, see . # ----------------------------------------------------------------------------- import tamaas as tm import numpy as np from numpy.linalg import norm def test_kelvin_volume_force(): N = 65 E = 1. nu = 0.3 mu = E / (2*(1+nu)) domain = np.array([1.] * 3) omega = 2 * np.pi * np.array([1, 1]) / domain[:2] omega_ = norm(omega) discretization = [N] * 3 model = tm.ModelFactory.createModel(tm.model_type.volume_2d, domain, discretization) model.E = E model.nu = nu engine = tm._tamaas._test_features.Kelvin(model) coords = [np.linspace(0, domain[i], discretization[i], endpoint=False) for i in range(2)]\ +[np.linspace(0, domain[2], discretization[2])] x, y = np.meshgrid(*coords[:2], indexing='ij') displacement = model.getDisplacement() source = np.zeros_like(displacement) source[N//2, :, :, 2] = np.sin(omega[0]*x) * np.sin(omega[1]*y) - engine.applyVolumeForcePotential(source, displacement) + engine.apply(source, displacement) z = coords[2] - 0.5 z, x, y = np.meshgrid(z, *coords[:2], indexing='ij') solution = np.zeros_like(source) solution[:, :, :, 0] = np.exp(-omega_*np.abs(z)) / (8*mu*(1-nu)*omega_) * omega[0]*z*np.cos(omega[0]*x)*np.sin(omega[1]*y) solution[:, :, :, 1] = np.exp(-omega_*np.abs(z)) / (8*mu*(1-nu)*omega_) * omega[1]*z*np.sin(omega[0]*x)*np.cos(omega[1]*y) solution[:, :, :, 2] = np.exp(-omega_*np.abs(z)) / (8*mu*(1-nu)*omega_) * (3-4*nu + omega_*np.abs(z))*np.sin(omega[0]*x)*np.sin(omega[1]*y) error = norm(displacement - solution) / norm(solution) assert error < 5e-2